JuliaDiffEq Vern9() solver taking too long

General Usage
1

I wrote a simple julia script to benchmark Verner9 solver from JuliaDiffEq against my fortran code and I am seeing the former is taking too long compared to the latter. Since I am very new to Julia, i thought someone here can help me figuring out why my julia script is unusually slow for this solver. The julia code (given below) takes 102.2 secs whereas my fortran code on the same computer does not take more than 1.5 secs.

using DifferentialEquations
using LinearAlgebra

const runs = 1000
const mu = 0.012277471

 struct Solver
    solver
    dense::Bool
    name::String
end

solvers = [Solver(Vern9(),
                  false,
                  "Verner9")]

function fire_ode(x0, ode, tspan, atol, rtol, method, prob)
    solve(prob, method, abstol=atol, reltol=rtol, timeseries_errors=false, dense_errors=false, dense=false)
end

struct ODETestResult
    final_state_error
    iom_error
    total_time
    ode_sol
end

function fire_ode_test(x0, ode, tspan, atol, rtol, method, final_state, iom, iom_func)
    prob = ODEProblem(ode, x0, tspan)
    local sol
    rtime = @elapsed for i in 1:runs
        sol = fire_ode(x0, ode, tspan, atol, rtol, method, prob)
    end
    serr = norm(sol[:, end] - final_state)
    iomsol = iom_func(sol)
    ierr = abs.(iomsol .- iom)
    ODETestResult(serr, ierr, rtime, sol)
end

function cr3bpeom(dX, X, p, t)
    @inbounds begin
        r1 = sqrt((X[1] + mu)^2 + X[2]^2)
        r2 = sqrt((X[1] - 1 + mu)^2 + X[2]^2)
        dX[1:3] .= @view X[4:6]
        dX[4] = X[1] + 2*X[5] - (1 - mu)*(X[1] + mu)/r1^3 - mu*(X[1] - 1 + mu)/r2^3
        dX[5] = X[2] - 2*X[4] - (1 - mu)*X[2]/r1^3 - mu*X[2]/r2^3
        dX[6] = 0.0
    end
end

R0 = [0.994, 0.0, 0.0]
V0 = [0.0, -2.00158510637908252240537862224, 0.0]
X0 = [R0; V0]

# Time period
Period = 17.0652165601579625588917206249

# propagate for bunch of orbits
t0 = 0.0
tf = 4.0*Period
tspan = (t0, tf)

for solver in solvers
    result = fire_ode_test(X0, cr3bpeom, tspan, atol, rtol, solver.solver, [R0; V0], 0.0, jacobi_constant)
    @show result.total_time
end
2

I think you may be measuring the compilation time of the julia code as well. There is a package called BechmarkTools.jl which runs code multiple times to get more accurate data for a benchmark.

You are summing up the time for 1000 runs(=compilation time + 1000 runs). Is fortran doing the same?

3

I am sure I am not including the compilation time. I run the script 2 or 3 times and the time I quoted is pretty consistent over all the runs.

I actually know the package BenchmarkTools.jl and it is doing (almost) the same thing that I am doing in this simple script. It is easy to isolate problems in simple scripts.

Of course time taken by my Fortran code includes the same number of runs and same number of orbit time-periods.

4

Which Fortran method? Is this equation stiff? If it’s stiff then you shouldn’t be using Vern9 and would expect this kind of difference.

Are you sure the model is implemented correctly? Plot the output and double check that it is the same. Most of the Fortran solvers have a wrapper, so you can just call for example solve(prob,radau()) and see if the values are exactly the same between the Julia and Fortran version.

5

There is something unusual going on. Even if I reduce the number of runs to just 1, the total time taken is still close to 100 secs. I ran the script in Jupyter 5 times without changing any thing in the script and the same result was obtained.

The differential equations are the same you use in your benchmarks under the name three-body orbits. Integration results are correct. Fortran timing results are also from Verner9 method. Here is my fortran library with which i am comparing: https://github.com/princemahajan/FLINT

6

@Prince_Mahajan I cannot recreate, and I am pretty sure you’re just measuring compile time. This method is pretty heavy on the compile time because of how it performs the kernel generation for GPUs, so that is expected (and there’s a fix we can do for this). When I run your script twice, I get:

result.total_time = 25.036839354
result.total_time = 1.347138324

The first time includes compile time. The second one, without compile time, does not take more than 1.5 seconds.

1 Like
7

That’s exactly what i expect to see. So somehow everytime I am trying to run the code, it is compiling first even though nothing has changed in the script.

