Update matrix block inplace from a sparse array

Hi,

Edit: sounds like it is related to BlockSparseMatrices.jl although the docs are quite minimal

I have a sparse matrix a which is made of blocks (a few).

using SparseArrays, LinearAlgebra
N = 100
# this is a completely made up example
a = kron(I(2), spdiagm(0 => ones(N), 1=>ones(N-1)))

I want to update the blocks as fast as possible. Let us say that the block (1,2) of size NxN has to be update by b. I can do

a[1:N,N+1:2N] .= b

Now, I want to speed this up. From

I,J,K = findnz(a)

I want to be able to extract the indices of I,J and K which correspond to a[1:N,N+1:2N]. I thought I could do the following but I get:

julia> findnz(view(a, 1:N,N+1:2N))
ERROR: MethodError: no method matching findnz(::SubArray{Float64,2,SparseMatrixCSC{Float64,Int64},Tuple{UnitRange{Int64},UnitRange{Int64}},false})
Closest candidates are:
  findnz(::SparseMatrixCSC{Tv,Ti}) where {Tv, Ti} at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/SparseArrays/src/sparsematrix.jl:1418
  findnz(::SparseVector{Tv,Ti}) where {Tv, Ti} at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/SparseArrays/src/sparsevector.jl:738
Stacktrace:
 [1] top-level scope at REPL[6]:100:

Do you have any idea on how to do this?

Thank you a lot for your suggestions…

Best regards,

I’ve been looking for this as well, and it looks like there’s some related work and interest from others. There is a related feature SparseArrays.sparse! (not exported) which allows you to get halfway there by specifying some pre-allocated arrays for the construction of sparse matrices.

I haven’t found a good solution for full in-place updates, however. As far as I can tell, it’s possible to do fully in-place sparse matrix updates by directly modifying rowptr,colptr,nzval, but more complicated when accounting for general sparse matrix functionality like combining repeated entries and sorted/unsorted representations.

1 Like

I have done this by looping at the non zero values. First I extract the indices of the sparse matrix in the OO format for each block, here. Then, I update the blocks here. For my use case, the problem is solved.

3 Likes