Repository files navigation

Sparse matrix class with efficient successive insertion of entries and entry update, supporting general number types.

The package allows for efficient assembly of a sparse matrix withoutthe need to handle intermediate arrays:

using ExtendableSparseA=ExtendableSparseMatrix(10,10)A[1,1]=1for i = 1:9   A[i + 1, i] += -1   A[i, i + 1] += -1   A[i, i] += 1   A[i + 1, i + 1] += 1endb=ones(10)x=A\b

While one could replace hereExtendableSparseMatrix(10,10) byspzeros(10,10), the later is inefficient for large problems. Without this package, the general advise is toconstruct a sparse matrix via the COO format.

Instead of\, the methods fromLinearSolve.jl can be used:

using LinearSolvep=LinearProblem(A,b)LinearSolve.solve(p)

With the help ofSparspak.jl, these examples work for general number types.

sparse(A) creates a standardSparseMatrixCSC from a filledExtendableSparseMatrix which can be used e.g. to create preconditioners. So one can instead perform e.g.

LinearSolve.solve(p, KrylovJL_CG(); Pl = ILUZero.ilu0(sparse(A)))

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 introducesExtendableSparseMatrix, a delegating wrapper containing a Julia standardSparseMatrixCSC struct for performing linear algebra operations and aSparseMatrixLNK struct realising a linked list based (but realised in vectors) format collecting new entries.

The later is modeled after the linked list sparse matrix format described in thewhitepaper by Y. Saad. See also exercise P.3-16 in hisbook.

Any linear algebra method onExtendableSparseMatrix starts with aflush! 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 toSparseMatrixCSC in finite element and finite volume codes especially 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.

The package provides a\ method forExtendableSparseMatrix which dispatches to Julia's standard\ method forSparseMatrixCSC where possible.It relies onSparspak.jl, P.Krysl's Julia MIT licensed re-implementation of Sparspak by George & Liu fornumber types not supported by Julia's standard implementation. Notably, this includesForwardDiff.Dual numbers, thus supporting for automatic differentiation. When used with a non-GPL version of the system image,\ is dispatched to Sparsepak.jl in all cases.

This package assumes that a$m \times n$ matrix is sparse ifeach row andeach column have less than$C$ entries with$C << n$ and$C << m$ . Adding a full matrix row will be a performance hit.

In particular, it cooperates withForwardDiff.jl when it comes to the assembly of a sparse jacobian. For a function 'f!(y,x)' returning it's result in a vectory, 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 thisdiscourse thread.

In addition, the package provides a methodupdateindex!(A,op,v,i,j) for bothSparseMatrixCSC and forExtendableSparse which allows to update a matrix element with one index search instead of two. It allows to replace e.g.A[i,j]+=v byupdateindex!(A,+,v,i,j). The former operation is lowered to

%1 = Base.getindex(A, 1, 2)%2 = %1 + 3Base.setindex!(A, %2, 1, 2)

triggering two index searches, one forgetindex! and another one forsetindex!.

SeeJulia issue #15630 for a discussion on this.

The package provides a common API for factorizations and preconditioners supportingseries of solutions of similar problem as they occur during nonlinear and transient solves.For details, see thecorresponding documentation.

With the advent of LinearSolve.jl, this functionality probably will be reduced to some core cases.

The package directly provides interfaces to other sparse matrix solvers and preconditioners. Dependencies on thesepackages are handled viaRequires.jl.Currently, support includes:

• Pardiso.jl (both"project Pardiso"andMKL Pardiso)
• IncompleteLU.jl
• AlgebraicMultigrid.jl (Ruge-Stüben AMG)

For a similar approach, seePreconditioners.jl

You may also evaluate alternatives to this package:

• DynamicSparseArrays.jl