I think my diffeq package is up to date but my julia kernel is 1.0.3 I think. I am gonna update every thing and try again.

8

A lot of our optimizations are 1.1+ due to some 1.0 issues. DiffEq will work on 1.0 but that was only for just a bit after our upgrade, and I think as you saw we haven’t backported all of the optimizations to 1.0 yet.

9

Yes after updating to kernel 1.2.0, I am seeing much better results. Thanks for all the replies.

If anyone is interested, here are the new results for a bunch of solvers (compared to FLINT). I ran the script 3 times and these are the numbers from the third run on my computer. Of course, these results apply to this particular set of differential equations and initial conditions.

tol: 1.0e-11
Computational time in seconds

  • DOPRI5: 1.38711e+01 (vs 2.19 Fortran)
  • DP5: 6.26217e+00 (-do-)
  • TSIT5: 7.09554e+00 (-do-)
  • DOP853: 3.98887e+00 (vs 0.73 Fortran)
  • DP8: 3.54820e+00 (-do-)
  • Vern6: 5.17847e+00 (vs 2.0 in Fortran)
  • Vern9: 5.01823e+00 (vs 0.98 in Fortran)
10

What is it doing on 1.1? I am curious why it’s doing so bad on your computer. Was the GC ran before the third solve and was the amount of saving the same between the two?

11

@ChrisRackauckas

What is it doing on 1.1?

I am not sure what you mean here.

I am curious why it’s doing so bad on your computer.

My computer configuration: Intel Core i5, 16 GB RAM, 64-bit Win10

Was the GC ran before the third solve and was the amount of saving the same between the two?

The numbers are pretty much the same on all three runs. I dont think GC was run before the 3rd run because I do not even know how to do that in julia. Basically, script was run 3 times one after the other in the same session without issuing any other command in between.

12

Julia 1.1 is what we’re still testing on since 1.2 came out last week. I don’t know if there’s a performance regression there.

Anything on the saving?

Julia will not necessarily deallocate the cache from the first or second run before starting the third. I know from running your script locally that the third time more than doubles if you the GC of the third run is counted as in the third run, so you might want to GC.gc() so the runs don’t have that cross talk.

Also, can you share the script you are running?

13

Correct me if I am wrong but it seems that the previous results was on 1.0.3:

14

@ChrisRackauckas The script is same as I had posted originally, just added a few more solvers. Here is it again in its entirety

using DifferentialEquations, ODEInterfaceDiffEq,  BenchmarkTools, Libdl, Statistics, PyCall, LinearAlgebra, ParameterizedFunctions, Plots,Printf,Sundials,LSODA;
gr();

###### @show rtol = 1.0./10.0.^(6:1:13);
@show rtol = 1.0./10.0.^(11);
@show atol = rtol*1e-3;



function cr3bpeom(dX, X, p, t)
    @inbounds begin
        r1 = sqrt((X[1] + mu)^2 + X[2]^2)
        r2 = sqrt((X[1] - 1 + mu)^2 + X[2]^2)

        dX[1:3] .= @view X[4:6]
        dX[4] = X[1] + 2*X[5] - (1 - mu)*(X[1] + mu)/r1^3 - mu*(X[1] - 1 + mu)/r2^3
        dX[5] = X[2] - 2*X[4] - (1 - mu)*X[2]/r1^3 - mu*X[2]/r2^3
        dX[6] = 0.0
    end
end



function JacobiConstant(odesol)
    
    X = odesol[1:3,:]
    Xdot = odesol[4:6,:]

    ~,n = size(X)
    
    # distances
    R1 = copy(X)
    R2 = copy(X)
    R1[1,:] = R1[1,:] .+ mu
    R2[1,:] = R2[1,:] .- 1 .+ mu
    r1 = [norm(R1[:,i]) for i = 1:n]
    r2 = [norm(R2[:,i]) for i = 1:n]
    
    # radius and velocity nagnitudes
    R =  [sqrt(sum((X.^2)[:,i])) for i = 1:n]
    V = [sqrt(sum((Xdot.^2)[:,i])) for i = 1:n]

    # pseudo-potential
    Omega = 1/2*(R.^2) .+ (1-mu)./r1 .+ mu./r2
    
    # Jacobi's Constant
    C = (V.^2)/2 .- Omega
    
    return C
end;




const Runs = 1000;

# The 2nd element of each solver sets the dense=true option
# The 3rd element is the name used in results

