1

I read this post and I have a similar problem. I want to in-place update sparse matrices multiple times, and it seems to me that SparseArrays.sparse! is my best bet. However, I don't quite understand the documentation, which says I may call

sparse!(I, J, V, csrrowptr, csrcolval, csrnzval, I, J, V)

My questions are:

  1. What is the first set of I, J, V?
  2. What is the second set of I, J, V?
  3. Suppose I have created a matrix A = sparse(I,J,V) and I want to in-place update A. How can I get csrrowptr, csrcolval, csrnzval? A by default has a CSC format, not CSR.
  4. Can I get a working example of sparse!?

Thank you.

user1691278
  • 1,751
  • 1
  • 15
  • 20
  • The docstring, together with the one from `sparse`, is quite exhaustive. What are information is missing exactly from there? – phipsgabler Dec 01 '22 at 08:00

1 Answers1

1

Having looked at the documentation, it definitely feels (and is) a package internal function of SparseArrays. In any case, per your questions (not in order):

Q4 (the easiest) do @edit sparse([1,2,3],[1,2,3],[5,6,7]) at the REPL, and search for sparse! in the edited file, this will show how sparse uses sparse! internally.

Q1 The first I, J, V are the inputs, giving the usual columns indices, row indices, values. Note, these are not the internal CSC storage vectors of the resulting matrix.

Q2 The second I, J, V are the arrays used in the resulting sparse matrix as the internal CSC storage vectors. Since the input is not necessary after sparse! call, instead of deallocating input and allocating internal CSC vectors, the input memory vectors can be reused (other vectors than I, J, K can be passed as output, especially if the input vectors have too little memory to hold CSC representation).

Q3 Getting the CSR-ish representation created inside sparse! is easy - they are in the vectors you pass as csrrowptr, csrcolval, csrnzval. The vectors need to have enough memory, and are populated by sparse!.

The whole point of sparse! is to make no allocations/deallocations. The memory management is done by a 'parent' interface (sparse in this case). This allows direct non-allocating processing when the memory is already allocated or recycled from another variable.

It took a bit of fiddling, but here is a demo of sparse! usage:

julia> I = sizehint!([1,2,3,1], 5)
4-element Vector{Int64}:
 1
 2
 3
 1

julia> J = sizehint!([3, 2, 1, 2], 5)
4-element Vector{Int64}:
 3
 2
 1
 2

julia> V = sizehint!([6,7,8,9],5)
4-element Vector{Int64}:
 6
 7
 8
 9

Now setting up working memory for sparse!:

julia> my_klasttouch = Vector{Int}(undef, 10);

julia> my_csrrowptr = Vector{Int}(undef, 10);

julia> my_csrcolval = Vector{Int}(undef, 10);

julia> my_csrnzval = Vector{Int}(undef, 10);

Finally:

julia> M = SparseArrays.sparse!(
  I, J, V, 3, 3, +, 
  my_klasttouch, my_csrrowptr, my_csrcolval, my_csrnzval, 
  I, J, V
)
3×3 SparseMatrixCSC{Int64, Int64} with 4 stored entries:
 ⋅  9  6
 ⋅  7  ⋅
 8  ⋅  ⋅

To see how the input memory becomes part of the returned object:

julia> M.colptr   # this is from returned object
4-element Vector{Int64}:
 1
 2
 4
 5

julia> I
4-element Vector{Int64}:
 1
 2
 4
 5

Note that I and M.colpotr use the same memory.

Given the above, I would suggest, before using this function, to go over the implementation code and understand the algorithm carefully.

Dan Getz
  • 17,002
  • 2
  • 23
  • 41