Scalar indexing when `mul!` into `SubArray`

[jonas-schulze]

Describe the bug

I can't mul!(::SubArray, ...).

To reproduce

The Minimal Working Example (MWE) for this bug:

using CUDA, CUDA.CUSPARSE
using LinearAlgebra
using SparseArrays

n = 100
A = sprandn(n, n, 2/n)
A = cu(A)

X = CUDA.randn(n, 2)
Y = CUDA.zeros(n, 2)
Yv = @view Y[1:n, 1:2]

@assert Yv isa SubArray
mul!(Yv, A, X) # boom
mul!(Yv, I, X) # only fails on Julia 1.12
mul!(Yv, 5, X) # only fails on Julia 1.12
# See https://github.com/JuliaGPU/CUDA.jl/issues/3041#issuecomment-3996816271

Project.zip

Expected behavior

It should just work.

Version info

Julia 1.10.9

Julia Version 1.10.9
Commit 5595d20a287 (2025-03-10 12:51 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × 12th Gen Intel(R) Core(TM) i5-12600K
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, alderlake)
Threads: 1 default, 0 interactive, 1 GC (on 16 virtual cores)
Environment:
  JULIA_PROJECT = @.

CUDA toolchain: 
- runtime 13.0, artifact installation
- driver 580.126.9 for 13.1
- compiler 13.1

CUDA libraries: 
- CUBLAS: 13.1.0
- CURAND: 10.4.0
- CUFFT: 12.0.0
- CUSOLVER: 12.0.4
- CUSPARSE: 12.6.3
- CUPTI: 2025.3.1 (API 13.0.1)
- NVML: 13.0.0+580.126.9

Julia packages: 
- CUDA: 5.9.7
- GPUArrays: 11.4.1
- GPUCompiler: 1.8.2
- KernelAbstractions: 0.9.40
- CUDA_Driver_jll: 13.1.0+2
- CUDA_Compiler_jll: 0.4.1+1
- CUDA_Runtime_jll: 0.19.2+0

Toolchain:
- Julia: 1.10.9
- LLVM: 15.0.7

1 device:
  0: NVIDIA RTX A4000 (sm_86, 15.359 GiB / 15.992 GiB available)

Julia 1.12.5

Julia Version 1.12.5
Commit 5fe89b8ddc1 (2026-02-09 16:05 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × 12th Gen Intel(R) Core(TM) i5-12600K
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, alderlake)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 16 virtual cores)
Environment:
  JULIA_PROJECT = @.

CUDA toolchain: 
- runtime 13.0, artifact installation
- driver 580.126.9 for 13.1
- compiler 13.1

CUDA libraries: 
- CUBLAS: 13.1.0
- CURAND: 10.4.0
- CUFFT: 12.0.0
- CUSOLVER: 12.0.4
- CUSPARSE: 12.6.3
- CUPTI: 2025.3.1 (API 13.0.1)
- NVML: 13.0.0+580.126.9

Julia packages: 
- CUDA: 5.9.6
- GPUArrays: 11.4.1
- GPUCompiler: 1.8.2
- KernelAbstractions: 0.9.40
- CUDA_Driver_jll: 13.1.0+2
- CUDA_Compiler_jll: 0.4.1+1
- CUDA_Runtime_jll: 0.19.2+0

Toolchain:
- Julia: 1.12.5
- LLVM: 18.1.7

1 device:
  0: NVIDIA RTX A4000 (sm_86, 15.359 GiB / 15.992 GiB available)

[kshyatt]

I think this is a problem related to JuliaGPU/AMDGPU.jl#875 (comment) and in fact it's occurring in GPUArrays.jl. I'll try to get a fix today.

[jonas-schulze]