struct Solver
    solver
    dense::Bool
    name::String
end

ODESolvers = [
                Solver(dopri5(),false,"DOPRI5"),                              # Hairer's DOPRI5
                Solver(DP5(),true,"DP5"),                                     # Julia's implementation of DOPRI5
                Solver(Tsit5(),true,"TSIT5"),                                 # Tsitouras 5/4 RK
                Solver(dop853(),false,"DOP853"),                              # Hairer's DOP853
                Solver(DP8(),true,"DP8"),                                     # Julia's implementation of DOP853
                #Solver(ARKODE(Sundials.Explicit(), order = 8),true,"ARKODE8"),# Sundial's explicit RK
                Solver(Vern6(),true,"Vern6"),                                 # Verner's 7/6 RK
                #Solver(Vern7(),true,"Vern7"),                                 # Verner's 7/6 RK
                #Solver(Vern8(),true,"Vern8"),                                 # Verner's 7/6 RK    
                Solver(Vern9(),true,"Vern9"),                                 # Verner's 9/8 RK
                Solver(ddeabm(),false,"DDEABM"),                              # Shampine's DDEABM
                #Solver(VCABM(),true,"JABM"),      # forcing dt_min           # Julia's implementation of Shampine's DDEABM
                #Solver(CVODE_Adams(max_order = 12),true,"SunAM12") # Too much error!   # Sundial's Adams-Moulton solver
            ]




# test result Struct
struct ODETestResult
    FinalStateErr
    IOMErr
    MeanTime
    OdeSol
end

# Integrate ODE
function FireODE(X0, ODE, tspan, atol, rtol, method, prob, DenseOn)
 if DenseOn == true
        solve(prob,method,abstol=atol, reltol=rtol,timeseries_errors = false, dense_errors=false, dense=false)  
    else 
        solve(prob,method,abstol=atol, reltol=rtol,timeseries_errors = false, dense_errors=false)  
    end
end
            
# Main ODE Integration Method Test function
function FireODEIntTest(X0, ODE, tspan, OdeMethod,  atol, rtol, FinalState, IOM, IOMFunc::Function, DenseOn)
    
    n = length(FinalState)
    # Define ODE problem
    prob = ODEProblem(ODE, X0, tspan);                
        
    # run solver bunch of times for benchmarking
    sol = 0
    rtime = @elapsed for i in 1:Runs
        sol = FireODE(X0, ODE, tspan, atol, rtol, OdeMethod, prob, DenseOn)
    end
    
    # final state error metric
    serr = norm(sol[1:n,end] - FinalState[1:n])
    
    # compute IOM variations
    IOMsol = IOMFunc(sol)
    ierr = abs.(IOMsol .- IOM)
    
    ODETestResult(serr, ierr, rtime, sol)
end
            
# mass-ratio
const mu = 0.012277471::Float64

# Initial conditions
R0 = [0.994::Float64, 0.0, 0.0]
V0 = [0.0, parse(Float64,"-2.00158510637908252240537862224"), 0.0]

# Time period
Period = parse(Float64,"17.0652165601579625588917206249")

# propagate for bunch of orbits
t0 = 0.0::Float64;
tf = 4.0*Period;

tspan = (t0, tf);


X0 = [R0; V0]

op1 = ODEProblem(cr3bpeom,X0,tspan);

nsol = length(ODESolvers)
ntol = length(rtol)

SolResults3B = Matrix(undef, nsol, ntol)

for itr in 1:ntol
    for ctr in 1:nsol
        print(string("Running: ",ODESolvers[ctr].name, "\n"))
        SolResults3B[ctr,itr] = FireODEIntTest(X0, cr3bpeom, tspan, ODESolvers[ctr].solver,  atol[itr], rtol[itr], [R0; V0], 0, JacobiConstant, ODESolvers[ctr].dense);
    end
    print(string("tol: ", rtol[itr],"\n"))
end

print("Computational time vs errors for the tightest tolerance\n=======================================================\n")
for ctr in 1:nsol
    @Printf.printf("%7s: %.5e, %.5e\n",ODESolvers[ctr].name,SolResults3B[ctr,end].MeanTime,SolResults3B[ctr,end].FinalStateErr)
end
15

I profiled it and the view in the differential equation was surprisingly unoptimized, so

