# Sparse Matrix with CUDA.jl

**URL:** <https://discourse.julialang.org/t/sparse-matrix-with-cuda-jl/62627>\
**Category:** Performance\
**Tags:** question\
**Created:** [June 9, 2021, 2:24pm UTC](https://discourse.julialang.org/t/sparse-matrix-with-cuda-jl/62627 "2021-06-09T14:24:01Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![boutor2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/boutor2/32/27869_2.png) [@boutor2](https://discourse.julialang.org/u/boutor2)\
**Post date:** [June 9, 2021, 2:24pm UTC](https://discourse.julialang.org/t/sparse-matrix-with-cuda-jl/62627/1 "2021-06-09T14:24:01Z")

</div>

Hi everyone, I am looking for the most performant way to create a `CuArray` where coefficients are `0` everywhere but `1` at specified indices.  
An easy way to do that with regular arrays would be

```nohighlight
a = randn(1000,1000)
imin = argmin(a,dims=1) # coefficients where we want b[i] =1
b = zeros(size(a))
b[imin] .= 1

```

But it gets trickier with `CuArray`s

```julia
using CUDA
CUDA.allowscalar(false)
a = CUDA.randn(1000,1000)
imin = argmin(a,dims=1)

```

One cannot broadcast in a similar way as above, since we have not allowed scalar because of performance issues. One could do

```nohighlight
imin = argmin(a,dims=1) |> Array
b = zeros(size(a))
b[imin] .= 1
b = b |> CuArray

```

but this involves some back and forth between the gpu and the cpu, that is not nice.  
Any idea of a trick to get around this problem? Cheers!

---

<div class="post-metadata">

**Author:** ![yixingfu](https://avatars.discourse-cdn.com/v4/letter/y/e95f7d/32.png) [@yixingfu](https://discourse.julialang.org/u/yixingfu)\
**Post date:** [June 9, 2021, 3:19pm UTC](https://discourse.julialang.org/t/sparse-matrix-with-cuda-jl/62627/2 "2021-06-09T15:19:44Z")

</div>

This may not be optimal, but seem to be faster on average:

```julia
b = CUDA.zeros(size(a))
b .= a.==minimum(a;dims=1);

```

I’d also be interested to know what is the best way to do this.

---

<div class="post-metadata">

**Author:** ![eliassno](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eliassno/32/18917_2.png) [@eliassno](https://discourse.julialang.org/u/eliassno)\
**Post date:** [June 9, 2021, 5:13pm UTC](https://discourse.julialang.org/t/sparse-matrix-with-cuda-jl/62627/3 "2021-06-09T17:13:25Z")

</div>

You may want to consider constructing the array in a CSC format.  
Note that the timings for `f_csc` are only relevant for the provided MWE, but you would have to rethink how you construct the sparsity pattern in practice.

```julia
using CUDA
using SparseArrays
CUDA.allowscalar(false)

function f_original(N)
    a = randn(N, N)
    imin = argmin(a,dims=1) |> Array
    b = zeros(size(a))
    b[imin] .= 1
    b = b |> CuArray
    return b
end

function f_csc(N)
    dcolptr = cu(collect(1:N+1))
    drowval = cu(rand(1:N, N))
    dnzval = CUDA.ones(Int64, N)
    dsa = CUSPARSE.CuSparseMatrixCSC(dcolptr, drowval, dnzval, (N, N))
    return dsa
end

for N_pow in 0:12
    N = 2^N_pow
    println("\nN = $N")
    @time CUDA.@sync f_original(N)
    @time CUDA.@sync f_csc(N)
end

```

> **@time results**
>
> ```julia
> N = 1
> 0.005567 seconds (24 allocations: 848 bytes)
> 0.000562 seconds (94 allocations: 3.422 KiB)
> 
> N = 2
> 0.000127 seconds (27 allocations: 1.281 KiB)
> 0.000565 seconds (94 allocations: 3.453 KiB)
> 
> N = 4
> 0.000173 seconds (27 allocations: 1.547 KiB)
> 0.000665 seconds (94 allocations: 3.547 KiB)
> 
> N = 8
> 0.000155 seconds (27 allocations: 2.516 KiB)
> 0.000512 seconds (94 allocations: 3.734 KiB)
> 
> N = 16
> 0.000131 seconds (27 allocations: 5.859 KiB)
> 0.000475 seconds (94 allocations: 4.109 KiB)
> 
> N = 32
> 0.000124 seconds (24 allocations: 18.156 KiB)
> 0.000488 seconds (94 allocations: 4.891 KiB)
> 
> N = 64
> 0.000157 seconds (26 allocations: 67.406 KiB)
> 0.000507 seconds (94 allocations: 6.500 KiB)
> 
> N = 128
> 0.000238 seconds (35 allocations: 262.094 KiB)
> 0.000495 seconds (94 allocations: 9.688 KiB)
> 
> N = 256
> 0.000495 seconds (26 allocations: 1.011 MiB)
> 0.000508 seconds (94 allocations: 15.797 KiB)
> 
> N = 512
> 0.001965 seconds (26 allocations: 4.020 MiB)
> 0.000523 seconds (94 allocations: 27.828 KiB)
> 
> N = 1024
> 0.013047 seconds (203 allocations: 16.043 MiB, 41.24% gc time)
> 0.001215 seconds (108 allocations: 55.883 KiB)
> 
> N = 2048
> 0.031436 seconds (44 allocations: 64.079 MiB, 7.08% gc time)
> 0.000608 seconds (101 allocations: 100.219 KiB)
> 
> N = 4096
> 0.115517 seconds (48 allocations: 256.158 MiB, 2.73% gc time)
> 0.000693 seconds (108 allocations: 196.422 KiB)
> 
> ```

---

<div class="post-metadata">

**Author:** ![boutor2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/boutor2/32/27869_2.png) [@boutor2](https://discourse.julialang.org/u/boutor2)\
**Post date:** [June 10, 2021, 4:14pm UTC](https://discourse.julialang.org/t/sparse-matrix-with-cuda-jl/62627/4 "2021-06-10T16:14:32Z")

</div>

Thank you very much for your comprehensive answer. However, it does not seem trivial to convert `CartesianIndex` in a CSC format. But may be I am wrong?

Here is a proposal to do so, that uses the `sparse` method from `SparseArrays`

> [@Cartesian indices, linear indices, and sparse matrices](https://discourse.julialang.org/t/cartesian-indices-linear-indices-and-sparse-matrices/13781/2):
>
> Thanks for writing this up as a comprehensive question! It’d probably be nice to add a sparse method that takes arrays of CartesianIndexes. Here’s a simple implementation that could serve you in the mean-time: function SparseArrays.sparse(IJ::Vector{\<:CartesianIndex}, v, m, n) IJ′ = reinterpret(Int, reshape(IJ, 1, :)) return sparse(view(IJ′, 1, :), view(IJ′, 2, :), v, m, n) end In the above I used reinterpret and lazy views to prevent needlessly allocating new vectors that contai…

but as you mentioned in an other post we miss the `sparse` method for `CuSparseMatricCSC`

> [@CuSparseMatrixCSC constructor](https://discourse.julialang.org/t/cusparsematrixcsc-constructor/62480/2):
>
> This only works for [sparse from SparseArrays](https://github.com/JuliaLang/julia/blob/master/stdlib/SparseArrays/src/sparsematrix.jl#L732-L799): julia\> using SparseArrays julia\> sparse(Array(I), Array(J), Array(V), 3, 3) 3×3 SparseMatrixCSC{Float64,Int64} with 2 stored entries: [1, 1] = 0.2 [1, 2] = 0.3 You need to convert your coordinate lists (I and J) to CSC format (colPtr and rowVal). julia\> colPtr = CuArray([1, 2, 3, 3]); julia\> rowVal = CuArray([1, 1]); julia\> CUSPARSE.CuSparseMatrixCSC(colPtr, rowVal, V, (3,3)) 3×3 CuSparseMatrixCSC{Float64} with 2 stored entries: [1, 1…