Sparse matrix handling

ExtendableSparseMatrix

Aliased to ExtendableSparseMatrixCSC

mutable struct ExtendableSparseMatrixCSC{Tv, Ti<:Integer} <: ExtendableSparse.AbstractExtendableSparseMatrixCSC{Tv, Ti<:Integer}

Extendable sparse matrix. A nonzero entry of this matrix is contained either in cscmatrix, or in lnkmatrix, never in both.

• cscmatrix::SparseArrays.SparseMatrixCSC: Final matrix data

• lnkmatrix::Union{Nothing, SparseMatrixLNK{Tv, Ti}} where {Tv, Ti<:Integer}: Linked list structure holding data of extension

• phash::UInt64: Pattern hash

ExtendableSparseMatrixCSC(A)

Create ExtendableSparseMatrixCSC from AbstractMatrix, dropping all zero entries. This is the equivalent to sparse(A).

ExtendableSparseMatrixCSC(I,J,V)
ExtendableSparseMatrixCSC(I,J,V,m,n)
ExtendableSparseMatrixCSC(I,J,V,combine::Function)
ExtendableSparseMatrixCSC(I,J,V,m,n,combine::Function)

Create ExtendableSparseMatrixCSC from triplet (COO) data.

ExtendableSparseMatrixCSC(D)

Create ExtendableSparseMatrixCSC from Diagonal

ExtendableSparseMatrixCSC(csc)

Create ExtendableSparseMatrixCSC from SparseMatrixCSC

copy(S)

getindex(ext, i, j)

Find index in CSC matrix and return value, if it exists. Otherwise, return value from extension.

setindex!(ext, v, i, j)

Find index in CSC matrix and set value if it exists. Otherwise, set index in extension if v is nonzero.

similar(m)

Create similar but emtpy extendableSparseMatrix

flush!(ext)

If there are new entries in extension, create new CSC matrix by adding the cscmatrix and linked list matrix and reset the linked list based extension.

pointblock(matrix,blocksize)

Create a pointblock matrix.

rawupdateindex!(ext, op, v, i, j)
rawupdateindex!(ext, op, v, i, j, part)

Like updateindex! but without checking if v is zero.

reset!(A)

Reset ExtenableSparseMatrix into state similar to that after creation.

lu(matrix)

Create LU factorization. It calls the LU factorization form Sparspak.jl, unless GPL components are allowed in the Julia sysimage and the floating point type of the matrix is Float64 or Complex64. In that case, Julias standard lu is called, which is realized via UMFPACK.

lu!(factorization, matrix)

Update LU factorization, possibly reusing information from the current state. This method is aware of pattern changes.

If nothing is passed as first parameter, factorize! is called.

ldiv!(r, ext, x)

flush! and ldiv with ext.cscmatrix

ldiv!(u,factorization,v)
ldiv!(factorization,v)

Solve factorization.

mark_dirichlet(A; penalty=1.0e20)

Return boolean vector marking Dirichlet nodes, known by A[i,i]>=penalty

eliminate_dirichlet!(A,dirichlet_marker)

Eliminate dirichlet nodes in matrix by setting

    A[:,i]=0; A[i,:]=0; A[i,i]=1

for a node i marked as Dirichlet.

Returns A.

eliminate_dirichlet(A,dirichlet_marker)

Create a copy B of A sharing the sparsity pattern. Eliminate dirichlet nodes in B by setting

    B[:,i]=0; B[i,:]=0; B[i,i]=1

for a node i marked as Dirichlet.

Returns B.

fdrand!(A, nx; ...)
fdrand!(A, nx, ny; ...)
fdrand!(A, nx, ny, nz; update, rand)

After setting all nonzero entries to zero, fill resp. update matrix with finite difference discretization data on a unit hypercube. See fdrand for documentation of the parameters.

It is required that size(A)==(N,N) where N=nx*ny*nz

fdrand(, nx; ...)
fdrand(, nx, ny; ...)
fdrand(, nx, ny, nz; matrixtype, update, rand, symmetric)
fdrand(nx; ...)

Create matrix for a mock finite difference operator for a diffusion problem with random coefficients on a unit hypercube $\Omega\subset\mathbb R^d$. with $d=1$ if nx>1 && ny==1 && nz==1, $d=2$ if nx>1 && ny>1 && nz==1 and $d=3$ if nx>1 && ny>1 && nz>1 . In the symmetric case it corresponds to

\[ \begin{align*} -\nabla a \nabla u &= f&& \text{in}\; \Omega \\ a\nabla u\cdot \vec n + bu &=g && \text{on}\; \partial\Omega \end{align*}\]

The matrix is irreducibly diagonally dominant, has positive main diagonal entries and nonpositive off-diagonal entries, hence it has the M-Property. Therefore, its inverse will be a dense matrix with positive entries, and the spectral radius of the Jacobi iteration matrix $ho(I-D(A)^{-1}A)<1$ .

Moreover, in the symmetric case, it is positive definite.

Parameters+ default values:

The sparsity structure is fixed to an orthogonal grid, resulting in a 3, 5 or 7 diagonal matrix depending on dimension. The entries are random unless e.g. rand=()->1 is passed as random number generator. Tested for Matrix, SparseMatrixCSC, ExtendableSparseMatrix, Tridiagonal, SparseMatrixLNK and :COO

sprand!(A, xnnz)

Fill empty sparse matrix A with random nonzero elements from interval [1,2] using incremental assembly.

sprand_sdd!(A; nnzrow)

Fill sparse matrix with random entries such that it becomes strictly diagonally dominant and thus invertible and has a fixed number nnzrow (default: 4) of nonzeros in its rows. The matrix bandwidth is bounded by sqrt(n) in order to resemble a typical matrix of a 2D piecewise linear FEM discretization.