function cr3bpeom(dX, X, p, t)
    @inbounds begin
        r1 = sqrt((X[1] + mu)^2 + X[2]^2)
        r2 = sqrt((X[1] - 1 + mu)^2 + X[2]^2)

        dX[1] = X[4]
        dX[2] = X[5]
        dX[3] = X[6]
        dX[4] = X[1] + 2*X[5] - (1 - mu)*(X[1] + mu)/r1^3 - mu*(X[1] - 1 + mu)/r2^3
        dX[5] = X[2] - 2*X[4] - (1 - mu)*X[2]/r1^3 - mu*X[2]/r2^3
        dX[6] = 0.0
    end
end

sped it up quite a bit. One thing that was happening was that the compilation was inside of the loop, so the ODE solver needs to be run once before the loop if you want to exclude the compile time. And as I expected Julia’s GC is what you’re hitting in the loop, so I get another reduction by doing:

function FireODEIntTest(X0, ODE, tspan, OdeMethod,  atol, rtol, FinalState, IOM, IOMFunc::Function, DenseOn)

    n = length(FinalState)
    # Define ODE problem
    prob = ODEProblem(ODE, X0, tspan);

    # run solver bunch of times for benchmarking
    sol = 0
    rtime = 0.0

    # Compile before timing
    sol = FireODE(X0, ODE, tspan, atol, rtol, OdeMethod, prob, DenseOn)

    for i in 1:Runs
        GC.gc()
        rtime += @elapsed sol = FireODE(X0, ODE, tspan, atol, rtol, OdeMethod, prob, DenseOn)
    end

    # final state error metric
    serr = norm(sol[1:n,end] - FinalState[1:n])

    # compute IOM variations
    IOMsol = IOMFunc(sol)
    ierr = abs.(IOMsol .- IOM)

    ODETestResult(serr, ierr, rtime, sol)
end

Before:

    DP5: 8.90705e+00, 1.71606e+00
  TSIT5: 9.15936e+00, 7.33951e-02
    DP8: 6.43439e+00, 6.66641e-01
  Vern6: 8.61046e+00, 3.99717e-01
  Vern9: 1.09995e+01, 5.92027e-02

After:

    DP5: 4.15515e+00, 1.71606e+00
  TSIT5: 4.82140e+00, 7.33951e-02
    DP8: 1.15518e+00, 6.66641e-01
  Vern6: 3.15056e+00, 3.99717e-01
  Vern9: 1.02262e+00, 5.92027e-02

That amount of time reduction probably puts it in line with Fortran, so that’s where the difference between the two likely were (I can’t find out how to run Flint on my computer: would you be willing to make a Julia wrapper? It would be nice to benchmark it!)

1 Like
16

Note that if you use StaticArrays you also get right down to there:

using DifferentialEquations, ODEInterfaceDiffEq,  BenchmarkTools, Libdl, Statistics,
      PyCall, LinearAlgebra, ParameterizedFunctions, Plots,Printf,Sundials,LSODA,
      StaticArrays
gr();

###### @show rtol = 1.0./10.0.^(6:1:13);
@show rtol = 1.0./10.0.^(11);
@show atol = rtol*1e-3;



function cr3bpeom(X, p, t)
    @inbounds begin
        r1 = sqrt((X[1] + mu)^2 + X[2]^2)
        r2 = sqrt((X[1] - 1 + mu)^2 + X[2]^2)

        dX4 = X[1] + 2*X[5] - (1 - mu)*(X[1] + mu)/r1^3 - mu*(X[1] - 1 + mu)/r2^3
        dX5 = X[2] - 2*X[4] - (1 - mu)*X[2]/r1^3 - mu*X[2]/r2^3

        @SVector [X[4],X[5],X[6],dX4,dX5,0.0]
    end
end



function JacobiConstant(odesol)

    X = odesol[1:3,:]
    Xdot = odesol[4:6,:]

    ~,n = size(X)

    # distances
    R1 = copy(X)
    R2 = copy(X)
    R1[1,:] = R1[1,:] .+ mu
    R2[1,:] = R2[1,:] .- 1 .+ mu
    r1 = [norm(R1[:,i]) for i = 1:n]
    r2 = [norm(R2[:,i]) for i = 1:n]

    # radius and velocity nagnitudes
    R =  [sqrt(sum((X.^2)[:,i])) for i = 1:n]
    V = [sqrt(sum((Xdot.^2)[:,i])) for i = 1:n]

    # pseudo-potential
    Omega = 1/2*(R.^2) .+ (1-mu)./r1 .+ mu./r2

    # Jacobi's Constant
    C = (V.^2)/2 .- Omega

    return C
end;




const Runs = 1000;

