I would like to farm out solving systems of ODEs onto GPUs, in a 'trivially parallelisable' setting. For example, doing a sensitivity analysis with 512 different parameter sets.

Ideally I want to do ODE solving with a smart adaptive timestep solver like CVODE, rather than a fixed timestep like Forward Euler, but running it on an NVIDIA GPU instead of CPU.

Has anyone done this? Are there libraries for it?

  • $\begingroup$ Hence the brackets! I'm considering an operator-splitting based technique (cardiac electrophysiology simulations), where you solve ODEs at nodes to get a source term for PDE, then change an ODE parameter for next iteration. $\endgroup$
    – mirams
    Commented Jan 12, 2014 at 9:34
  • 1
    $\begingroup$ Maybe related? What's the state of the art in parallel ODE methods? $\endgroup$
    – Kirill
    Commented Sep 24, 2015 at 18:52
  • 2
    $\begingroup$ It's important to specify whether you want to use the same time-stepping for every ODE or not. $\endgroup$ Commented Sep 30, 2015 at 16:15
  • $\begingroup$ Also, if you're specifically interested in the bidomain (or monodomain) equations, you might want to take a look at how CARP does it. $\endgroup$ Commented Sep 30, 2015 at 16:17
  • $\begingroup$ Different timesteps, if the method is adaptive then it will need them for different parameter sets... thanks for the link to what CARP is doing - a fixed timestep Rush Larsen ODE solver if I read it correctly. $\endgroup$
    – mirams
    Commented Sep 30, 2015 at 22:27

2 Answers 2


You may want to look into Boost's odeint library and Thrust. They can be combined as discussed here.

  • $\begingroup$ This seems to be a bit different - solving massive ODE systems on the GPU in parallel (with communication). That link says "We have experienced that the vector size over which is parallelized should be of the order of 10^6 to make full use of the GPU.". I'm looking for a nice way of farming out O(10) or O(100) vector sized trivially parallelisable ODE solves... $\endgroup$
    – mirams
    Commented Jan 19, 2014 at 13:05
  • $\begingroup$ Have you thought on writing directly in cuda or openCL? If I undertood right, what you are doing is iterating over some ODE equation in each thread separately, it shouldn't be difficult to write it directly. $\endgroup$
    – Hydro Guy
    Commented Sep 24, 2015 at 13:22
  • $\begingroup$ I imagine it would be possible to code a Forward Euler or other fixed timestep method, where every GPU process uses the same timestep, fairly easily, I'd like to know whether anyone has managed to get adaptive timestepping like CVODE working, or whether this is impossible to make efficient on a GPGPU? $\endgroup$
    – mirams
    Commented Sep 25, 2015 at 7:24
  • $\begingroup$ the problem with gpu is that you need to write data-parallel code. If you write the same adaptive routine but absorbing all that flexibility on the values of some parameters, probably it's possible to code it efficiently on gpu. That also means that you can't use branching on instructions, which is probably what you think that would make it impossible to do it. $\endgroup$
    – Hydro Guy
    Commented Sep 29, 2015 at 3:30
  • 1
    $\begingroup$ @mirams there's an example for odeint that covers exactly what you are looking for: boost.org/doc/libs/1_59_0/libs/numeric/odeint/doc/html/…, see also github.com/boostorg/odeint/blob/master/examples/thrust/…. Also, odeint supports adaptive multistep methods as in CVODE: github.com/boostorg/odeint/blob/master/examples/… $\endgroup$ Commented Sep 30, 2015 at 16:06

