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.