Broadcast of sparse vectors

[ChrisRackauckas]

Describe the bug

Some cases of appropriate sparse matrix broadcasts error on CUDA sparse matrices

To reproduce

The Minimal Working Example (MWE) for this bug:

using CUDA
A = zeros(Float32, N, N)
A[:,1] .= ones(N)
W_d = CuArray(W)
W_sparse = sparse(W)
A = CUDA.CUSPARSE.CuSparseMatrixCSC(sparse(A))
false .* A[1,:]

julia> false .* A[1,:]
ERROR: broadcast with sparse arrays is currently only implemented for CSR and CSC matrices
Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:35
 [2] copy(bc::Base.Broadcast.Broadcasted{CUDA.CUSPARSE.CuSparseVecStyle, Tuple{…}, typeof(*), Tuple{…}})
   @ CUDA.CUSPARSE C:\Users\accou\.julia\packages\CUDA\YIj5X\lib\cusparse\broadcast.jl:473
 [3] materialize(bc::Base.Broadcast.Broadcasted{CUDA.CUSPARSE.CuSparseVecStyle, Nothing, typeof(*), Tuple{…}})
   @ Base.Broadcast .\broadcast.jl:903
 [4] top-level scope
   @ REPL[14]:1
Some type information was truncated. Use `show(err)` to see complete types.

[maleadt]

Your MWE does not work (no definitions for N, W; W_sparse unused) and performs scalar iteration (slicing A).

Looking at the backtrace though, it looks like you're trying to:

julia> x = CuSparseVector(sprand(1, 1.0))
1-element CuSparseVector{Float64, Int32} with 1 stored entry:
  [1]  =  0.753707

julia> false .* x
ERROR: broadcast with sparse arrays is currently only implemented for CSR and CSC matrices

Is this what you need? In any case, this isn't a bug, but just a missing feature. I'm also not sure it's worth the effort, seeing how you're doing scalar iteration first? In that case, it's probably much better to just copy to the CPU first.

[ChrisRackauckas]

Yeah that seems to be the a simpler version of the same issue. It's mixing broadcast of sparse matrices or views/slices with scalars.

[maleadt]

It's mixing broadcast of sparse matrices or views/slices with scalars.

You'll have to be more specific then. Is it scalar .* CuSparseVector you need? Does it involve multiple sparse vectors? Or also involve sparse matrices? If so, which format, and how exactly? Since there isn't a way (I know of) to implement this generically, we need specialized kernels to handle exactly what you need.

[ChrisRackauckas]

Yes, specifically it's just trying to zero a sparse vector. false .* W[1,:] . fill!(W[1,:],false) also fails:

julia> fill!(W_d_sparse[1,:],false)
ERROR: CanonicalIndexError: setindex! not defined for CUDA.CUSPARSE.CuSparseVector{Float32, Int32}
Stacktrace:
 [1] error_if_canonical_setindex(::IndexCartesian, A::CUDA.CUSPARSE.CuSparseVector{Float32, Int32}, ::Int64)
   @ Base .\abstractarray.jl:1403
 [2] setindex!
   @ Base .\abstractarray.jl:1392 [inlined]
 [3] fill!(A::CUDA.CUSPARSE.CuSparseVector{Float32, Int32}, x::Bool)
   @ Base .\multidimensional.jl:1114
 [4] top-level scope
   @ REPL[2]:1