Hi,
I was writing a quadrature code with inplace integrands and I noticed that Integrals.jl has the same problem I noticed. Depending on how the quadrature library allocates arrays, it is possible that multiple integrand evaluations are made before they are summed together. When the result is overwritten to the same array this means some of the quadrature nodes are overwritten by the values of another. For example, integrating sin(x) on [-1,1] as below should give zero on the initial segment, but the inplace evaluator seems to integrate sin(abs(x)) on that segment.
julia> using Integrals
julia> prob = IntegralProblem((x,p) -> sin(p*x), -1, 1, 3.0)
IntegralProblem. In-place: false
julia> inplaceprob = IntegralProblem((y,x,p) -> y .= sin(p*x), -1.0, 1.0, 3.0)
IntegralProblem. In-place: true
julia> solve(prob, HCubatureJL(), maxiters=15)
u: 0.0
julia> solve(inplaceprob, HCubatureJL(), maxiters=15)
u: 1-element Vector{Float64}:
1.3044979103237169The problem is in the definition of the integrand wrapper
dx = zeros(eltype(lb), prob.nout)
f = x -> (prob.f(dx, x, prob.p); dx)
An issue with dx is that it may have the wrong type, since the user's function could return complex values at real nodes, for example. A fix is to allow the user to pass in the output array. The second issue in the example above can be fixed by returning a new array at each integrand evaluation, which can be done by returning one(eltype(lb))*dx. The downside to this last fix is that most of the array operations aren't inplace, but the only way to achieve this is if the library supports them (e.g. see quadgk! in QuadGK.jl for an implementation).
Hi,
I was writing a quadrature code with inplace integrands and I noticed that Integrals.jl has the same problem I noticed. Depending on how the quadrature library allocates arrays, it is possible that multiple integrand evaluations are made before they are summed together. When the result is overwritten to the same array this means some of the quadrature nodes are overwritten by the values of another. For example, integrating
sin(x)on[-1,1]as below should give zero on the initial segment, but the inplace evaluator seems to integratesin(abs(x))on that segment.The problem is in the definition of the integrand wrapper
An issue with
dxis that it may have the wrong type, since the user's function could return complex values at real nodes, for example. A fix is to allow the user to pass in the output array. The second issue in the example above can be fixed by returning a new array at each integrand evaluation, which can be done by returningone(eltype(lb))*dx. The downside to this last fix is that most of the array operations aren't inplace, but the only way to achieve this is if the library supports them (e.g. seequadgk!in QuadGK.jl for an implementation).