Add C RHS using ccall

[sdwfrost]

The idea here is to use a vector field defined in C. The following works:

using OrdinaryDiffEq
using Libdl

C_code = """
void sir_ode(double *du, double *u, double *p, double *t){
    double β = p[0];
    double c = p[1];
    double γ = p[2];
    double S = u[0];
    double I = u[1];
    double R = u[2];
    double N = S + I + R;
    du[0] = -β*c*S*I/N;
    du[1] = β*c*S*I/N - γ*I;
    du[2] = γ*I;
}
""";

const Clib = tempname();

open(`gcc -fPIC -O3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code)
end

function sir_ode_jl(du,u,p,t)
    ccall((:sir_ode,Clib,), Cvoid,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}), du, u, p, Ref(t))
end

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
u0 = [990.0,10.0,0.0] # S,I,R
p = [0.05,10.0,0.25] # β,c,γ
prob = ODEProblem{true}(sir_ode_jl, u0, tspan, p)
sol = solve(prob, Tsit5(), dt = δt)

An out-of-place version is probably too fiddly given that in-place is preferable in any case.

[sdwfrost]

Closed with 8045d5c020b47bbb8443622edff787f749cf1f92