Different solutions for in-place and not in-place functions with Vern7

[devmotion]

During testing of PR SciML/DelayDiffEq.jl#17 I discovered that Vern7 produces slightly different results for in-place and not in-place functions already for ODEs:

julia> using OrdinaryDiffEq, DiffEqProblemLibrary

julia> prob_inplace = ODEProblem(DiffEqProblemLibrary.f_2dlinear, [0.5, 1.0], (0., 10.))

julia> prob_notinplace = ODEProblem(DiffEqProblemLibrary.f_2dlinear_notinplace, [0.5, 1.0], (0., 10.))

julia> sol_inplace = solve(prob_inplace, Vern7())

julia> sol_notinplace = solve(prob_notinplace, Vern7())

julia> sol_inplace.t == sol_notinplace.t
false

julia> sol_inplace.t[3]
0.6807461486962336

julia> sol_notinplace.t[3]
0.6807461622279954

julia> sol_inplace.errors[:l2]
0.006139913326522092

julia> sol_notinplace.errors[:l2]
0.006139913425004702

[devmotion]

I found the error, it's the first @muladd in line https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/verner_rk_integrators.jl#L320, or rather the fact that there's a subtraction inside the expression following @muladd but not in the corresponding line of the not in-place method https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/verner_rk_integrators.jl#L218. Moving @muladd to the same position as in the not in-place method fixes the problem, just removing it still leads to a small difference between both methods. So apparently one has to position @muladd at exactly the same place for corresponding methods to get the same results.

Interestingly, for Vern8 and Vern9 the expressions following @muladd in the lines in which EEst is calculated always contain the subtraction as well (Vern6 is handled differently). So should one rather change line https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/verner_rk_integrators.jl#L218 for consistency? Or change it in those lines, too?

[ChrisRackauckas]

The best thing to check would be

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/verner_rk_integrators.jl#L320

macroexpand that: does it create the correct expression? If it does then it's slightly more correct since using muladd reduces numerical error. It must because Vern7 works pretty well? But I have cases where it does expand to something subtly wrong. It always works on linear combinations well though, so we can always fallback to that.

[devmotion]

I checked it with macroexpand. Currently it creates the expression

julia> macroexpand(:(@inbounds atmp[i] = (@muladd(update[i] - dt*(bhat1*k1[i] + bhat4*k4[i] + bhat5*k5[i] + bhat6*k6[i] + bhat7*k7[i] + bhat10*k10[i]))./@muladd(atol+max(abs(uprev[i]),abs(u[i])).*rtol))))
quote 
    $(Expr(:inbounds, true))
    atmp[i] = (update[i] - dt * (bhat1 * k1[i] + bhat4 * k4[i] + bhat5 * k5[i] + bhat6 * k6[i] + bhat7 * k7[i] + bhat10 * k10[i])) ./ (atol + max(abs(uprev[i]), abs(u[i])) .* rtol)
    $(Expr(:inbounds, :pop))
end

whereas moving the first @muladd and removing the last unnecessary dot leads to

julia> macroexpand(:(@inbounds atmp[i] = (update[i] - dt*(@muladd(bhat1*k1[i] + bhat4*k4[i] + bhat5*k5[i] + bhat6*k6[i] + bhat7*k7[i] + bhat10*k10[i])))./@muladd(atol+max(abs(uprev[i]),abs(u[i]))*rtol)))
quote 
    $(Expr(:inbounds, true))
    atmp[i] = (update[i] - dt * (muladd)(bhat1, k1[i], (muladd)(bhat4, k4[i], (muladd)(bhat5, k5[i], (muladd)(bhat6, k6[i], (muladd)(bhat7, k7[i], bhat10 * k10[i])))))) ./ (muladd)(max(abs(uprev[i]), abs(u[i])), rtol, atol)
    $(Expr(:inbounds, :pop))
end

