Use GPU to generate known sparse matrix

Specific Domains GPU
1

Suppose we have a known matrix A:

\begin{equation} A=\left[\begin{array}{lll} {1} & {0} & {0} \\ {0} & {2} & {4} \\ {0} & {0} & {3} \end{array}\right] \end{equation}

At this point, I want to express A as a sparse matrix only by using SparseArrays in Julia. The code is as follows:

using SparseArrays

I = [1, 2, 2, 3]
J = [1, 2, 3, 3]
V = [1, 2, 4, 3]

A = sparse(I, J, V)

Result:

3×3 SparseMatrixCSC{Int64,Int64} with 4 stored entries:
  [1, 1]  =  1
  [2, 2]  =  2
  [2, 3]  =  4
  [3, 3]  =  3

Now I also have known I, J, and V. I want to generate sparse matrix in CSR or CSC format in GPU. The code is as follows:

using CuArrays

D_I = CuArray{Int64,1}(I)
D_J = CuArray{Int64,1}(J)
D_V = CuArray{Float32,1}(V)

CUSPARSE.sparse(D_I, D_J, D_V)

Result:

┌ Warning: Performing scalar operations on GPU arrays: This is very slow, consider disallowing these operations with `allowscalar(false)`
└ @ GPUArrays C:\Users\zenan\.julia\packages\GPUArrays\1wgPO\src\indexing.jl:16
3×3 SparseMatrixCSC{Float32,Int64} with 4 stored entries:
  [1, 1]  =  1.0
  [2, 2]  =  2.0
  [2, 3]  =  4.0
  [3, 3]  =  3.0

The way to generate sparse matrix in CSC format is very slow as shown in the result. Did I do something wrong? :sweat_smile:

2

This function does not exist in CUSPARSE, and falls back to Base where scalar iteration happens:

julia> @which CUSPARSE.sparse(D_I, D_J, D_V)
sparse(I, J, V::AbstractArray{T,1} where T) in SparseArrays at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.3/SparseArrays/src/sparsematrix.jl:845

julia> CuArrays.allowscalar(false)

julia> CUSPARSE.sparse(D_I, D_J, D_V)
ERROR: scalar getindex is disallowed
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] assertscalar(::String) at /home/tim/Julia/pkg/GPUArrays/src/indexing.jl:14
 [3] getindex at /home/tim/Julia/pkg/GPUArrays/src/indexing.jl:54 [inlined]
 [4] sparse!(::CuArray{Int64,1,Nothing}, ::CuArray{Int64,1,Nothing}, ::CuArray{Float32,1,Nothing}, ::Int64, ::Int64, ::typeof(+), ::Array{Int64,1}, ::Array{Int64,1}, ::Array{Int64,1}, ::Array{Float32,1}, ::Array{Int64,1}, ::Array{Int64,1}, ::Array{Float32,1}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.3/SparseArrays/src/sparsematrix.jl:730
 [5] sparse(::CuArray{Int64,1,Nothing}, ::CuArray{Int64,1,Nothing}, ::CuArray{Float32,1,Nothing}, ::Int64, ::Int64, ::Function) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.3/SparseArrays/src/sparsematrix.jl:660
 [6] sparse(::CuArray{Int64,1,Nothing}, ::CuArray{Int64,1,Nothing}, ::CuArray{Float32,1,Nothing}, ::Int64, ::Int64) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.3/SparseArrays/src/sparsematrix.jl:849
 [7] sparse(::CuArray{Int64,1,Nothing}, ::CuArray{Int64,1,Nothing}, ::CuArray{Float32,1,Nothing}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.3/SparseArrays/src/sparsematrix.jl:845
 [8] top-level scope at REPL[28]:1

You can however just call the CuSparseMatrixCSC constructor on your CPU SparseMatrixCSC, but that would reveal another nasty error:

julia> using CuArrays.CUSPARSE

julia> CuSparseMatrixCSC(A)
ERROR: MethodError: no method matching CuSparseMatrixCSC(::Array{Int64,1}, ::Array{Int64,1}, ::Array{Int64,1}, ::Tuple{Int64,Int64})

The problem here is that CUSPARSE does not support sparse matrices with integer types. Using Float32 (or any other supported type) fixes that:

julia> A = sparse(Float32.(I), Float32.(J), Float32.(V))
3×3 SparseMatrixCSC{Float32,Int64} with 4 stored entries:
  [1, 1]  =  1.0
  [2, 2]  =  2.0
  [2, 3]  =  4.0
  [3, 3]  =  3.0

julia> CuSparseMatrixCSC(A)
3×3 CuSparseMatrixCSC{Float32}:
Error showing value of type CuSparseMatrixCSC{Float32}

Everything works now, except displaying the sparse matrix. This is a known issue: Showing sparse arrays goes wrong · Issue #146 · JuliaGPU/CUDA.jl · GitHub (the sparse functionality has not been getting as much love as the rest of the package).

1 Like
3

Thank you very much. :grinning: I think I know what to do.