# DIfferentialEquations and GPU

Hello,

According to the following blog post http://www.stochasticlifestyle.com/solving-systems-stochastic-pdes-using-gpus-julia/ it is very easy to turn on GPU parallelism with DifferentialEquations.jl. The only thing one needs to do is to replace ordinary arrays with CLArrays or CuArrays. However, the following code shows that the GPU version is orders of magnitude slower:

``````using DifferentialEquations
using CuArrays

function func(du, u, p, t)
@inbounds @. du = u
return nothing
end

tspan = (0., 10.)
alg = BS3()

# ******************************************************************************
println("CPU:")

u0 = fill(1., 500)
prob = ODEProblem(func, u0, tspan)

@time sol = solve(prob, alg, saveat=0.01)
@time sol = solve(prob, alg, saveat=0.01)
@time sol = solve(prob, alg, saveat=0.01)

# ******************************************************************************
println("GPU:")

# CuArrays.allowscalar(false)

u0_gpu = CuArray(convert(Array{Float32, 1}, u0))
prob_gpu = ODEProblem(func, u0_gpu, tspan)

@time sol_gpu = solve(prob_gpu, alg, saveat=0.01)
@time sol_gpu = solve(prob_gpu, alg, saveat=0.01)
@time sol_gpu = solve(prob_gpu, alg, saveat=0.01)
``````

The output on my machine with Nvidia GTX 1070 is the following:

``````CPU:
3.627708 seconds (14.59 M allocations: 756.548 MiB, 5.28% gc time)
0.002929 seconds (1.11 k allocations: 4.082 MiB)
0.003089 seconds (1.11 k allocations: 4.082 MiB)
GPU:
32.020281 seconds (43.64 M allocations: 2.092 GiB, 1.71% gc time)
24.563494 seconds (20.42 M allocations: 951.305 MiB, 0.49% gc time)
27.890973 seconds (20.42 M allocations: 951.711 MiB, 0.48% gc time)
``````

Can you please explain me what am I doing wrong?

Thank you.

Did you run with that line uncommented? If that errors, thereâ€™s a scalar operation that isnâ€™t supported by CuArrays.jl (yet) and will completely ruin performance.

Yes, it gives the error.

Well, then I will try to reformulate the question. Are there any examples which show how DifferentialEquations can be speed up using GPU? So far I have found only the above blog post which discusses this possibility.

Maybe not completely what you request but see

not to my knowledge. A few points though:

1. it will be difficult to speed up real problems on gpu, even big ones, without having to write a custom kernel somewhere
differential equations solvers are typically iterative in nature, so the speedup from the GPU needs to be the individual array operations themselves. This is a hard problem to speed up on GPU unless the operations involved are very large and very computationally intensive
2. for your example, the GPU probably wonâ€™t beat the CPU by much even if `func` is written as a custom kernel. youâ€™re pretty much just testing the L2 cache bandwidth of a single CPU core vs L2 cache on a GPU.

Thank you all!

Youâ€™re timing the worst case scenario for GPUs, which is where all of your kernel calls are trivial. Make your ODE something with a matrix multiplication or something else with a lot of linear algebra and youâ€™ll see a different result. BLAS1 computations (â€śelement-wiseâ€ť or broadcast) is usually not generally enough to warrant GPU usage.

Weâ€™re working on this. We will have a GSoC that optimizes GPU code and gets it running with methods for stiff ODEs, focusing on getting solvers for PDEs.

1 Like

Hey @ChrisRackauckas,

can we now solve this issue?

``````
CPU:
2.448770 seconds (8.22 M allocations: 409.291 MiB, 7.44% gc time)
0.002015 seconds (1.09 k allocations: 4.070 MiB)
0.001831 seconds (1.09 k allocations: 4.070 MiB)
GPU:
26.472584 seconds (41.66 M allocations: 2.079 GiB, 6.59% gc time)
0.221883 seconds (140.06 k allocations: 15.828 MiB)
0.269016 seconds (170.76 k allocations: 16.724 MiB, 4.14% gc time)
``````

You had time as Float64, and the kernels werenâ€™t filled. Hereâ€™s an example:

``````using OrdinaryDiffEq
using CuArrays
CuArrays.allowscalar(false)

function func(du, u, p, t)
@inbounds @. du = u
return nothing
end

tspan = (0f0, 10f0)
alg = BS3()
u0 = fill(1f0, Int(5e5))

# ******************************************************************************
println("CPU:")
prob = ODEProblem(func, u0, tspan)

@time sol = solve(prob, alg, saveat=1f-2)
@time sol = solve(prob, alg, saveat=1f-2)
@time sol = solve(prob, alg, saveat=1f-2)

# ******************************************************************************
println("GPU:")

# CuArrays.allowscalar(false)

u0_gpu = CuArray(convert(Array{Float32, 1}, u0))
prob_gpu = ODEProblem(func, u0_gpu, tspan)

sol_gpu = nothing; GC.gc(true); CuArrays.reclaim()
@time sol_gpu = solve(prob_gpu, alg, saveat=1f-2)
sol_gpu = nothing; GC.gc(true); CuArrays.reclaim()
@time sol_gpu = solve(prob_gpu, alg, saveat=1f-2)
sol_gpu = nothing; GC.gc(true); CuArrays.reclaim()
@time sol_gpu = solve(prob_gpu, alg, saveat=1f-2)
``````
``````CPU:
3.288380 seconds (8.77 M allocations: 2.316 GiB, 12.26% gc time)
0.912710 seconds (2.10 k allocations: 1.891 GiB)
0.973323 seconds (2.10 k allocations: 1.891 GiB, 8.27% gc time)
GPU:
9.955818 seconds (27.20 M allocations: 1.347 GiB, 3.39% gc time)
0.462713 seconds (113.43 k allocations: 6.869 MiB)
0.452901 seconds (113.44 k allocations: 6.878 MiB)
``````
3 Likes

Now itâ€™s solved, Thanks!

Hello Chris,

Thank you very much for making it work.

Can you, please, comment on this line:

``````sol_gpu = nothing; GC.gc(true); CuArrays.reclaim()
``````

What exactly happens here and do I need to do the same in the actual code before each run of solve?

Oh, that was just for clearing memory when testing

I notice that you use `OrdinaryDiffEq` instead of `DifferentialEquations`. Does that make any difference? Is that means `DifferentialEquations` cant be used this way?

Hi,