$ julia +release mwe.jl 
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 ~/.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
    @ ~/.julia/packages/GPUArrays/3a5jB/src/host/indexing.jl:50 [inlined]
  [6] scalar_getindex
    @ ~/.julia/packages/GPUArrays/3a5jB/src/host/indexing.jl:36 [inlined]
  [7] _getindex
    @ ~/.julia/packages/GPUArrays/3a5jB/src/host/indexing.jl:19 [inlined]
  [8] getindex
    @ ~/.julia/packages/GPUArrays/3a5jB/src/host/indexing.jl:17 [inlined]
  [9] macro expansion
    @ /scratch/jschulze/dotjulia/juliaup/julia-1.12.5+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/generic.jl:100 [inlined]
 [10] _generic_matmatmul_nonadjtrans!(C::SubArray{Float32, 2, CuArray{Float32, 2, CUDA.DeviceMemory}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, A::CuSparseMatrixCSC{Float32, Int32}, B::CuArray{Float32, 2, CUDA.DeviceMemory}, alpha::Bool, beta::Bool)
    @ LinearAlgebra /scratch/jschulze/dotjulia/juliaup/julia-1.12.5+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:1033
 [11] __generic_matmatmul!
    @ /scratch/jschulze/dotjulia/juliaup/julia-1.12.5+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:1025 [inlined]
 [12] _generic_matmatmul!(C::SubArray{Float32, 2, CuArray{Float32, 2, CUDA.DeviceMemory}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, A::CuSparseMatrixCSC{Float32, Int32}, B::CuArray{Float32, 2, CUDA.DeviceMemory}, alpha::Bool, beta::Bool)
    @ LinearAlgebra /scratch/jschulze/dotjulia/juliaup/julia-1.12.5+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:1019
 [13] generic_matmatmul!
    @ /scratch/jschulze/dotjulia/juliaup/julia-1.12.5+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:1009 [inlined]
 [14] generic_matmatmul_wrapper!
    @ /scratch/jschulze/dotjulia/juliaup/julia-1.12.5+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:343 [inlined]
 [15] _mul!
    @ /scratch/jschulze/dotjulia/juliaup/julia-1.12.5+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:328 [inlined]
 [16] mul!
    @ /scratch/jschulze/dotjulia/juliaup/julia-1.12.5+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:297 [inlined]
 [17] mul!(C::SubArray{Float32, 2, CuArray{Float32, 2, CUDA.DeviceMemory}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, A::CuSparseMatrixCSC{Float32, Int32}, B::CuArray{Float32, 2, CUDA.DeviceMemory})
    @ LinearAlgebra /scratch/jschulze/dotjulia/juliaup/julia-1.12.5+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:265
 [18] top-level scope
    @ ~/git/obsidian-vault/0 Daily Notes/CUDA.jl#xyz mul! into SubArray/mwe.jl:14
 [19] include(mod::Module, _path::String)
    @ Base ./Base.jl:306
 [20] exec_options(opts::Base.JLOptions)
    @ Base ./client.jl:317
 [21] _start()
    @ Base ./client.jl:550
in expression starting at /home/jschulze/git/obsidian-vault/0 Daily Notes/CUDA.jl#xyz mul! into SubArray/mwe.jl:14

$ julia +lts mwe.jl 
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/3a5jB/src/host/indexing.jl:50
  [6] getindex(A::CuSparseMatrixCSC{Float32, Int32}, i0::Int64, i1::Int64)
    @ CUDA.CUSPARSE ~/.julia/packages/CUDA/YUHHP/lib/cusparse/array.jl:472
  [7] _generic_matmatmul!(C::SubArray{Float32, 2, CuArray{Float32, 2, CUDA.DeviceMemory}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, tA::Char, tB::Char, A::CuSparseMatrixCSC{Float32, Int32}, B::CuArray{Float32, 2, CUDA.DeviceMemory}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
    @ LinearAlgebra /scratch/jschulze/dotjulia/juliaup/julia-1.10.10+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:816
  [8] generic_matmatmul!(C::SubArray{Float32, 2, CuArray{Float32, 2, CUDA.DeviceMemory}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, tA::Char, tB::Char, A::CuSparseMatrixCSC{Float32, Int32}, B::CuArray{Float32, 2, CUDA.DeviceMemory}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
    @ LinearAlgebra /scratch/jschulze/dotjulia/juliaup/julia-1.10.10+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:785
  [9] mul!
    @ /scratch/jschulze/dotjulia/juliaup/julia-1.10.10+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:263 [inlined]
 [10] mul!(C::SubArray{Float32, 2, CuArray{Float32, 2, CUDA.DeviceMemory}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, A::CuSparseMatrixCSC{Float32, Int32}, B::CuArray{Float32, 2, CUDA.DeviceMemory})
    @ LinearAlgebra /scratch/jschulze/dotjulia/juliaup/julia-1.10.10+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:237
 [11] top-level scope
    @ ~/git/obsidian-vault/0 Daily Notes/CUDA.jl#xyz mul! into SubArray/mwe.jl:14
in expression starting at /home/jschulze/git/obsidian-vault/0 Daily Notes/CUDA.jl#xyz mul! into SubArray/mwe.jl:14

[kshyatt]

Indeed I think it's dispatching to the wrong matmatmul because GPUArrays doesn't really have a StridedCuMatrix equivalent.

[jonas-schulze]

Wow, I did not expect such a quick response. Thanks! 🚀

I have another, related issue for which I pirated defined the following. Should I open another issue?

# Adapted from uniformscaling.jl:287 as of Julia 1.10.9
# TODO: Discuss this change upstream:
function LinearAlgebra.mul!(out::SubArray{T, 2, <: CuMatrix{T}}, a::Number, B::UniformScaling, α::Number, β::Number) where {T}
    checksquare(out)
    if iszero(β)  # zero contribution of the out matrix
        fill!(out, zero(T))
    elseif !isone(β)
        rmul!(out, β)
    end
    s = convert(T, a*B.λ*α)
    if !iszero(s)
        ids = diagind(out)
        out[ids] .+= s # allocates :-(
    end
    return out
end

[kshyatt]

IMO yes, please open a separate issue for that so this one stays focused

[kshyatt]

Wow, I did not expect such a quick response. Thanks! 🚀

Same timezone + it's blocking my PR too 😉

[jonas-schulze]

Can I sneek mul!(Yv, I, X) and mul!(Yv, 5, X) into this issue? This popped up while creating a separate issue from my previous comment. Interestingly, these only fail on 1.12 but not on 1.10.

[jonas-schulze]

CC @vchuravy ☝️