sundials doesn't GPU

[vpuri3]

using CUDA, Sundials

N  = 4
f  = (u, p, t) -> zero(u)
u0 = CUDA.rand(N) 
ts = (0f0, 1f0)

prob = ODEProblem(f, u0, ts)
alg  = CVODE_BDF()
sol  = solve(prob, alg)

julia> include("examples/ad/sundials_gpu.jl")                                     [42/1828]
ERROR: LoadError: MethodError: no method matching Ptr{Sundials._generic_N_Vector}(::Vector{
Float32})                                                                                  
Closest candidates are:                                                                    
  Ptr{Sundials._generic_N_Vector}(::Sundials.NVector) at ~/.julia/packages/Sundials/k9hc3/s
rc/nvector_wrapper.jl:29                                                                   
  Ptr{T}() where T at boot.jl:781                                                          
  Ptr{T}(::Union{Int64, UInt64, Ptr}) where T at boot.jl:780                               
Stacktrace:                                                                                
  [1] convert(#unused#::Type{Ptr{Sundials._generic_N_Vector}}, v::Vector{Float32})         
    @ Sundials ~/.julia/packages/Sundials/k9hc3/src/nvector_wrapper.jl:70                  
  [2] __init(prob::ODEProblem{CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, Tuple{Float32, Fl
oat32}, false, SciMLBase.NullParameters, ODEFunction{false, var"#196#197", UniformScaling{B
ool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Noth
ing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Uni
on{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::CVODE_BDF{:New
ton, :Dense, Nothing, Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}; 
verbose::Bool, callback::Nothing, abstol::Float64, reltol::Float64, saveat::Vector{Float64}
, tstops::Vector{Float64}, maxiters::Int64, dt::Nothing, dtmin::Float64, dtmax::Float64, ti
meseries_errors::Bool, dense_errors::Bool, save_everystep::Bool, save_idxs::Nothing, save_o
n::Bool, save_start::Bool, save_end::Bool, dense::Bool, progress::Bool, progress_name::Stri
ng, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), save_timeseries::Nothing
, advance_to_tstop::Bool, stop_at_next_tstop::Bool, userdata::Nothing, alias_u0::Bool, kwar
gs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})  
    @ Sundials ~/.julia/packages/Sundials/k9hc3/src/common_interface/solve.jl:232          
  [3] __init(prob::ODEProblem{CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, Tuple{Float32, Fl
oat32}, false, SciMLBase.NullParameters, ODEFunction{false, var"#196#197", UniformScaling{B
ool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Noth
ing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Uni
on{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::CVODE_BDF{:New
ton, :Dense, Nothing, Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any})
    @ Sundials ~/.julia/packages/Sundials/k9hc3/src/common_interface/solve.jl:87           
  [4] __solve(prob::ODEProblem{CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, Tuple{Float32, F
loat32}, false, SciMLBase.NullParameters, ODEFunction{false, var"#196#197", UniformScaling{
Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Not
hing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Un
ion{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::CVODE_BDF{:Ne
wton, :Dense, Nothing, Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any},
 recompile::Type{Val{true}}; calculate_error::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tup
le{}, NamedTuple{(), Tuple{}}})                                                            
    @ Sundials ~/.julia/packages/Sundials/k9hc3/src/common_interface/solve.jl:14           
  [5] __solve(prob::ODEProblem{CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, Tuple{Float32, F
loat32}, false, SciMLBase.NullParameters, ODEFunction{false, var"#196#197", UniformScaling{
Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Not
hing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Un
ion{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::CVODE_BDF{:Ne
wton, :Dense, Nothing, Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any},
 recompile::Type{Val{true}}) (repeats 5 times)                      
    @ Sundials ~/.julia/packages/Sundials/k9hc3/src/common_interface/solve.jl:3            
  [6] solve_call(_prob::ODEProblem{CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, Tuple{Float3
2, Float32}, false, SciMLBase.NullParameters, ODEFunction{false, var"#196#197", UniformScal
ing{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing,
 Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol
, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::CVODE_BD
F{:Newton, :Dense, Nothing, Nothing}; merge_callbacks::Bool, kwargshandle::KeywordArgError,
 kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/xOkrD/src/solve.jl:439
  [7] solve_call(_prob::ODEProblem{CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, Tuple{Float3
2, Float32}, false, SciMLBase.NullParameters, ODEFunction{false, var"#196#197", UniformScal
ing{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing,
 Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol
, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::CVODE_BD
F{:Newton, :Dense, Nothing, Nothing})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/xOkrD/src/solve.jl:409

[ChrisRackauckas]

That is expected. Sundials is a different repository from OrdinaryDiffEq.jl. It's a completely different set of solvers and in fact it's a set of wrappers to the C/Fortran SUNDIALS library. As such, it won't work with Julia GPU kernels. In theory we could wrapper the CUDA N_Vector and pass pointers or something, but at this time it's not in scope for the people developing the library.

[ViralBShah]

There is significant improvement in building binaries that can use the CUDA stack, which should pave the path for this if someone wants to work on it.

[vpuri3]

@abillscmu if you are still interested.

@ViralBShah is support for distributed computing with Sundials within scope?

[ViralBShah]

I believe there is now a recipe for building MPI libraries in BinaryBuilder - so it is certainly possible. Should be a separate JLL whenever we decide to do it.

[ViralBShah]

GPU build is in JuliaPackaging/Yggdrasil#12067.

I do believe this will need perhaps more wrapping of the nvector GPU stuff in Sundials.jl, to make sure everything stays on the GPU. That needs some work, and then one should be able to do RHS that are GPU kernels and such.