Approximate equality in not-in-place stepping

[ChrisRackauckas]

This commit d3c5bae turned exact checks to approximate checks because the inplace vs outofplace versions now only step within floating point error. This is a minor point would could be corrected again in the future since on v0.6 it was completely error free.

[ChrisRackauckas]

d0a7c60

[ChrisRackauckas]

8c8885d

[ChrisRackauckas]

Could potentially have been fixed by SciML/OrdinaryDiffEq.jl@b86aae7

[devmotion]

I observed that (at least) one approximate check is required due to the fact that interpolating an ODE solution returns different values for no, scalar, and vector indices in some cases (on Julia 1.0.2, current master branch of DiffEqBase, OrdinaryDiffEq, and DelayDiffEq):

using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary
DDEProblemLibrary.importddeproblems()

sol = solve(prob_dde_1delay, MethodOfSteps(Tsit5()))

sol(7) # -0.05731721456328998
sol(7; idxs=1) # -0.05731721456328997
sol(7; idxs=[1])[1] # -0.05731721456328998

Hence (some of) the approximate checks might be caused by inconsistent interpolations in OrdinaryDiffEq, I guess. In that case it should be possible to observe the same approximate equalities when solving ODEs (which I haven't tried).

[devmotion]

After removing @muladd in https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/dense/interpolants.jl#L449, https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/dense/interpolants.jl#L457, and https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/dense/interpolants.jl#L466 I obtain the same result -0.05731721456340372 in all three use cases. However, I don't understand why adding @muladd leads to inconsistent results since

julia> @macroexpand @muladd @inbounds for (j,i) in enumerate(idxs)
                  out[j] = y₀[i] + dt*(k[1][i]*b1Θ + k[2][i]*b2Θ + k[3][i]*b3Θ + k[4][i]*b4Θ + k[5][i]*b5Θ + k[6][i]*b6Θ + k[7][i]*b7Θ)
                end
quote
    $(Expr(:inbounds, true))
    local #1554#val = for (j, i) = enumerate(idxs)
                #= REPL[79]:2 =#
                out[j] = (muladd)(dt, (muladd)((k[7])[i], b7Θ, (muladd)((k[6])[i], b6Θ, (muladd)((k[5])[i], b5Θ, (muladd)((k[4])[i], b4Θ, (muladd)((k[3])[i], b3Θ, (muladd)((k[2])[i], b2Θ, (k[1])[i] * b1Θ)))))), y₀[i])
            end
    $(Expr(:inbounds, :pop))
    #1554#val
end

and

julia> @macroexpand @muladd y₀[idxs] + dt*(k[1][idxs]*b1Θ + k[2][idxs]*b2Θ + k[3][idxs]*b3Θ +
                                      k[4][idxs]*b4Θ + k[5][idxs]*b5Θ + k[6][idxs]*b6Θ + k[7][idxs]*b7Θ)
:((muladd)(dt, (muladd)((k[7])[idxs], b7Θ, (muladd)((k[6])[idxs], b6Θ, (muladd)((k[5])[idxs], b5Θ, (muladd)((k[4])[idxs], b4Θ, (muladd)((k[3])[idxs], b3Θ, (muladd)((k[2])[idxs], b2Θ, (k[1])[idxs] * b1Θ)))))), y₀[idxs]))

indicate that muladd is applied in the same way.