Jacobians for arrays of static arrays

[ChrisRackauckas]

using StaticArrays, OrdinaryDiffEq

@inline @inbounds function loop(u, p, t)
    σ = p[1]; ρ = p[2]; β = p[3]
    du1 = σ*(u[2]-u[1])
    du2 = u[1]*(ρ-u[3]) - u[2]
    du3 = u[1]*u[2] - β*u[3]
    return SVector{3}(du1, du2, du3)
end

p = [10, 28, 8/3]

state = [0, 10.0, 0]

states = [state, state + 1e-9rand(3)]

st = [SVector{3}(s) for s in states]
L = length(st)
# The following may be inneficient
paralleleom = (du, u, p, t) -> begin
    for i in 1:L
        @inbounds du[i] = loop(u[i], p, t)
    end
end


pprob = ODEProblem(paralleleom, st, (0.0, Inf), p)

init(pprob, Vern9(); save_everystep = false)
init(pprob, Tsit5(); save_everystep = false)
init(pprob, KenCarp5(); save_everystep = false)

This might just need generic linalg tools to really be handled appropriately.

[ChrisRackauckas]

Using reinterpret can be an option to solve this.