Selectdim returns view - CuSparseMatrixCSC defaults to scalar operations on views

[AaronGhost]

Describe the bug

• The selectdim function systematically returns a SubArray when dense view calls return a CuArray.
  This leads to a discrepency in behaviours when writing view(X, :, 2) and selectdim(X, 2, 2).

• Multiplying a dense SubArray and a CuSparseMatrixCSC gets dispatched to _generic_matmatmul! and scalar indexing.

To reproduce

using CUDA
using SparseArrays
using CUDA.CUSPARSE

spm = sparse([1, 2, 3], [1, 2, 3], [1.0f0, 2.0f0, 3.0f0], 3, 3)
spm_d = CuSparseMatrixCSC(spm)

x_d = CUDA.rand(3, 50)
# Works fine with CuArray
spm_d * x_d
# Works fine with view (as view returns a CuArray)
spm_d * view(x_d, :, 1:20)
# Fails with selectdim (as selectdim returns a SubArray)
spm_d * selectdim(x_d, 2, 1:20)

Expected behavior

I would expect view(X, :, 2) and selectdim(X, 2, 2) to return the same object as the documentation for selectdim states that the results should be equivalent with a view. CuSparseMatrixCSC supporting sub-arrays would also be great but I understand that there is no easy way to dispatch on "dense" arrays.

Version info

Details on Julia:

Julia Version 1.12.2
Commit ca9b6662be (2025-11-20 16:25 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 16 × 11th Gen Intel(R) Core(TM) i9-11900K @ 3.50GHz
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, rocketlake)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 16 virtual cores)
Environment:
  JULIA_PKG_USE_CLI_GIT = true

Details on CUDA:

CUDA toolchain: 
- runtime 12.9, artifact installation
- driver 573.44.0 for 12.8
- compiler 12.9

CUDA libraries:
- CUBLAS: 12.9.1
- CURAND: 10.3.10
- CUFFT: 11.4.1
- CUSOLVER: 11.7.5
- CUSPARSE: 12.5.10
- CUPTI: 2025.2.1 (API 12.9.1)
- NVML: 12.0.0+573.44

Julia packages:
- CUDA: 5.9.5
- CUDA_Driver_jll: 13.0.2+0
- CUDA_Compiler_jll: 0.3.0+0
- CUDA_Runtime_jll: 0.19.2+0

Toolchain:
- Julia: 1.12.2
- LLVM: 18.1.7

1 device:
  0: NVIDIA RTX A5000 (sm_86, 21.120 GiB / 23.988 GiB available)

Additional context

Error trace:

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:44
  [2] errorscalar(op::String)
    @ GPUArraysCore C:\Users\removed\.julia\packages\GPUArraysCore\aNaXo\src\GPUArraysCore.jl:151
  [3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
    @ GPUArraysCore C:\Users\removed\.julia\packages\GPUArraysCore\aNaXo\src\GPUArraysCore.jl:124
  [4] assertscalar(op::String)
    @ GPUArraysCore C:\Users\removed\.julia\packages\GPUArraysCore\aNaXo\src\GPUArraysCore.jl:112
  [5] getindex
    @ C:\Users\removed\.julia\packages\GPUArrays\0F4Dn\src\host\indexing.jl:50 [inlined]
  [6] scalar_getindex
    @ C:\Users\removed\.julia\packages\GPUArrays\0F4Dn\src\host\indexing.jl:36 [inlined]
  [7] _getindex
    @ C:\Users\removed\.julia\packages\GPUArrays\0F4Dn\src\host\indexing.jl:19 [inlined]
  [8] getindex
    @ C:\Users\removed\.julia\packages\GPUArrays\0F4Dn\src\host\indexing.jl:17 [inlined]
  [9] getindex
    @ .\subarray.jl:316 [inlined]
 [10] macro expansion
    @ C:\Users\removed\.julia\juliaup\julia-1.12.2+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\generic.jl:100 [inlined]
 [11] _generic_matmatmul_nonadjtrans!(C::CuArray{Float32, 2, CUDA.DeviceMemory}, A::CuSparseMatrixCSC{Float32, Int32}, B::SubArray{Float32, 2, CuArray{…}, Tuple{…}, true}, alpha::Bool, beta::Bool)
    @ LinearAlgebra C:\Users\removed\.julia\juliaup\julia-1.12.2+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\matmul.jl:1035
 [12] __generic_matmatmul!
    @ C:\Users\removed\.julia\juliaup\julia-1.12.2+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\matmul.jl:1027 [inlined]
 [13] _generic_matmatmul!(C::CuArray{Float32, 2, CUDA.DeviceMemory}, A::CuSparseMatrixCSC{Float32, Int32}, B::SubArray{Float32, 2, CuArray{…}, Tuple{…}, true}, alpha::Bool, beta::Bool)
    @ LinearAlgebra C:\Users\removed\.julia\juliaup\julia-1.12.2+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\matmul.jl:1021
 [14] generic_matmatmul!
    @ C:\Users\removed\.julia\juliaup\julia-1.12.2+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\matmul.jl:1011 [inlined]
 [15] generic_matmatmul_wrapper!
    @ C:\Users\removed\.julia\juliaup\julia-1.12.2+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\matmul.jl:343 [inlined]
 [16] _mul!
    @ C:\Users\removed\.julia\juliaup\julia-1.12.2+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\matmul.jl:328 [inlined]
 [17] mul!
    @ C:\Users\removed\.julia\juliaup\julia-1.12.2+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\matmul.jl:297 [inlined]
 [18] mul!
    @ C:\Users\removed\.julia\juliaup\julia-1.12.2+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\matmul.jl:265 [inlined]
 [19] mul
    @ C:\Users\removed\.julia\juliaup\julia-1.12.2+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\matmul.jl:118 [inlined]
 [20] *(A::CuSparseMatrixCSC{Float32, Int32}, B::SubArray{Float32, 2, CuArray{Float32, 2, CUDA.DeviceMemory}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true})
    @ LinearAlgebra C:\Users\removed\.julia\juliaup\julia-1.12.2+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\matmul.jl:114
 [21] top-level scope
    @ D:\Documents\Julia\CompressedSense.jl\CUDA_bug.jl:14
 [22] include(mapexpr::Function, mod::Module, _path::String)
    @ Base .\Base.jl:307
 [23] top-level scope
    @ REPL[3]:1