Default algs for CUSPARSE failing

[hexaeder]

Currently, the default algorithm choice is not really able to handle sparse cuda matrices. In short:

• CuSparseMatrixCSC -> defaults always to GenericLU which is not GPU compatible (scalar indexing)
• CuSparseMatrixCSR
  • CUDSS not available -> defaults to Krylov which not GPU compatibe (scalar indexing). It gives a warning however
  • CUDSS available + asym matrix -> works
  • CUDSS available + symmetric matrix -> Cholesky which is not GPU compatible (scalar indexing)

I am not sure what the right thing here is. Maybe just give a clearer error message? Or maybe convert to dense?

Code to play around with:

using LinearSolve
using SparseArrays
using CUDA
using CUDSS

b = Float64[1, 2, 3, 4]
b_gpu = adapt(CuArray, b)

A = Float64[1 1 0 0
            0 1 1 0
            0 0 3 1
            0 0 0 4]
A_sym = Float64[1 1 0 0
                1 0 0 2
                0 0 3 0
                0 2 0 0]

# asymetric matrix on GPU (dense)
A_gpu = adapt(CuArray, A)
prob_gpu = LinearProblem(A_gpu, b_gpu)
solve(prob_gpu) # works

# asymetric sparse matrix
A_gpu_sparse1 = CUDA.CUSPARSE.CuSparseMatrixCSR(sparse(A))
prob_gpu_sparse1 = LinearProblem(A_gpu_sparse1, b_gpu)
solve(prob_gpu_sparse1) # works

A_gpu_sparse2 = CUDA.CUSPARSE.CuSparseMatrixCSC(sparse(A))
prob_gpu_sparse2 = LinearProblem(A_gpu_sparse2, b_gpu)
solve(prob_gpu_sparse2) # fails because of generic LU


# symmetric matrix on GPU (dense)
A_gpu_sym = adapt(CuArray, A_sym)
prob_gpu_sym = LinearProblem(A_gpu_sym, b_gpu)
solve(prob_gpu_sym) # works

# symetric sparse on GPU
A_gpu_sym_sparse1 = CUDA.CUSPARSE.CuSparseMatrixCSR(sparse(A_sym))
prob_gpu_sym_sparse1 = LinearProblem(A_gpu_sym_sparse1, b_gpu)
solve(prob_gpu_sym_sparse1) # errors in colesky

A_gpu_sym_sparse2 = CUDA.CUSPARSE.CuSparseMatrixCSC(sparse(A_sym))
prob_gpu_sym_sparse2 = LinearProblem(A_gpu_sym_sparse2, b_gpu)
solve(prob_gpu_sym_sparse2) # errors in generic LU

[ChrisRackauckas]

Interesting... I thought CUDSS would work on these. If you force CUDSS is it fine?

[hexaeder]

Ah yes I missed that! Explicitly requsting LUFacotorization does not fall back to Cholesky.

A_gpu_sym_sparse1 = CUDA.CUSPARSE.CuSparseMatrixCSR(sparse(A_sym))
prob_gpu_sym_sparse1 = LinearProblem(A_gpu_sym_sparse1, b_gpu)
solve(prob_gpu_sym_sparse1, LUFactorization()) # works

[ChrisRackauckas]

Okay so the default just needs to be set.

[hexaeder]

unfortunatley, does not fix it. Your robot should have used the MWE :D

• CSC -> does not work at all, I think thats expected since its not supported by CUDSS, but would be nice if it errors with something else than scalar indexing error
• CSR -> still does not work on symmetric matrixes, only if you pass LUFactorization() explicitly. Otherwise it falls back to Cholesky and falls.

using LinearSolve
using SparseArrays
using CUDA
using CUDSS
using Adapt

b = Float64[1, 2, 3, 4]
b_gpu = adapt(CuArray, b)

A = Float64[1 1 0 0
            0 1 1 0
            0 0 3 1
            0 0 0 4]
A_sym = Float64[1 1 0 0
                1 0 0 2
                0 0 3 0
                0 2 0 0]

# asymetric sparse matrix
A_gpu_sparse1 = CUDA.CUSPARSE.CuSparseMatrixCSR(sparse(A))
prob_gpu_sparse1 = LinearProblem(A_gpu_sparse1, b_gpu)
solve(prob_gpu_sparse1) # works

