Scalar indexing when using CuSparseMatrixCSR

[albertomercurio]

Hello, the following example seems to work

N = 1000
A = sprand(N, N, 1/N) + Diagonal(rand(N))
x = rand(N)

A_gpu = CUSPARSE.CuSparseMatrixCSC(A)
x_gpu = CuVector(x)

Pr_gpu = KrylovPreconditioners.kp_ilu0(A_gpu)

gmres(A_gpu, x_gpu, N=Pr_gpu, ldiv=true)

But not when I using CuSparseMatrixCSR, where I get the scalar indexing error.

[amontoison]

I am checking that now!

[amontoison]

I don' have any issue when I use the CSR format:

N = 1000
A = sprand(N, N, 1/N) + Diagonal(rand(N))
b = rand(N)

A_gpu = CUSPARSE.CuSparseMatrixCSR(A)
b_gpu = CuVector(b)

Pr_gpu = KrylovPreconditioners.kp_ilu0(A_gpu)

x_gpu, stats = gmres(A_gpu, b_gpu, N=Pr_gpu, ldiv=true)

[albertomercurio]

I get instead the following error

LoadError: 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/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/uiVyU/src/host/indexing.jl:50
  [6] getindex(A::CUDA.CUSPARSE.CuSparseMatrixCSR{Float64, Int32}, i0::Int64, i1::Int64)
    @ CUDA.CUSPARSE ~/.julia/packages/CUDA/1kIOw/lib/cusparse/array.jl:398
  [7] getindex
    @ ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:280 [inlined]
  [8] getindex
    @ ./subarray.jl:320 [inlined]
  [9] iterate
    @ ./abstractarray.jl:1209 [inlined]
 [10] iterate
    @ ./abstractarray.jl:1207 [inlined]
 [11] _all
    @ ./reduce.jl:1309 [inlined]
 [12] all
    @ ./reducedim.jl:995 [inlined]
 [13] _istriu
    @ ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1317 [inlined]
 [14] istriu(A::UnitLowerTriangular{Float64, CUDA.CUSPARSE.CuSparseMatrixCSR{Float64, Int32}}, k::Int64)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1310
 [15] istriu
    @ ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1309 [inlined]
 [16] ldiv!
    @ ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:972 [inlined]
 [17] ldiv!
    @ ~/.julia/packages/KrylovPreconditioners/dc3YZ/ext/CUDA/ilu0.jl:78 [inlined]
 [18] macro expansion
    @ ~/.julia/packages/KrylovPreconditioners/dc3YZ/ext/CUDA/ilu0.jl:86 [inlined]
 [19] macro expansion
    @ ./timing.jl:421 [inlined]
 [20] ldiv!(y::CuArray{Float64, 1, CUDA.DeviceMemory}, ilu::KrylovPreconditionersCUDAExt.NVIDIA_ILU0{CUDA.CUSPARSE.CuSparseMatrixCSR{Float64, Int32}}, x::CuArray{Float64, 1, CUDA.DeviceMemory})
    @ KrylovPreconditionersCUDAExt ~/.julia/packages/KrylovPreconditioners/dc3YZ/ext/CUDA/ilu0.jl:85
 [21] mulorldiv!
    @ ~/.julia/packages/Krylov/wgV9p/src/krylov_utils.jl:307 [inlined]
 [22] gmres!(solver::GmresSolver{Float64, Float64, CuArray{Float64, 1, CUDA.DeviceMemory}}, A::CUDA.CUSPARSE.CuSparseMatrixCSR{Float64, Int32}, b::CuArray{Float64, 1, CUDA.DeviceMemory}; M::UniformScaling{Bool}, N::KrylovPreconditionersCUDAExt.NVIDIA_ILU0{CUDA.CUSPARSE.CuSparseMatrixCSR{Float64, Int32}}, ldiv::Bool, restart::Bool, reorthogonalization::Bool, atol::Float64, rtol::Float64, itmax::Int64, timemax::Float64, verbose::Int64, history::Bool, callback::Krylov.var"#977#982", iostream::Core.CoreSTDOUT)
    @ Krylov ~/.julia/packages/Krylov/wgV9p/src/gmres.jl:241
 [23] gmres!
    @ ~/.julia/packages/Krylov/wgV9p/src/gmres.jl:108 [inlined]
 [24] #gmres#975
    @ ~/.julia/packages/Krylov/wgV9p/src/krylov_solve.jl:87 [inlined]
 [25] top-level scope
    @ ~/GitHub/Research/Undef/LinearSolve GPU/linearsolve_GPU.jl:110
 [26] eval
    @ ./boot.jl:430 [inlined]
 [27] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:2734
 [28] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::@Kwargs{})
    @ Base ./essentials.jl:1055
 [29] invokelatest(::Any, ::Any, ::Vararg{Any})
    @ Base ./essentials.jl:1052
 [30] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
    @ VSCodeServer ~/.vscode-server/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/eval.jl:271
 [31] (::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode-server/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/eval.jl:181
 [32] withpath(f::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
    @ VSCodeServer ~/.vscode-server/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/repl.jl:276
 [33] (::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode-server/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/eval.jl:179
 [34] hideprompt(f::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
    @ VSCodeServer ~/.vscode-server/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/repl.jl:38
 [35] #67
    @ ~/.vscode-server/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/eval.jl:150 [inlined]
 [36] with_logstate(f::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, logstate::Base.CoreLogging.LogState)
    @ Base.CoreLogging ./logging/logging.jl:522
 [37] with_logger
    @ ./logging/logging.jl:632 [inlined]
 [38] (::VSCodeServer.var"#66#71"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode-server/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/eval.jl:263
 [39] #invokelatest#2
    @ ./essentials.jl:1055 [inlined]
 [40] invokelatest(::Any)
    @ Base ./essentials.jl:1052
 [41] (::VSCodeServer.var"#64#65")()
    @ VSCodeServer ~/.vscode-server/extensions/julialang.language-julia-1.127.2/scripts/packages/VSCodeServer/src/eval.jl:34
in expression starting at /home/alberto/GitHub/Research/Undef/LinearSolve GPU/linearsolve_GPU.jl:110

This is my environment

Julia Version 1.11.3
Commit d63adeda50d (2025-01-21 19:42 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 32 × 13th Gen Intel(R) Core(TM) i9-13900KF
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, alderlake)
Threads: 16 default, 0 interactive, 8 GC (on 32 virtual cores)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 16

And CUDA info

CUDA runtime 12.6, artifact installation
CUDA driver 12.4
NVIDIA driver 550.120.0

CUDA libraries: 
- CUBLAS: 12.6.4
- CURAND: 10.3.7
- CUFFT: 11.3.0
- CUSOLVER: 11.7.1
- CUSPARSE: 12.5.4
- CUPTI: 2024.3.2 (API 24.0.0)
- NVML: 12.0.0+550.120

Julia packages: 
- CUDA: 5.6.1
- CUDA_Driver_jll: 0.10.4+0
- CUDA_Runtime_jll: 0.15.5+0

Toolchain:
- Julia: 1.11.3
- LLVM: 16.0.6

1 device:
  0: NVIDIA GeForce RTX 4090 (sm_89, 23.181 GiB / 23.988 GiB available)

[amontoison]

@albertomercurio What is your version of GPUArrays.jl? Is it v11 or not?

[albertomercurio]

I have version v11

(LinearSolve GPU) pkg> st --manifest GPUArrays
Status `~/GitHub/Research/Undef/LinearSolve GPU/Manifest.toml`
  [0c68f7d7] GPUArrays v11.2.2

[albertomercurio]

Here also the environment I'm using

(LinearSolve GPU) pkg> st
Status `~/GitHub/Research/Undef/LinearSolve GPU/Project.toml`
  [6e4b80f9] BenchmarkTools v1.6.0
  [052768ef] CUDA v5.6.1
  [ba0b0d4f] Krylov v0.9.10
  [45d422c2] KrylovPreconditioners v0.3.1
  [7ed4a6bd] LinearSolve v3.3.1

[albertomercurio]

I noticed that if I use the master branch of CUDA.jl, everything works. Let's wait the next stable release and see if this fixes it.

[amontoison]

It seems related to the issue with istriu (JuliaGPU/CUDA.jl#2585) but it was only related to Julia 1.11.2 (to the best of my knowledge).
Keep me updated if the next stable release fixes the problem.

[albertomercurio]

The new version of CUDA.jl fixes the problem. We can close this issue.