NonlinearSolve.jl

NonlinearSolve.jl copied to clipboard

It seems that (parts of) the library won't work with 1D CUDA arrays, but with higher dimensional arrays it does.

using NonlinearSolve, CUDA
CUDA.allowscalar(false) 
f(u,p) = u .* u .- 2
u0 = [1.0, 1.0]
probN = NonlinearProblem{false}(f, u0)
solver = solve(probN, NewtonRaphson(), tol = 1e-9)

works

u0 = reshape(cu(u0),:,1)
probN = NonlinearProblem{false}(f, u0)
solver = @time solve(probN, NewtonRaphson(), tol = 1e-9)

works

u0 = cu([1.0, 1.0])
probN = NonlinearProblem{false}(f, u0)
solver = solve(probN, NewtonRaphson(), tol = 1e-9)

doesn't work, with error:

MethodError: no method matching copy_cublasfloat(::CuArray{Float32, 1, CUDA.Mem.DeviceBuffer})
Closest candidates are:
  copy_cublasfloat(::CuArray{T, 2}) where T at ~/.julia/packages/CUDA/Uurn4/lib/cusolver/linalg.jl:6

Stacktrace:
 [1] \(_A::CuArray{Float32, 2, CUDA.Mem.DeviceBuffer}, _B::CuArray{Float32, 1, CUDA.Mem.DeviceBuffer})
   @ CUDA.CUSOLVER ~/.julia/packages/CUDA/Uurn4/lib/cusolver/linalg.jl:24
 [2] macro expansion
   @ ~/.julia/packages/Setfield/AS2xF/src/sugar.jl:197 [inlined]
 [3] perform_step(solver::NonlinearSolve.NewtonImmutableSolver{NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}, CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, SciMLBase.NullParameters, typeof(NonlinearSolve.DEFAULT_NORM), Float64, Nothing, NonlinearProblem{CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, false, SciMLBase.NullParameters, NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}}, alg::NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}, #unused#::Val{false})
   @ NonlinearSolve ~/.julia/packages/NonlinearSolve/ZxbEJ/src/raphson.jl:58
 [4] solve!(solver::NonlinearSolve.NewtonImmutableSolver{NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}, CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, SciMLBase.NullParameters, typeof(NonlinearSolve.DEFAULT_NORM), Float64, Nothing, NonlinearProblem{CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, false, SciMLBase.NullParameters, NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}})
   @ NonlinearSolve ~/.julia/packages/NonlinearSolve/ZxbEJ/src/solve.jl:63
 [5] #solve#6
   @ ~/.julia/packages/NonlinearSolve/ZxbEJ/src/solve.jl:5 [inlined]
 [6] top-level scope
   @ In[7]:3
 [7] eval
   @ ./boot.jl:373 [inlined]
 [8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1196

Thanks! I'll take a look at this when I get around to reworking this library to use LinearSolve.jl