A_gpu_sparse2 = CUDA.CUSPARSE.CuSparseMatrixCSC(sparse(A))
prob_gpu_sparse2 = LinearProblem(A_gpu_sparse2, b_gpu)
solve(prob_gpu_sparse2) # fails because of generic LU (scalar indexing)
solve(prob_gpu_sparse2, LUFactorization()) # also fails because of generic LU

# symetric sparse matrix
A_gpu_sym_sparse1 = CUDA.CUSPARSE.CuSparseMatrixCSR(sparse(A_sym))
prob_gpu_sym_sparse1 = LinearProblem(A_gpu_sym_sparse1, b_gpu)
solve(prob_gpu_sym_sparse1) # fails in colesky
solve(prob_gpu_sym_sparse1, LUFactorization()) # works, epxlicitly stating LU fixes cholesky fallback

A_gpu_sym_sparse2 = CUDA.CUSPARSE.CuSparseMatrixCSC(sparse(A_sym))
prob_gpu_sym_sparse2 = LinearProblem(A_gpu_sym_sparse2, b_gpu)
solve(prob_gpu_sym_sparse2) # errors in generic LU
solve(prob_gpu_sym_sparse2, LUFactorization()) # errors in generic LU

[ChrisRackauckas]

I don't have an NVIDIA GPU on me nor can it run it. I don't see why it would go to Cholesky if its symmetric in value but not Symmetric as a matrix type. The error messages are still useful though.

[hexaeder]

For context, here are the full stacktraces from the MWE shown above

Asymmetric + CSC Style, fails regardless of alg specification

✻ julia> A_gpu_sparse2 = CUDA.CUSPARSE.CuSparseMatrixCSC(sparse(A));

✻ julia> prob_gpu_sparse2 = LinearProblem(A_gpu_sparse2, b_gpu);

✻ julia> solve(prob_gpu_sparse2) # fails because of generic LU (scalar indexing)

