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
using Test

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

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

# Have the Y matrix be larger than Yv, such that Yv is non-contiguous.
# Only then we get an object of type SubArray in the end.
# See https://github.com/JuliaGPU/GPUArrays.jl/pull/704
@assert Yv isa SubArray

@testset "mul!(::Subarray, _, _)" begin
    @test mul!(Yv, A, X) isa AbstractMatrix # scalar indexing

    # These fail on Julia 1.12 but are fine on Julia 1.10:
    @test mul!(Yv, I, X) isa AbstractMatrix # scalar indexing
    @test mul!(Yv, 5, X) isa AbstractMatrix # scalar indexing
    # See https://github.com/JuliaGPU/CUDA.jl/issues/3041#issuecomment-3996816271
end

Project.zip

Expected behavior

It should just work.

Version info

Julia 1.10.11

Julia Version 1.10.11
Commit a2b11907d7b (2026-03-09 14:59 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × Intel(R) Core(TM) Ultra 7 255H
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, goldmont)
Threads: 1 default, 0 interactive, 1 GC (on 16 virtual cores)
Environment:
  JULIA_CPU_TARGET = goldmont
  JULIA_PKG_PRECOMPILE_AUTO = 0
  JULIA_PROJECT = @.

CUDA toolchain: 
- runtime 13.3.0, artifact installation
- driver 595.84.0 for 13.3
- compiler 13.3.33, artifact installation

CUDA libraries: 
- cuBLAS: 13.6.0
- cuSPARSE: 12.8.2
- cuSOLVER: 12.2.6
- cuFFT: 12.3.0
- cuRAND: 10.4.3
- CUPTI: 2026.2.1 (API 13.3.1)
- NVML: 13.0.0+595.84

Julia packages: 
- CUDACore: 6.2.1
- GPUArrays: 11.5.10
- GPUCompiler: 1.23.0
- KernelAbstractions: 0.9.42
- CUDA_Driver_jll: 13.3.1+0
- CUDA_Compiler_jll: 0.4.4+1
- CUDA_Runtime_jll: 0.23.0+1
- NVPTX_LLVM_Backend_jll: 22.1.7+1

Toolchain:
- Julia: 1.10.11
- LLVM: 15.0.7

1 device:
  0: NVIDIA RTX PRO 500 Blackwell Generation Laptop GPU (sm_120, 5.525 GiB / 5.970 GiB available)
     compiles to sm_120a / PTX 9.3 (LLVM: sm_120a / PTX 9.0)

Julia 1.12.6

Julia Version 1.12.6
Commit 15346901f00 (2026-04-09 19:20 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × Intel(R) Core(TM) Ultra 7 255H
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, arrowlake)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 16 virtual cores)
Environment:
  JULIA_CPU_TARGET = goldmont
  JULIA_PKG_PRECOMPILE_AUTO = 0
  JULIA_PROJECT = @.

CUDA toolchain: 
- runtime 13.3.0, artifact installation
- driver 595.84.0 for 13.3
- compiler 13.3.33, artifact installation

CUDA libraries: 
- cuBLAS: 13.6.0
- cuSPARSE: 12.8.2
- cuSOLVER: 12.2.6
- cuFFT: 12.3.0
- cuRAND: 10.4.3
- CUPTI: 2026.2.1 (API 13.3.1)
- NVML: 13.0.0+595.84

Julia packages: 
- CUDACore: 6.2.1
- GPUArrays: 11.5.10
- GPUCompiler: 1.23.0
- KernelAbstractions: 0.9.42
- CUDA_Driver_jll: 13.3.1+0
- CUDA_Compiler_jll: 0.4.4+1
- CUDA_Runtime_jll: 0.23.0+1
- NVPTX_LLVM_Backend_jll: 22.1.7+1

Toolchain:
- Julia: 1.12.6
- LLVM: 18.1.7

1 device:
  0: NVIDIA RTX PRO 500 Blackwell Generation Laptop GPU (sm_120, 5.525 GiB / 5.970 GiB available)
     compiles to sm_120a / PTX 9.3 (LLVM: sm_120a / PTX 9.0)

[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 ☝️

[kshyatt]

Also can't reproduce the MWE failure here, Yv is no longer a SubArray. Feel free to reopen with another MWE if you run into this again.

[jonas-schulze]

After upgrading both environments to CUDA.jl version 5.11.3 via

pkg> compat CUDA 5

pkg> up

the MWE above passes for me as well. 🚀

[jonas-schulze]

It seems that the underlying issue persists (see #3236). How do I create a SubArray now? @vchuravy, you were helping me during SIAM PP on the same problem ... 🙈

Edit: never mind, I found a new way.