Scalar indexing with isapprox when we compare sparse matrices

[amontoison]

  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/x86_64/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:151
    [3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
      @ GPUArraysCore ~/.julia/x86_64/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:124
    [4] assertscalar(op::String)
      @ GPUArraysCore ~/.julia/x86_64/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:112
    [5] getindex
      @ ~/.julia/x86_64/packages/GPUArrays/w335n/src/host/indexing.jl:50 [inlined]
    [6] getindex(A::ROCSparseMatrixCSR{Float64, Int32}, i0::Int64, i1::Int64)
      @ AMDGPU.rocSPARSE ~/.julia/x86_64/packages/AMDGPU/TqRG0/src/sparse/array.jl:353
    [7] _getindex
      @ ./abstractarray.jl:1358 [inlined]
    [8] getindex
      @ ./abstractarray.jl:1312 [inlined]
    [9] iterate
      @ ./abstractarray.jl:1209 [inlined]
   [10] iterate
      @ ./abstractarray.jl:1207 [inlined]
   [11] _foldl_impl
      @ ./reduce.jl:56 [inlined]
   [12] foldl_impl
      @ ./reduce.jl:48 [inlined]
   [13] mapfoldl_impl
      @ ./reduce.jl:44 [inlined]
   [14] mapfoldl
      @ ./reduce.jl:175 [inlined]
   [15] _mapreduce
      @ ./reduce.jl:450 [inlined]
   [16] _mapreduce_dim
      @ ./reducedim.jl:337 [inlined]
   [17] mapreduce
      @ ./reducedim.jl:329 [inlined]
   [18] generic_normInf
      @ ~/Applications/julia/x86_64/julia-1.11.6/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:459 [inlined]
   [19] normInf
      @ ~/Applications/julia/x86_64/julia-1.11.6/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:533 [inlined]
   [20] generic_norm2(x::ROCSparseMatrixCSR{Float64, Int32})
      @ LinearAlgebra ~/Applications/julia/x86_64/julia-1.11.6/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:469
   [21] norm2
      @ ~/Applications/julia/x86_64/julia-1.11.6/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:535 [inlined]
   [22] norm(itr::ROCSparseMatrixCSR{Float64, Int32}, p::Int64)
      @ LinearAlgebra ~/Applications/julia/x86_64/julia-1.11.6/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:604
   [23] norm
      @ ~/Applications/julia/x86_64/julia-1.11.6/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:602 [inlined]
   [24] isapprox(x::ROCSparseMatrixCSR{Float64, Int32}, y::ROCSparseMatrixCSR{Float64, Int32}; atol::Int64, rtol::Float64, nans::Bool, norm::typeof(norm))
      @ LinearAlgebra ~/Applications/julia/x86_64/julia-1.11.6/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1806
   [25] isapprox(x::ROCSparseMatrixCSR{Float64, Int32}, y::ROCSparseMatrixCSR{Float64, Int32})
      @ LinearAlgebra ~/Applications/julia/x86_64/julia-1.11.6/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1802

[luraess]

Could you post a MWE?

[amontoison]

Yes @luraess

using AMDGPU
using AMDGPU.rocSPARSE

n = 100
d = 0.02
A = sprand(n,n,d)
B = sprand(n,n,d)
A ≈ B

A_gpu = ROCSparseMatrixCSR(A)
B_gpu = ROCSparseMatrixCSR(B)
A_gpu ≈ B_gpu

A_gpu = ROCSparseMatrixCSC(A)
B_gpu = ROCSparseMatrixCSC(B)
A_gpu ≈ B_gpu

A_gpu = ROCSparseMatrixCOO(A)
B_gpu = ROCSparseMatrixCOO(B)
A_gpu ≈ B_gpu

[luraess]

We may miss a method. I'll have a look. Else, feel free to draft a PR possibly getting inspiration on how other GPU backend handle this.

[gdalle]

Hi, gentle bump on this :) Apparently it is needed for AMD sparse matrices to work in SparseMatrixColorings.jl

[simeonschaub]

I believe the original issue is fixed. The only one that still errors for me is ROCSparseMatrixCOO due to isapprox requiring broadcasting, and I don't believe CUDA supports broadcasting for COO either (see e.g. https://github.com/JuliaGPU/CUDA.jl/blob/fbb90981cbde21d979087ad518a510f5b38f95b3/lib/cusparse/src/interfaces.jl#L363). @gdalle Is COO support required for SparseMatrixColorings?

[gdalle]

No, we can probably get started with CSC / CSR!

[luraess]

Closing this for now as CSC and CSR works. The remaining COO limitations may be brought up in another or follow-up issue.