So I guess we should change it to the latter expression. However, since the corresponding line of the not in-place method creates the expression

julia> macroexpand(:(integrator.EEst = integrator.opts.internalnorm( ((update - dt*(@muladd(bhat1*k1 + bhat4*k4 + bhat5*k5 + bhat6*k6 + bhat7*k7 + bhat10*k10)))./@muladd(integrator.opts.abstol+max.(abs.(uprev),abs.(u)).*integrator.opts.reltol)))))
:(integrator.EEst = integrator.opts.internalnorm((update - dt * (muladd)(bhat1, k1, (muladd)(bhat4, k4, (muladd)(bhat5, k5, (muladd)(bhat6, k6, (muladd)(bhat7, k7, bhat10 * k10)))))) ./ (integrator.opts.abstol + max.(abs.(uprev), abs.(u)) .* integrator.opts.reltol)))

and hence does not apply muladd to the last part because of the necessary dot, this might again lead to different results in some cases. However, it produces the same results for my toy example but I do not know whether either the second method needs to be changed to a for loop (muladd does not work with vector multiplication) or the second @muladd should be removed in both cases to be completely sure.

By the way, the same problems occur for Vern8 and Vern9 in lines

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/verner_rk_integrators.jl#L393
https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/verner_rk_integrators.jl#L515
https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/verner_rk_integrators.jl#L595
https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/verner_rk_integrators.jl#L736

There are two @muladd in each of those lines, but there are no calls of muladd in any of the created expressions. It just consistently does not work, and hence it could not be discovered by any of my comparisons. Probably the problem that because of the dots the second @muladd in each of those lines does not create a call to muladd does occur for some other methods as well (at least, I remember some other integrators with the same expression).

[ChrisRackauckas]

I thought I fixed the @muladd macro for this before, but it didn't carry over to v0.6 because it looked for .+ and .* which no longer exist. I'm not sure how to handle this when broadcast fusion can occur. But this isn't a correctness issue per se, it's very small performance issue. We should try to get this cleaned up as much as possible but it's not breaking at all.

[devmotion]

So do you suggest removing the second @muladd in cases where broadcast fusion can occur, such as

@muladd(integrator.opts.abstol+max.(abs.(uprev),abs.(u)).*integrator.opts.reltol)))

? And removing dots in lines with only scalar operations, such as

@muladd(atol+max(abs(uprev[i]),abs(u[i])).*rtol))

? And regarding the original problem with subtractions (which is caused by

https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/utils.jl#L34

) - should @muladd just be moved as shown above or do you rather want to change the definition of to_muladd such that it can handle substractions as well by not changing subtractions and just adding calls to muladd at other parts of the expression?

[ChrisRackauckas]

So do you suggest removing the second @muladd in cases where broadcast fusion can occur, such as

Well, we don't want to remove @muladd in the equivalent inplace version places since using muladd is numerically strictly better than not. I think here we may just want to find/replace in the more complicated version:

muladd.(integrator.opts.abstol,max.(abs.(uprev),abs.(u)),integrator.opts.reltol)

That should then fuse correctly and be everything we need.

And removing dots in lines with only scalar operations, such as

Yeah. That'll make nested arrays never work though, but I am not convinced that it's something we could support without a VectorOfArray wrapper in the future. Reference test which shows a failure in this area already:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/test/static_array_tests.jl#L9

So sure, delete these extra dots is a good fix.

should @muladd just be moved as shown above or do you rather want to change the definition of to_muladd such that it can handle substractions as well by not changing subtractions and just adding calls to muladd at other parts of the expression?

It's probably best to just extend it to handle subtraction, though @muladd may be getting replaced by @fastmath down the line (though there are some issues with that).

[ChrisRackauckas]

Since @muladd is applied to the full functions by #89, are we good now?

[devmotion]

Just tested the example above on the master branch, and it seems to work now (although I still have to update the Vern methods to the syntax of the new @muladd).