Many of the implicit SDE solvers fail when adaptive=true with no method matching +(::Array{Float64,1}, ::Float64)

[Krastanov]

f(u,p,t) = 0.4.*u
g(u,p,t) = 0.0.*u
f_jac(u,p,t) = [0.4 0;0 0.4]
u0 = [1., 1.];
sde = SDEFunction(f,g,jac=f_jac)
prob = SDEProblem(sde,g,u0,(0.0,1.0),noise=W);
solve(prob, SKenCarp(), saveat=dt_save)

fails with:

MethodError: no method matching +(::Array{Float64,1}, ::Float64)

Stacktrace:
 [1] muladd(::Array{Float64,1}, ::Float64, ::Float64) at ./math.jl:997
 [2] perform_step! ([...]) at /home/stefan/.julia/packages/StochasticDiffEq/zInIS/src/perform_step/sdirk.jl:77
 [3] perform_step! at /home/stefan/.julia/packages/StochasticDiffEq/zInIS/src/perform_step/sdirk.jl:6 [inlined]
 [4] macro expansion at /home/stefan/.julia/packages/StochasticDiffEq/zInIS/src/integrators/integrator_utils.jl:66 [inlined]
 [5] solve!

The failure happens with other implicit solvers as well (presumably, all the solvers from sdirk.jl).

If I use adaptive=false it works.

[Krastanov]

Take what I say with a large grain of salt, but I think the issue is here:

https://github.com/JuliaDiffEq/StochasticDiffEq.jl/blob/master/src/perform_step/sdirk.jl#L80

In the if branch En is an array, while in the else branch it is a scalar... and then it is added to an array...

[devmotion]

I'm pretty sure that the error is caused by https://github.com/JuliaDiffEq/StochasticDiffEq.jl/blob/master/src/perform_step/sdirk.jl#L88. All operations (apart from the outer internalnorm) should be dotted: the @muladd macro rewrites abstol + something * reltol as muladd(something, reltol, abstol), which fails since in your case abstol and reltol are scalars whereas something is an array. This also explains your error message.

Actually, the same problem occurs with other algorithms as well (e.g. a similar error that you mentioned in your issue in DifferentialEquations.jl is caused by the same problem), hence I was about to propose changing the computations of residuals in a similar way as it is done in OrdinaryDiffEq (as I mentioned before in your last PR).

[ChrisRackauckas]

Yes, this is due to changes in Julia v1.0. In v0.6, [1,2,3] + 1 worked. Now it doesn't. That hurts a lot of out of place usage where we didn't consciously add some dots, so we need more tests for that. Generally, the answer for your problem here is that, since you are using Array types, you should be using the in-place methods (du,u,p,t) and that would fix these issues. However, we should get out of place working as well.

There is a confounding factor that with the v1.0 upgrade, I found SciML/OrdinaryDiffEq.jl#571 and JuliaLang/julia#28126 and so I ended up not taking the time to broadcast everything correctly in the updated versions because doing so would lead to some major performance regressions on code that is not designed for arrays in the first place. This kind of put us in a bind because I want it to work anyways since out of place with arrays does have its uses, including simplicity, but I don't want the SVector usage to be hit hard when arrays should be using in-place.

However, we are tracking this JuliaArrays/StaticArrays.jl#560 and hoping to get it fixed directly.