DifferentialEquations.jl library is a library for a high level language (Julia) which has tools for automatically transforming the ODE system to an optimized version for parallel solution on GPUs. There are two forms of parallelism that can be employed: array-based parallelism for large ODE systems and parameter parallelism for parameter studies on relatively small (<100) ODE systems. It supports high order implicit and explicit methods and routinely outperforms or matches other systems in benchmarks (at the very least, it wraps the others so it's easy to check and use them!)

For this specific functionality, you may want to take a look at DiffEqGPU.jl which is the module for automated parameter parallelism. The DifferentialEquations.jl library has functionality for parallel parameter studies, and this module augments the existing configurations to make the study happen automatically in parallel. What one does is transforms their existing ODEProblem (or other DEProblem like SDEProblem) into an EnsembleProblem and specify with a prob_func how the other problems are generated from the prototype. The following solves 10,000 trajectories of the Lorenz equation on the GPU with a high order explicit adaptive method:

using OrdinaryDiffEq, DiffEqGPU
function lorenz(du,u,p,t)
 @inbounds begin
     du[1] = p[1]*(u[2]-u[1])
     du[2] = u[1]*(p[2]-u[3]) - u[2]
     du[3] = u[1]*u[2] - p[3]*u[3]

u0 = Float32[1.0;0.0;0.0]
tspan = (0.0f0,100.0f0)
p = (10.0f0,28.0f0,8/3f0)
prob = ODEProblem(lorenz,u0,tspan,p)
prob_func = (prob,i,repeat) -> remake(prob,p=rand(Float32,3).*p)
monteprob = EnsembleProblem(prob, prob_func = prob_func)
@time sol = solve(monteprob,Tsit5(),EnsembleGPUArray(),trajectories=10_000,saveat=1.0f0)

Notice the user needs to write no GPU code, and with a single RTX 2080 this benchmarks as a 5x improvement over using a 16 core Xeon machine with multithreaded parallelism. One can then check out the README for how to do things like utilize multiple GPUs and doing multiprocessing+GPUs for utilizing a full cluster of GPUs simultaneously. Note that switching to multithreading instead of GPUs is one line change: EnsembleThreads() instead of EnsembleGPUArray().

Then for implicit solvers, the same interface holds. For example, the following uses high order Rosenbrock and implicit Runge-Kutta methods:

function lorenz_jac(J,u,p,t)
 @inbounds begin
     σ = p[1]
     ρ = p[2]
     β = p[3]
     x = u[1]
     y = u[2]
     z = u[3]
     J[1,1] = -σ
     J[2,1] = ρ - z
     J[3,1] = y
     J[1,2] = σ
     J[2,2] = -1
     J[3,2] = x
     J[1,3] = 0
     J[2,3] = -x
     J[3,3] = -β

function lorenz_tgrad(J,u,p,t)

func = ODEFunction(lorenz,jac=lorenz_jac,tgrad=lorenz_tgrad)
prob_jac = ODEProblem(func,u0,tspan,p)
monteprob_jac = EnsembleProblem(prob_jac, prob_func = prob_func)

@time solve(monteprob_jac,Rodas5(linsolve=LinSolveGPUSplitFactorize()),EnsembleGPUArray(),dt=0.1,trajectories=10_000,saveat=1.0f0)
@time solve(monteprob_jac,TRBDF2(linsolve=LinSolveGPUSplitFactorize()),EnsembleGPUArray(),dt=0.1,trajectories=10_000,saveat=1.0f0)

While this form requires that you give a Jacobian in order to be used on the GPU (currently, will be fixed soon), the DifferentialEquations.jl documentation demonstrates how to do automatic symbolic Jacobian calculations on numerically defined functions, so there is still no manual labor here. I would highly recommend these algorithms because the branching logic of a method like CVODE generally causes thread desync and doesn't seem to perform as well as a Rosenbrock method in these types of scenarios anyways.

By using DifferentialEquations.jl, you also get access to the full library, which includes functionality like global sensitivity analysis which can make use of this GPU-acceleration. It is also compatible with dual numbers for fast local sensitivity analysis. The GPU-based code gets all of the features of DifferentialEquations.jl, like the event handling and the large set of ODE solvers which are optimized for different types of problems, meaning it's not just a simple one-off GPU ODE solver but instead a part of a fully featured system which also happens to have efficient GPU support.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.