Sparse matrix handling

mutable struct ExtendableSparseMatrix{Tv, Ti<:Integer} <: SparseArrays.AbstractSparseMatrixCSC{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

ExtendableSparseMatrix(A)

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

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

Create ExtendableSparseMatrix from triplet (COO) data.

ExtendableSparseMatrix(D)

Create ExtendableSparseMatrix from Diagonal

SparseMatrixCSC(A)

Create SparseMatrixCSC from ExtendableSparseMatrix

*(d, ext)

*(ext, d)

+(ext, csc)

Add SparseMatrixCSC matrix and ExtendableSparseMatrix ext.

-(ext, csc)

Subtract SparseMatrixCSC matrix from ExtendableSparseMatrix ext.

-(csc, ext)

Subtract ExtendableSparseMatrix ext from SparseMatrixCSC.

\(symm_ext, b)

\ for Hermitian{ExtendableSparse}

\(symm_ext, b)

\ for Symmetric{ExtendableSparse}

A

\ for ExtendableSparse. 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 `` is called, which is realized via UMFPACK.

copy(S)

eltype(_)

Return element type.

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.

show(io, _, ext)

Show matrix

similar(m)

Create similar but emtpy extendableSparseMatrix

size(ext)

Size of 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)

Like updateindex! but without checking if v is zero.

cond(A)
cond(A, p)

flush! and calculate cond from cscmatrix

issymmetric(A)

flush! and check for symmetry of cscmatrix

ldiv!(r, ext, x)

flush! and ldiv with ext.cscmatrix

mul!(r, ext, x)

flush! and multiply with ext.cscmatrix

norm(A)
norm(A, p)

flush! and calculate norm from cscmatrix

opnorm(A)
opnorm(A, p)

flush! and calculate opnorm from cscmatrix

dropzeros!(ext)

findnz(ext)

flush! and return findnz(ext.cscmatrix).

getcolptr(ext)

flush! and return colptr of in ext.cscmatrix.

nnz(ext)

flush! and return number of nonzeros in ext.cscmatrix.

nonzeros(ext)

flush! and return nonzeros in ext.cscmatrix.

rowvals(ext)

flush! and return rowvals in ext.cscmatrix.

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.