Approximate equality in not-in-place stepping

[ChrisRackauckas]

This commit https://github.com/JuliaDiffEq/DelayDiffEq.jl/pull/72/commits/d3c5baeaf35e94433e3a5fd8075d0d3c15067f8b 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]

https://github.com/JuliaDiffEq/DelayDiffEq.jl/pull/72/commits/d0a7c60598a1634332978143c05c9e3c157ac01d

[ChrisRackauckas]

https://github.com/JuliaDiffEq/DelayDiffEq.jl/pull/72/commits/8c8885d30f2a92875608252ea8451932bf681364

[ChrisRackauckas]

Could potentially have been fixed by https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/commit/b86aae7eb0cddaa1e8fb76c75a67eac366e26ebb

[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.