A version of sparse! that updates a previously existing SparseMatrixCSC?

Internals & Design
1

I will start with some context:

I am dealing with a problem where I have to diagonalize a matrix H for diferent values of a parameter k. I know that the matrix H is sparse, but I do not now a priori either the structure of the matrix or the number of non-zero entries. I only know the exact sparcity structure of the matrix once I try to build it. Initialy, I simply defined a function that for each value of k, allocated and built vector I, J & V with the non-zero elements and then called sparse(I, J, V). Obviously this is wasteful, since for each value of k I have to allocate a bunch of vectors and matrices.

So I was looking for a way to preallocate all the vectors and matrices before hand and then reuse them for each value of k. Preallocting the I``J``V is easy, and if I do not want to use the push! to incrementally build those vectors I can use PushVectors.jl.

Now for the preallocation of the sparse matrix itself. Looking through the source code, I found out that the function sparse actually calls a lower level function

sparse!(I, J, V, m, n, combine, klasttouch, csrrowptr, csrcolval, csrnzval,csccolptr, cscrowval, cscnzval)

which works directly with the underlying structure of a SparseMatrixCSC. This essentially solves my problem! The only thing, is that at the end, the function sparse! returns

SparseMatrixCSC(m, n, csccolptr, cscrowval, cscnzval)

instead os simply updating the vectors. I think it would be benifial to have a function updatesparsestructure!(I, J, V, m, n, combine, klasttouch, csrrowptr, csrcolval, csrnzval,csccolptr, cscrowval, cscnzval) that does the same as the current sparse! but then returns nothing. Then the current sparse! function would be simply:

function sparse!(I, J, V, m, n, combine, klasttouch, csrrowptr, csrcolval, csrnzval,csccolptr, cscrowval, cscnzval)

    updatesparsestructure!(I, J, V, m, n, combine, klasttouch, csrrowptr, csrcolval, csrnzval,csccolptr, cscrowval, cscnzval)

    return SparseMAtrixCSC(m, n, csccolptr, cscrowval, cscnzval)
end

A preexisting sparse matrix could be updated simply by defining a function:

function sparse!(spmat, klasttouch, csrrowptr, csrcolval, csrnzval, I, J, V, combine) 
    updatesparsestructure!(I, J, V, spmat.m, spmat.n, combine, klasttouch, csrrowptr, csrcolval, csrnzval, spmat.colptr, spmat.rowval, spmat.nzval)
end

For ease of use, the vectors klasttouch, csrrowptr, csrcolval and csrnzval could be wrapped in a struct.

This is what I did to solve my problem. And I think this would be a useful feature to have in Julia. Would anyone else find this useful? I don’t have much experience in collaborative devellopment, but could such a proposal be included in Julia? Or there are good reasons why this is not included?

Thank you for your time and sorry for the long post!

1 Like
2

Isn’t the return elided if you don’t use it?

3

I wonder if that is true.

I define a function sparse_return_spmatcsc! (which is just sparse!) and a function sparse_return_nothing! (which is sparse! but modified so that it does not call SparseMatrixCSC at the end and returns nothing). I then wrapped this functions into another one.

function update_csc_return_spmatcsc!(m, n, csccolptr, cscrowval, cscnzval, 
    klasttouch, csrrowptr, csrcolval, csrnzval, 
    I, J, V)
    
    sparse_return_spmatcsc!(m, n, csccolptr, cscrowval, cscnzval, 
        klasttouch, csrrowptr, csrcolval, csrnzval, 
        I, J, V, +)
    
    return csccolptr, cscrowval, cscnzval
    
end

and

function update_csc_return_nothing!(m, n, csccolptr, cscrowval, cscnzval, 
    klasttouch, csrrowptr, csrcolval, csrnzval, 
    I, J, V)
    
    sparse_return_nothing!(m, n, csccolptr, cscrowval, cscnzval, 
        klasttouch, csrrowptr, csrcolval, csrnzval, 
        I, J, V, +)
    
    return csccolptr, cscrowval, cscnzval
    
end

Then I tested:

m, n, nnzval = 1000, 800, 20

I = rand(1:m, nnzval)
J = rand(1:n, nnzval)
V = rand(Float64, nnzval)

klasttouch = zeros(Int, n)
csrrowptr = zeros(Int, m + 1)
csrcolval = zeros(Int, nnzval)
csrnzval = zeros(Float64, nnzval)

csccolptr = zeros(Int, n+1)
cscrowval = zeros(Int, nnzval)
cscnzval = zeros(Float64, nnzval)

@btime update_csc_return_spmatcsc!(m, n, csccolptr, cscrowval, cscnzval, 
    klasttouch, csrrowptr, csrcolval, csrnzval, 
    I, J, V)

m, n, nnzval = 1000, 800, 20

I = rand(1:m, nnzval)
J = rand(1:n, nnzval)
V = rand(Float64, nnzval)

klasttouch = zeros(Int, n)
csrrowptr = zeros(Int, m + 1)
csrcolval = zeros(Int, nnzval)
csrnzval = zeros(Float64, nnzval)

csccolptr = zeros(Int, n+1)
cscrowval = zeros(Int, nnzval)
cscnzval = zeros(Float64, nnzval)

@btime update_csc_return_nothing!(m, n, csccolptr, cscrowval, cscnzval, 
    klasttouch, csrrowptr, csrcolval, csrnzval, 
    I, J, V)

For update_csc_return_spmatcsc! I obtained:
3.871 μs (2 allocations: 80 bytes)
while for the update_csc_return_nothing! I got:
3.553 μs (1 allocation: 32 bytes)

So it seems that the call of SparseMatrixCSC is doing something.