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

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

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

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.