Sparse matrix class with efficient successive insertion of entries and entry update.
Without an intermediate data structure, efficient successive insertion/update of possibly duplicate entries in random order into a standard compressed column storage structure appears to be not possible. The package introduces
ExtendableSparseMatrix, a delegating wrapper containing a Julia standard
SparseMatrixCSC struct for performing linear algebra operations and a
SparseMatrixLNK struct realising a linked list based (but realised in vectors) format collecting new entries.
Any linear algebra method on
ExtendableSparseMatrix starts with a
flush! method which adds the LNK entries and the existing CSC entries into a new CSC struct and resets the LNK struct.
ExtendableSparseMatrix is aimed to work as a drop-in replacement to
SparseMatrixCSC in finite element and finite volume codes especally in those cases where the sparsity structure is hard to detect a priori and where working with an intermediadte COO representation appears to be not convenient.
This package assumes that a $m \times n$ matrix is sparse if each row and each column have less than $C$ entries with $C << n$ and $C <<m$ . Adding a full matrix row will be a performance hit.
Working with ForwardDiff
In particular, it cooperates with ForwardDiff.jl when it comes to the assembly of a sparse jacobian. For a function 'f!(y,x)' returning it's result in a vector
y, one can use e.g.
x=... y=zeros(n) dresult=DiffResults.DiffResult(zeros(n),ExtendableSparseMatrix(n,n)) x=ForwardDiff.jacobian!(dresult,f!,y,x) jac=DiffResults.jacobian(dresult) h=jac\x
However, without a priori information on sparsity, ForwardDiff calls element insertion for the full range of n^2 indices, leading to a O(n^2) scaling behavior due to the nevertheless necessary search operations, see this discourse thread.
In addition, the package provides a method
updateindex!(A,op,v,i,j) for both
SparseMatrixCSC and for
ExtendableSparse which allows to update a matrix element with one index search instead of two. It allows to replace e.g.
updateindex!(A,+,v,i,j). The former operation is lowered to
%1 = Base.getindex(A, 1, 2) %2 = %1 + 3 Base.setindex!(A, %2, 1, 2)
triggering two index searches, one for
getindex! and another one for
See Julia issue #15630 for a discussion on this.