# The 2nd element of each solver sets the dense=true option
# The 3rd element is the name used in results

struct Solver
    solver
    dense::Bool
    name::String
end

ODESolvers = [
                #Solver(dopri5(),false,"DOPRI5"),                              # Hairer's DOPRI5
                Solver(DP5(),true,"DP5"),                                     # Julia's implementation of DOPRI5
                Solver(Tsit5(),true,"TSIT5"),                                 # Tsitouras 5/4 RK
                #Solver(dop853(),false,"DOP853"),                              # Hairer's DOP853
                Solver(DP8(),true,"DP8"),                                     # Julia's implementation of DOP853
                #Solver(ARKODE(Sundials.Explicit(), order = 8),true,"ARKODE8"),# Sundial's explicit RK
                Solver(Vern6(),true,"Vern6"),                                 # Verner's 7/6 RK
                #Solver(Vern7(),true,"Vern7"),                                 # Verner's 7/6 RK
                #Solver(Vern8(),true,"Vern8"),                                 # Verner's 7/6 RK
                Solver(Vern9(),true,"Vern9"),                                 # Verner's 9/8 RK
                #Solver(ddeabm(),false,"DDEABM"),                              # Shampine's DDEABM
                #Solver(VCABM(),true,"JABM"),      # forcing dt_min           # Julia's implementation of Shampine's DDEABM
                #Solver(CVODE_Adams(max_order = 12),true,"SunAM12") # Too much error!   # Sundial's Adams-Moulton solver
            ]




# test result Struct
struct ODETestResult
    FinalStateErr
    IOMErr
    MeanTime
    OdeSol
end

# Integrate ODE
function FireODE(X0, ODE, tspan, atol, rtol, method, prob, DenseOn)
 if DenseOn == true
        solve(prob,method,abstol=atol, reltol=rtol,timeseries_errors = false, dense_errors=false, dense=false)
    else
        solve(prob,method,abstol=atol, reltol=rtol,timeseries_errors = false, dense_errors=false)
    end
end

# Main ODE Integration Method Test function
function FireODEIntTest(X0, ODE, tspan, OdeMethod,  atol, rtol, FinalState, IOM, IOMFunc::Function, DenseOn)

    n = length(FinalState)
    # Define ODE problem
    prob = ODEProblem(ODE, X0, tspan);

    # run solver bunch of times for benchmarking
    sol = 0

    # Remove compilation run
    sol = FireODE(X0, ODE, tspan, atol, rtol, OdeMethod, prob, DenseOn)

    rtime = @elapsed for i in 1:Runs
        sol = FireODE(X0, ODE, tspan, atol, rtol, OdeMethod, prob, DenseOn)
    end

    # final state error metric
    serr = norm(sol[1:n,end] - FinalState[1:n])

    # compute IOM variations
    IOMsol = IOMFunc(sol)
    ierr = abs.(IOMsol .- IOM)

    ODETestResult(serr, ierr, rtime, sol)
end

# mass-ratio
const mu = 0.012277471::Float64

# Initial conditions
R0 = @SVector [0.994::Float64, 0.0, 0.0]
V0 = @SVector [0.0, parse(Float64,"-2.00158510637908252240537862224"), 0.0]

# Time period
Period = parse(Float64,"17.0652165601579625588917206249")

# propagate for bunch of orbits
t0 = 0.0::Float64;
tf = 4.0*Period;

tspan = (t0, tf);


X0 = [R0; V0]

op1 = ODEProblem(cr3bpeom,X0,tspan);

nsol = length(ODESolvers)
ntol = length(rtol)

SolResults3B = Matrix(undef, nsol, ntol)

for itr in 1:ntol
    for ctr in 1:nsol
        print(string("Running: ",ODESolvers[ctr].name, "\n"))
        SolResults3B[ctr,itr] = FireODEIntTest(X0, cr3bpeom, tspan, ODESolvers[ctr].solver,  atol[itr], rtol[itr], [R0; V0], 0, JacobiConstant, ODESolvers[ctr].dense);
    end
    print(string("tol: ", rtol[itr],"\n"))
end

print("Computational time vs errors for the tightest tolerance\n=======================================================\n")
for ctr in 1:nsol
    @Printf.printf("%7s: %.5e, %.5e\n",ODESolvers[ctr].name,SolResults3B[ctr,end].MeanTime,SolResults3B[ctr,end].FinalStateErr)
end