⚠️ Stacktrace

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:44
  [2] errorscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:151
  [3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:124
  [4] assertscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:112
  [5] getindex(A::CuArray{Int32, 1, CUDA.DeviceMemory}, I::Int64)
    @ GPUArrays ~/.julia/packages/GPUArrays/w335n/src/host/indexing.jl:50
  [6] getindex(A::CuSparseMatrixCSC{Float64, Int32}, i0::Int64, i1::Int64)
    @ CUDA.CUSPARSE ~/.julia/packages/CUDA/wPaDI/lib/cusparse/array.jl:468
  [7] generic_lufact!(A::CuSparseMatrixCSC{Float64, Int32}, pivot::LinearAlgebra.RowMaximum; check::Bool, allowsingular::Bool)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:167
  [8] generic_lufact!
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:152 [inlined]
  [9] lu! (repeats 2 times)
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:151 [inlined]
 [10] _lu
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:347 [inlined]
 [11] lu
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:341 [inlined]
 [12] lu_instance
    @ ~/.julia/packages/ArrayInterface/U8kYm/src/ArrayInterface.jl:515 [inlined]
 [13] init_cacheval
    @ ~/.julia/dev/LinearSolve/ext/LinearSolveSparseArraysExt.jl:144 [inlined]
 [14] macro expansion
    @ ~/.julia/dev/LinearSolve/src/default.jl:479 [inlined]
 [15] init_cacheval(alg::LinearSolve.DefaultLinearSolver, A::CuSparseMatrixCSC{…}, b::CuArray{…}, u::CuArray{…}, Pl::IdentityOperator, Pr::IdentityOperator, maxiters::Int64, abstol::Float64, reltol::Float64, verbose::LinearVerbosity{…}, assump::OperatorAssumptions{…})
    @ LinearSolve ~/.julia/dev/LinearSolve/src/default.jl:462
 [16] __init(::LinearProblem{…}, ::LinearSolve.DefaultLinearSolver; alias::LinearAliasSpecifier, abstol::Float64, reltol::Float64, maxiters::Int64, verbose::LinearVerbosity{…}, Pl::Nothing, Pr::Nothing, assumptions::OperatorAssumptions{…}, sensealg::LinearSolveAdjoint{…}, kwargs::@Kwargs{})
    @ LinearSolve ~/.julia/dev/LinearSolve/src/common.jl:381
 [17] __init
    @ ~/.julia/dev/LinearSolve/src/common.jl:279 [inlined]
 [18] #init#15
    @ ~/.julia/dev/LinearSolve/src/common.jl:276 [inlined]
 [19] init
    @ ~/.julia/dev/LinearSolve/src/common.jl:275 [inlined]
 [20] #solve#20
    @ ~/.julia/dev/LinearSolve/src/common.jl:437 [inlined]
 [21] solve
    @ ~/.julia/dev/LinearSolve/src/common.jl:435 [inlined]
 [22] #solve#19
    @ ~/.julia/dev/LinearSolve/src/common.jl:432 [inlined]
 [23] solve
    @ ~/.julia/dev/LinearSolve/src/common.jl:430 [inlined]
 [24] #solve#18
    @ ~/.julia/dev/LinearSolve/src/common.jl:427 [inlined]
 [25] solve(::LinearProblem{Nothing, true, CuSparseMatrixCSC{Float64, Int32}, CuArray{Float64, 1, CUDA.DeviceMemory}, SciMLBase.NullParameters, Nothing, @Kwargs{}})
    @ LinearSolve ~/.julia/dev/LinearSolve/src/common.jl:426
 [26] top-level scope
    @ REPL[147]:1
 [27] top-level scope
    @ REPL:1
Some type information was truncated. Use `show(err)` to see complete types.

Same system, this time epxlicitly passing alg

✻ julia> solve(prob_gpu_sparse2, LUFactorization()) # also fails because of generic LU

⚠️ Stacktrace

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:44
  [2] errorscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:151
  [3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:124
  [4] assertscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:112
  [5] getindex(A::CuArray{Int32, 1, CUDA.DeviceMemory}, I::Int64)
    @ GPUArrays ~/.julia/packages/GPUArrays/w335n/src/host/indexing.jl:50
  [6] getindex(A::CuSparseMatrixCSC{Float64, Int32}, i0::Int64, i1::Int64)
    @ CUDA.CUSPARSE ~/.julia/packages/CUDA/wPaDI/lib/cusparse/array.jl:468
  [7] generic_lufact!(A::CuSparseMatrixCSC{Float64, Int32}, pivot::LinearAlgebra.RowMaximum; check::Bool, allowsingular::Bool)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:167
  [8] generic_lufact!
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:152 [inlined]
  [9] lu! (repeats 2 times)
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:151 [inlined]
 [10] _lu
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:347 [inlined]
 [11] lu
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:341 [inlined]
 [12] lu_instance
    @ ~/.julia/packages/ArrayInterface/U8kYm/src/ArrayInterface.jl:515 [inlined]
 [13] init_cacheval
    @ ~/.julia/dev/LinearSolve/ext/LinearSolveSparseArraysExt.jl:144 [inlined]
 [14] __init(::LinearProblem{…}, ::LUFactorization{…}; alias::LinearAliasSpecifier, abstol::Float64, reltol::Float64, maxiters::Int64, verbose::LinearVerbosity{…}, Pl::Nothing, Pr::Nothing, assumptions::OperatorAssumptions{…}, sensealg::LinearSolveAdjoint{…}, kwargs::@Kwargs{})
    @ LinearSolve ~/.julia/dev/LinearSolve/src/common.jl:381
 [15] __init
    @ ~/.julia/dev/LinearSolve/src/common.jl:279 [inlined]
 [16] #init#15
    @ ~/.julia/dev/LinearSolve/src/common.jl:276 [inlined]
 [17] init
    @ ~/.julia/dev/LinearSolve/src/common.jl:275 [inlined]
 [18] #solve#20
    @ ~/.julia/dev/LinearSolve/src/common.jl:437 [inlined]
 [19] solve(::LinearProblem{Nothing, true, CuSparseMatrixCSC{Float64, Int32}, CuArray{Float64, 1, CUDA.DeviceMemory}, SciMLBase.NullParameters, Nothing, @Kwargs{}}, ::LUFactorization{LinearAlgebra.RowMaximum})
    @ LinearSolve ~/.julia/dev/LinearSolve/src/common.jl:435
 [20] top-level scope
    @ REPL[148]:1
 [21] top-level scope
    @ REPL:1
Some type information was truncated. Use `show(err)` to see complete types.

Symmetric + CSR - Style, fails on default algo, passes on expicit LU

✻ julia> A_gpu_sym_sparse1 = CUDA.CUSPARSE.CuSparseMatrixCSR(sparse(A_sym))
4×4 CuSparseMatrixCSR{Float64, Int32} with 6 stored entries:
 1.0  1.0   ⋅    ⋅ 
 1.0   ⋅    ⋅   2.0
  ⋅    ⋅   3.0   ⋅ 
  ⋅   2.0   ⋅    ⋅ 

✻ julia> prob_gpu_sym_sparse1 = LinearProblem(A_gpu_sym_sparse1, b_gpu);

✻ julia> solve(prob_gpu_sym_sparse1) # fails in colesky

⚠️ Stacktrace

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:44
  [2] errorscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:151
  [3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:124
  [4] assertscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:112
  [5] getindex(A::CuArray{Int32, 1, CUDA.DeviceMemory}, I::Int64)
    @ GPUArrays ~/.julia/packages/GPUArrays/w335n/src/host/indexing.jl:50
  [6] getindex(A::CuSparseMatrixCSR{Float64, Int32}, i0::Int64, i1::Int64)
    @ CUDA.CUSPARSE ~/.julia/packages/CUDA/wPaDI/lib/cusparse/array.jl:478
  [7] _chol!(A::CuSparseMatrixCSR{Float64, Int32}, ::Type{LinearAlgebra.UpperTriangular})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/cholesky.jl:202
  [8] #cholesky!#163
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/cholesky.jl:417 [inlined]
  [9] cholesky!
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/cholesky.jl:416 [inlined]
 [10] cholesky!(A::CuSparseMatrixCSR{Float64, Int32}, ::LinearAlgebra.NoPivot; check::Bool)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/cholesky.jl:450
 [11] cholesky! (repeats 2 times)
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/cholesky.jl:444 [inlined]
 [12] _cholesky
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/cholesky.jl:552 [inlined]
 [13] cholesky
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/cholesky.jl:542 [inlined]
 [14] cholesky_instance
    @ ~/.julia/packages/ArrayInterface/U8kYm/src/ArrayInterface.jl:435 [inlined]
 [15] init_cacheval
    @ ~/.julia/dev/LinearSolve/src/factorization.jl:398 [inlined]
 [16] macro expansion
    @ ~/.julia/dev/LinearSolve/src/default.jl:479 [inlined]
 [17] init_cacheval(alg::LinearSolve.DefaultLinearSolver, A::CuSparseMatrixCSR{…}, b::CuArray{…}, u::CuArray{…}, Pl::IdentityOperator, Pr::IdentityOperator, maxiters::Int64, abstol::Float64, reltol::Float64, verbose::LinearVerbosity{…}, assump::OperatorAssumptions{…})
    @ LinearSolve ~/.julia/dev/LinearSolve/src/default.jl:462
 [18] __init(::LinearProblem{…}, ::LinearSolve.DefaultLinearSolver; alias::LinearAliasSpecifier, abstol::Float64, reltol::Float64, maxiters::Int64, verbose::LinearVerbosity{…}, Pl::Nothing, Pr::Nothing, assumptions::OperatorAssumptions{…}, sensealg::LinearSolveAdjoint{…}, kwargs::@Kwargs{})
    @ LinearSolve ~/.julia/dev/LinearSolve/src/common.jl:381
 [19] __init
    @ ~/.julia/dev/LinearSolve/src/common.jl:279 [inlined]
 [20] #init#15
    @ ~/.julia/dev/LinearSolve/src/common.jl:276 [inlined]
 [21] init
    @ ~/.julia/dev/LinearSolve/src/common.jl:275 [inlined]
 [22] #solve#20
    @ ~/.julia/dev/LinearSolve/src/common.jl:437 [inlined]
 [23] solve
    @ ~/.julia/dev/LinearSolve/src/common.jl:435 [inlined]
 [24] #solve#19
    @ ~/.julia/dev/LinearSolve/src/common.jl:432 [inlined]
 [25] solve
    @ ~/.julia/dev/LinearSolve/src/common.jl:430 [inlined]
 [26] #solve#18
    @ ~/.julia/dev/LinearSolve/src/common.jl:427 [inlined]
 [27] solve(::LinearProblem{Nothing, true, CuSparseMatrixCSR{Float64, Int32}, CuArray{Float64, 1, CUDA.DeviceMemory}, SciMLBase.NullParameters, Nothing, @Kwargs{}})
    @ LinearSolve ~/.julia/dev/LinearSolve/src/common.jl:426
 [28] top-level scope
    @ REPL[151]:1
 [29] top-level scope
    @ REPL:1
Some type information was truncated. Use `show(err)` to see complete types.

Works wnn passing LUFactorization as this seems to suppress the auto cholesky on sym values

✻ julia> solve(prob_gpu_sym_sparse1, LUFactorization()) # works, epxlicitly stating LU fixes cholesky fallback
retcode: Success
u: 4-element CuArray{Float64, 1, CUDA.DeviceMemory}:
 -0.999999999999925
  1.999999999999925
  1.0
  1.5010215292932116

Symmetric + CSC - Style, fails regardless of alg

✻ julia> A_gpu_sym_sparse2 = CUDA.CUSPARSE.CuSparseMatrixCSC(sparse(A_sym))
4×4 CuSparseMatrixCSC{Float64, Int32} with 6 stored entries:
 1.0  1.0   ⋅    ⋅
 1.0   ⋅    ⋅   2.0
  ⋅    ⋅   3.0   ⋅
  ⋅   2.0   ⋅    ⋅

✻ julia> prob_gpu_sym_sparse2 = LinearProblem(A_gpu_sym_sparse2, b_gpu);

✻ julia> solve(prob_gpu_sym_sparse2) # errors in generic LU

⚠️ Stacktrace

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:44
  [2] errorscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:151
  [3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:124
  [4] assertscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:112
  [5] getindex(A::CuArray{Int32, 1, CUDA.DeviceMemory}, I::Int64)
    @ GPUArrays ~/.julia/packages/GPUArrays/w335n/src/host/indexing.jl:50
  [6] getindex(A::CuSparseMatrixCSC{Float64, Int32}, i0::Int64, i1::Int64)
    @ CUDA.CUSPARSE ~/.julia/packages/CUDA/wPaDI/lib/cusparse/array.jl:468
  [7] generic_lufact!(A::CuSparseMatrixCSC{Float64, Int32}, pivot::LinearAlgebra.RowMaximum; check::Bool, allowsingular::Bool)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:167
  [8] generic_lufact!
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:152 [inlined]
  [9] lu! (repeats 2 times)
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:151 [inlined]
 [10] _lu
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:347 [inlined]
 [11] lu
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:341 [inlined]
 [12] lu_instance
    @ ~/.julia/packages/ArrayInterface/U8kYm/src/ArrayInterface.jl:515 [inlined]
 [13] init_cacheval
    @ ~/.julia/dev/LinearSolve/ext/LinearSolveSparseArraysExt.jl:144 [inlined]
 [14] macro expansion
    @ ~/.julia/dev/LinearSolve/src/default.jl:479 [inlined]
 [15] init_cacheval(alg::LinearSolve.DefaultLinearSolver, A::CuSparseMatrixCSC{…}, b::CuArray{…}, u::CuArray{…}, Pl::IdentityOperator, Pr::IdentityOperator, maxiters::Int64, abstol::Float64, reltol::Float64, verbose::LinearVerbosity{…}, assump::OperatorAssumptions{…})
    @ LinearSolve ~/.julia/dev/LinearSolve/src/default.jl:462
 [16] __init(::LinearProblem{…}, ::LinearSolve.DefaultLinearSolver; alias::LinearAliasSpecifier, abstol::Float64, reltol::Float64, maxiters::Int64, verbose::LinearVerbosity{…}, Pl::Nothing, Pr::Nothing, assumptions::OperatorAssumptions{…}, sensealg::LinearSolveAdjoint{…}, kwargs::@Kwargs{})
    @ LinearSolve ~/.julia/dev/LinearSolve/src/common.jl:381
 [17] __init
    @ ~/.julia/dev/LinearSolve/src/common.jl:279 [inlined]
 [18] #init#15
    @ ~/.julia/dev/LinearSolve/src/common.jl:276 [inlined]
 [19] init
    @ ~/.julia/dev/LinearSolve/src/common.jl:275 [inlined]
 [20] #solve#20
    @ ~/.julia/dev/LinearSolve/src/common.jl:437 [inlined]
 [21] solve
    @ ~/.julia/dev/LinearSolve/src/common.jl:435 [inlined]
 [22] #solve#19
    @ ~/.julia/dev/LinearSolve/src/common.jl:432 [inlined]
 [23] solve
    @ ~/.julia/dev/LinearSolve/src/common.jl:430 [inlined]
 [24] #solve#18
    @ ~/.julia/dev/LinearSolve/src/common.jl:427 [inlined]
 [25] solve(::LinearProblem{Nothing, true, CuSparseMatrixCSC{Float64, Int32}, CuArray{Float64, 1, CUDA.DeviceMemory}, SciMLBase.NullParameters, Nothing, @Kwargs{}})
    @ LinearSolve ~/.julia/dev/LinearSolve/src/common.jl:426
 [26] top-level scope
    @ REPL[155]:1
 [27] top-level scope
    @ REPL:1
Some type information was truncated. Use `show(err)` to see complete types.

problem persists when explicitly passing alg

✻ julia> solve(prob_gpu_sym_sparse2, LUFactorization()) # errors in generic LU

⚠️ Stacktrace

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:44
  [2] errorscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:151
  [3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:124
  [4] assertscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:112
  [5] getindex(A::CuArray{Int32, 1, CUDA.DeviceMemory}, I::Int64)
    @ GPUArrays ~/.julia/packages/GPUArrays/w335n/src/host/indexing.jl:50
  [6] getindex(A::CuSparseMatrixCSC{Float64, Int32}, i0::Int64, i1::Int64)
    @ CUDA.CUSPARSE ~/.julia/packages/CUDA/wPaDI/lib/cusparse/array.jl:468
  [7] generic_lufact!(A::CuSparseMatrixCSC{Float64, Int32}, pivot::LinearAlgebra.RowMaximum; check::Bool, allowsingular::Bool)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:167
  [8] generic_lufact!
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:152 [inlined]
  [9] lu! (repeats 2 times)
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:151 [inlined]
 [10] _lu
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:347 [inlined]
 [11] lu
    @ ~/.julia/juliaup/julia-1.12.1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:341 [inlined]
 [12] lu_instance
    @ ~/.julia/packages/ArrayInterface/U8kYm/src/ArrayInterface.jl:515 [inlined]
 [13] init_cacheval
    @ ~/.julia/dev/LinearSolve/ext/LinearSolveSparseArraysExt.jl:144 [inlined]
 [14] __init(::LinearProblem{…}, ::LUFactorization{…}; alias::LinearAliasSpecifier, abstol::Float64, reltol::Float64, maxiters::Int64, verbose::LinearVerbosity{…}, Pl::Nothing, Pr::Nothing, assumptions::OperatorAssumptions{…}, sensealg::LinearSolveAdjoint{…}, kwargs::@Kwargs{})
    @ LinearSolve ~/.julia/dev/LinearSolve/src/common.jl:381
 [15] __init
    @ ~/.julia/dev/LinearSolve/src/common.jl:279 [inlined]
 [16] #init#15
    @ ~/.julia/dev/LinearSolve/src/common.jl:276 [inlined]
 [17] init
    @ ~/.julia/dev/LinearSolve/src/common.jl:275 [inlined]
 [18] #solve#20
    @ ~/.julia/dev/LinearSolve/src/common.jl:437 [inlined]
 [19] solve(::LinearProblem{Nothing, true, CuSparseMatrixCSC{Float64, Int32}, CuArray{Float64, 1, CUDA.DeviceMemory}, SciMLBase.NullParameters, Nothing, @Kwargs{}}, ::LUFactorization{LinearAlgebra.RowMaximum})
    @ LinearSolve ~/.julia/dev/LinearSolve/src/common.jl:435
 [20] top-level scope
    @ REPL[156]:1
 [21] top-level scope
    @ REPL:1
Some type information was truncated. Use `show(err)` to see complete types.