Compilation does not terminate (>3 hours) with nested numerical integration

I am trying to do two-dimensional numerical integration. For this my first attempt was something like the following, using DoubleExponentialFormulas.jl. This example is significantly reduced and no longer does anything mathematically useful:

using DoubleExponentialFormulas
#using QuadGK

function f4(c, t, S, α)
    t^(c - 1) * S - (c - 1) * α
end

function f3(j, t, S, α)
    S * prod((f4(c, t, S, α) / c) for c = 0:j)
end

function f2(n, a, b, r, ξ, η)
    sum(f3(j, (r - b) / a, r, b / a) for j = 0:(n+2))
end

function inner_integral(n, a, b, δ, ξ)
    f(η) = f2(n, a, b, a, ξ, η)
    I, E = quadde(f, 0, 2π, rtol = δ)
    I
end

function outer_integral(n, a, b, δ)
    f(ξ) = inner_integral(n, a, b, δ, ξ) * cos(b + ξ)
    I, E = quadde(f, -π / 2 - b, 0, π / 2 - b, rtol = δ)
    I
end

function sum_of_integrals(N, a, b, δ)
    sum(outer_integral(n, a, b, δ) for n = 0:N)
end

@time println(sum_of_integrals(3, 100, 0.2, 10^-6))

However, I cannot get this code to compile within a reasonable timeframe. I left a variant of this running for 3.5 hours before I ran out of system memory. I have verified that this time is actually spent in (some early step of) compilation: inserting a statement that should lead to a compilation error, e.g. a call to a non-existent function, does not cause an error during all this time.

Using QuadGK.jl instead (by switching out the using statements and replacing the quadde calls with quadgk) makes the compilation work but it nevertheless takes around 90 seconds, which still seems extremely slow to me. (Regardless, for mathematical reasons I have to use double-exponential quadrature for integration, so I cannot use quadgk) There are some other things that speed up the compilation some more, most importantly inlining the f3 and f4 functions as loops within f2 or making the integration intervals smaller, but it never reaches a usable point with quadde for me.

I am running Julia 1.6.1 on Arch Linux; my CPU is an i5-8350U.

I am aware that nested 1D integration is not the best way to do multidimensional integration, but at the same time I feel the compiler should not behave in this way. In a sense, I am not really looking for a solution to make the compilation run faster (the real solution would probably be “use an actual multidimensional integrator”), but more asking whether this is expected behaviour from the compiler, and if not, whether it would make sense to investigate this further on my part so that the compiler may be improved in the future.

3 Likes