Your old code before:

    DP5: 8.90705e+00, 1.71606e+00
  TSIT5: 9.15936e+00, 7.33951e-02
    DP8: 6.43439e+00, 6.66641e-01
  Vern6: 8.61046e+00, 3.99717e-01
  Vern9: 1.09995e+01, 5.92027e-02

This code:

    DP5: 3.88171e+00, 1.72335e+00
  TSIT5: 5.00398e+00, 6.24401e-02
    DP8: 1.08576e+00, 6.16423e-01
  Vern6: 2.71387e+00, 5.33841e-01
  Vern9: 1.14183e+00, 5.07811e-02
1 Like
17

@ChrisRackauckas Good work on your part!

I tried your changes and removing view and running the solver once before the timing loop has made a difference. Seems like to match compiled code, the lesson is stay away from array slicing as much as possible and compile before hand by running the same code once (which may not be possible in more realistic scenarios where initial conditions are changing between the runs). On my computer:

Before

 DOPRI5: 1.43511e+01, 1.71222e+00
    DP5: 6.35068e+00, 1.71606e+00
  TSIT5: 7.22626e+00, 7.33951e-02
 DOP853: 4.05822e+00, 7.84506e-01
    DP8: 3.58036e+00, 6.66641e-01
  Vern6: 5.21433e+00, 3.99717e-01
  Vern9: 4.98851e+00, 5.92027e-02
 DDEABM: 8.01595e+00, 2.16018e+00

After

DOPRI5: 1.22066e+01, 1.71222e+00 (2.19 s)
    DP5: 3.78245e+00, 1.71606e+00 (-do-) 
  TSIT5: 4.48538e+00, 7.33951e-02 (-do-)
 DOP853: 3.16085e+00, 7.84506e-01 (0.734 s)
    DP8: 9.13293e-01, 6.66641e-01 (-do-)
  Vern6: 2.52114e+00, 3.99717e-01 (2.06 s)
  Vern9: 9.01457e-01, 5.92027e-02 (0.984 s) *
 DDEABM: 7.53781e+00, 2.16018e+00

(I can’t find out how to run Flint on my computer: would you be willing to make a Julia wrapper? It would be nice to benchmark it!)

The travis file on github has all the command that can be replicated to compile and run FLINT on Linux/Win. You only need cmake and gfortran.

A wrapper in Julia and Python is actually in my plans. FLINT is fairly a new code and needs to be tested thoroughly, which is my priority at this time.

I did not have garbage collection command in my original script, so adding GC() between the runs and even outside the loop made it a lot slower in my trials.

For static arrays, I got this warning while running DOPRI5. You may want to look at it. (I am using kernel 1.2.0)

@ ODEInterfaceDiffEq C:\Users\bmahajan.julia\packages\ODEInterfaceDiffEq\ZJ8z2\src\solve.jl:169
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ ODEInterfaceDiffEq

18

Initial conditions don’t cause a recompile. It’s just the first time you solve with a method in a new Julia session. You can change initial conditions, parameters, etc (run a parameter estimation). The only thing that matters is that every time a new code is generated (new types are used), then you have to recompile.

And the issue was that the view was allocating for some reason. Usually Julia’s compiler catches this and optimizes it away. I’m not sure why it wasn’t in this case :man_shrugging:. But yeah, at the end of the day Julia will get to the speed of compiled code if you write it to be like the compiled code, which isn’t always the nicest form. I will say that Fortran does look pretty high level for being a compiled language, so it is a fun one to code in.

Yes, only the pure-Julia methods can use static arrays.

BTW, it looks like you’re doing Verner robust and efficient? How’re you seeing the robust ones doing? I did a bunch of tests with a tableau-based implementation awhile back but never did too many benchmarks after optimizing everything else. I’ve wanted to revisit them.

19

Initial conditions don’t cause a recompile. It’s just the first time you solve with a method in a new Julia session. You can change initial conditions, parameters, etc (run a parameter estimation). The only thing that matters is that every time a new code is generated (new types are used), then you have to recompile.

That is indeed a plus point. I will try that with random initial conditions when I get time.

BTW, it looks like you’re doing Verner robust and efficient? How’re you seeing the robust ones doing? I did a bunch of tests with a tableau-based implementation awhile back but never did too many benchmarks after optimizing everything else. I’ve wanted to revisit them.

Actually, I used efficient version of the coefficients for Verner6 and robust coefficients for Verner9 for no particular reason. I have not compared the same algorithms with two different sets of coefficients to see the difference. Even though adding new coefficients in FLINT is fairly easy. Some day…