ChrisRackauckas
commented
1 week ago Therefor the solver shouldn't compute the full matrix exponential of the given matrix operator.
Did you set krylov=true?
Open mcarmesin opened 1 week ago
ChrisRackauckas
commented
1 week ago Therefor the solver shouldn't compute the full matrix exponential of the given matrix operator.
Did you set krylov=true?
mcarmesin
commented
1 week ago solve(prob, LinearExponential(krylov=:simple)) fails with a different error:
ERROR: This object is not a GPU array
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] backend(::Type)
@ GPUArraysCore ~/.julia/packages/GPUArraysCore/GMsgk/src/GPUArraysCore.jl:225
[3] backend(::Type{SubArray{Float32, 1, Matrix{Float32}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}})
@ GPUArraysCore ~/.julia/packages/GPUArraysCore/GMsgk/src/GPUArraysCore.jl:229
[4] backend(x::SubArray{Float32, 1, Matrix{Float32}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true})
@ GPUArraysCore ~/.julia/packages/GPUArraysCore/GMsgk/src/GPUArraysCore.jl:226
[5] _copyto!
@ ~/.julia/packages/GPUArrays/qt4ax/src/host/broadcast.jl:78 [inlined]
[6] materialize!
@ ~/.julia/packages/GPUArrays/qt4ax/src/host/broadcast.jl:38 [inlined]
[7] materialize!
@ ./broadcast.jl:911 [inlined]
[8] firststep!(Ks::KrylovSubspace{…}, V::SubArray{…}, H::SubArray{…}, b::CuArray{…})
@ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/xLH9y/src/arnoldi.jl:171
[9] lanczos!(Ks::KrylovSubspace{…}, A::CUDA.CUSPARSE.CuSparseMatrixCSC{…}, b::CuArray{…}; tol::Float64, m::Int64, opnorm::Nothing, init::Int64, t::Float64, mu::Float64, l::Int64)
@ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/xLH9y/src/arnoldi.jl:331
[10] lanczos!
@ ~/.julia/packages/ExponentialUtilities/xLH9y/src/arnoldi.jl:318 [inlined]
[11] arnoldi!(Ks::KrylovSubspace{…}, A::CUDA.CUSPARSE.CuSparseMatrixCSC{…}, b::CuArray{…}; tol::Float64, m::Int64, ishermitian::Bool, opnorm::Function, iop::Int64, init::Int64, t::Float64, mu::Float64, l::Int64)
@ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/xLH9y/src/arnoldi.jl:246
[12] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.LinearExponentialCache{…}, repeat_step::Bool)
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/tAI61/src/perform_step/linear_perform_step.jl:801
[13] perform_step!
@ ~/.julia/packages/OrdinaryDiffEq/tAI61/src/perform_step/linear_perform_step.jl:790 [inlined]
[14] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/tAI61/src/solve.jl:553
[15] __solve(::ODEProblem{…}, ::LinearExponential; kwargs::@Kwargs{})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/tAI61/src/solve.jl:7
[16] __solve
@ ~/.julia/packages/OrdinaryDiffEq/tAI61/src/solve.jl:1 [inlined]
[17] #solve_call#44
@ ~/.julia/packages/DiffEqBase/frOsk/src/solve.jl:612 [inlined]
[18] solve_call
@ ~/.julia/packages/DiffEqBase/frOsk/src/solve.jl:569 [inlined]
[19] #solve_up#53
@ ~/.julia/packages/DiffEqBase/frOsk/src/solve.jl:1092 [inlined]
[20] solve_up
@ ~/.julia/packages/DiffEqBase/frOsk/src/solve.jl:1078 [inlined]
[21] solve(prob::ODEProblem{…}, args::LinearExponential; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
@ DiffEqBase ~/.julia/packages/DiffEqBase/frOsk/src/solve.jl:1015
[22] solve(prob::ODEProblem{…}, args::LinearExponential)
@ DiffEqBase ~/.julia/packages/DiffEqBase/frOsk/src/solve.jl:1005
[23] top-level scope
@ ~/CN/test_linear-exponential_solver.jl:12while solve(prob, LinearExponential(krylov=:adaptive)) yields again
ERROR: Scalar indexing is disallowed.
Invocation of getindex resulted in scalar indexing of a GPU array.
This is typically caused by calling an iterating implementation of a method.
Such implementations *do not* execute on the GPU, but very slowly on the CPU,
and therefore should be avoided.
If you want to allow scalar iteration, use `allowscalar` or `@allowscalar`
to enable scalar iteration globally or for the operations in question.
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] errorscalar(op::String)
@ GPUArraysCore ~/.julia/packages/GPUArraysCore/GMsgk/src/GPUArraysCore.jl:155
[3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
@ GPUArraysCore ~/.julia/packages/GPUArraysCore/GMsgk/src/GPUArraysCore.jl:128
[4] assertscalar(op::String)
@ GPUArraysCore ~/.julia/packages/GPUArraysCore/GMsgk/src/GPUArraysCore.jl:116
[5] getindex
@ ~/.julia/packages/GPUArrays/qt4ax/src/host/indexing.jl:50 [inlined]
[6] copyto_unaliased!(deststyle::IndexLinear, dest::SubArray{…}, srcstyle::IndexLinear, src::CuArray{…})
@ Base ./abstractarray.jl:1088
[7] copyto!
@ ./abstractarray.jl:1068 [inlined]
[8] phiv_timestep!(U::CuArray{…}, ts::Vector{…}, A::CUDA.CUSPARSE.CuSparseMatrixCSC{…}, B::CuArray{…}; tau::Float64, m::Int64, tol::Float32, opnorm::typeof(LinearAlgebra.opnorm), iop::Int64, correct::Bool, caches::Tuple{…}, adaptive::Bool, delta::Float64, ishermitian::Bool, gamma::Float64, NA::Int64, verbose::Bool)
@ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/xLH9y/src/krylov_phiv_adaptive.jl:170
[9] #expv_timestep!#39
@ ~/.julia/packages/ExponentialUtilities/xLH9y/src/krylov_phiv_adaptive.jl:54 [inlined]
[10] expv_timestep!
@ ~/.julia/packages/ExponentialUtilities/xLH9y/src/krylov_phiv_adaptive.jl:51 [inlined]
[11] #expv_timestep!#38
@ ~/.julia/packages/ExponentialUtilities/xLH9y/src/krylov_phiv_adaptive.jl:48 [inlined]
[12] expv_timestep!
@ ~/.julia/packages/ExponentialUtilities/xLH9y/src/krylov_phiv_adaptive.jl:46 [inlined]
[13] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.LinearExponentialCache{…}, repeat_step::Bool)
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/tAI61/src/perform_step/linear_perform_step.jl:805
[14] perform_step!
@ ~/.julia/packages/OrdinaryDiffEq/tAI61/src/perform_step/linear_perform_step.jl:790 [inlined]
[15] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/tAI61/src/solve.jl:553
[16] __solve(::ODEProblem{…}, ::LinearExponential; kwargs::@Kwargs{})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/tAI61/src/solve.jl:7
[17] __solve
@ ~/.julia/packages/OrdinaryDiffEq/tAI61/src/solve.jl:1 [inlined]
[18] #solve_call#44
@ ~/.julia/packages/DiffEqBase/frOsk/src/solve.jl:612 [inlined]
[19] solve_call
@ ~/.julia/packages/DiffEqBase/frOsk/src/solve.jl:569 [inlined]
[20] #solve_up#53
@ ~/.julia/packages/DiffEqBase/frOsk/src/solve.jl:1092 [inlined]
[21] solve_up
@ ~/.julia/packages/DiffEqBase/frOsk/src/solve.jl:1078 [inlined]
[22] solve(prob::ODEProblem{…}, args::LinearExponential; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
@ DiffEqBase ~/.julia/packages/DiffEqBase/frOsk/src/solve.jl:1015
[23] solve(prob::ODEProblem{…}, args::LinearExponential)
@ DiffEqBase ~/.julia/packages/DiffEqBase/frOsk/src/solve.jl:1005
[24] top-level scope
@ ~/CN/test_linear-exponential_solver.jl:12
Some type information was truncated. Use `show(err)` to see complete types.Okay so it looks like the root of the issue is that phiv! needs to support this better.
ChrisRackauckas
commented
1 week ago Let's narrow this down to something that's just ExponentialUtilties.jl? There are some basic expv tests:
https://github.com/SciML/ExponentialUtilities.jl/blob/master/test/gpu/gputests.jl#L40-L57
so something must be missed by the tests.
mcarmesin
commented
5 days ago At least for krylov=:simple, the problem is that the V component of the KrylovSubspace is not constructed on the GPU. We have to alter the function alg_cache(alg::LinearExponential, …) in order to take care of constructing the KrylovSubspace in the right way, see my PR [https://github.com/SciML/OrdinaryDiffEq.jl/pull/2538]().
Describe the bug 🐞
The LinearExponential() fails, if the ODEProblem is based on a MatrixOperator with a CuSparseMatrixCSC.
Expected behavior
The ODE should be successfully solved. Therefor the solver shouldn't compute the full matrix exponential of the given matrix operator.
Minimal Reproducible Example 👇
Error & Stacktrace ⚠️ The solver computes directly the matrix exponential, which fails for a CuSparseMatrixCSC. See my bug report at Exponentiation utilities.
Environment (please complete the following information):
using Pkg; Pkg.status()using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)versioninfo()