ExtendableSparse.AbstractFactorization — Type
abstract type AbstractFactorization

Abstract type for a factorization with ExtandableSparseMatrix.

This type is meant to be a "type flexible" (with respect to the matrix element type) and lazily construcdet (can be constructed without knowing the matrix, and updated later) LU factorization or preconditioner. It wraps different concrete, type fixed factorizations which shall provide the usual ldiv! methods.

Any such preconditioner/factorization MyFact should have the following fields

  A::ExtendableSparseMatrix
  factorization
  phash::UInt64

and provide methods

  MyFact(;kwargs...) 
  update!(precon::MyFact)

The idea is that, depending if the matrix pattern has changed, different steps are needed to update the preconditioner.

ExtendableSparse.AbstractSparseMatrixExtension — Type
abstract type AbstractSparseMatrixExtension{Tv, Ti} <: SparseArrays.AbstractSparseArray{Tv, Ti, 2}

Abstract type for sparse matrix extension.

Subtypes T_ext must implement:

  • Constructor T_ext(m,n)
  • SparseArrays.nnz(ext::T_ext)
  • Base.size(ext::T_ext)
  • Base.sum(extmatrices::Vector{T_ext}, csx): add csr or csc matrix and extension matrices (one per partition) and return csx matrix
  • Base.+(ext::T_ext, csx) (optional) - Add extension matrix and csc/csr matrix, return csx matrix
  • rawupdateindex!(ext::Text, op, v, i, j, tid) where {Tv, Ti}: Set ext[i,j]op=v, possibly insert new entry into matrix. tid is a

task or partition id

ExtendableSparse.BlockPreconditioner — Type
 BlockPreconditioner(;partitioning, factorization=LUFactorization)

Create a block preconditioner from partition of unknowns given by partitioning, a vector of AbstractVectors describing the indices of the partitions of the matrix. For a matrix of size n x n, e.g. partitioning could be [ 1:n÷2, (n÷2+1):n] or [ 1:2:n, 2:2:n]. Factorization is a callable (Function or struct) which allows to create a factorization (with ldiv! methods) from a submatrix of A.

ExtendableSparse.ExtendableSparseMatrixCSC — Type
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
ExtendableSparse.ExtendableSparseMatrixCSC — Method
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.

ExtendableSparse.LUFactorization — Type
LUFactorization()
LUFactorization(matrix)

Default LU Factorization. Maps to Sparspak.jl for non-GPL builds, otherwise to UMFPACK.

ExtendableSparse.SparseMatrixDILNKC — Type
mutable struct SparseMatrixDILNKC{Tv, Ti<:Integer} <: ExtendableSparse.AbstractSparseMatrixExtension{Tv, Ti<:Integer}
Modification of SparseMatrixLNK where the pointer to first index of

column j is stored in a dictionary.

ExtendableSparse.SparseMatrixLNK — Type
mutable struct SparseMatrixLNK{Tv, Ti<:Integer} <: ExtendableSparse.AbstractSparseMatrixExtension{Tv, Ti<:Integer}

Struct to hold sparse matrix in the linked list format.

Modeled after the linked list sparse matrix format described in the whitepaper and the SPARSEKIT2 source code by Y. Saad. He writes "This is one of the oldest data structures used for sparse matrix computations."

The relevant source formats.f is also available in the debian/science gitlab.

The advantage of the linked list structure is the fact that upon insertion of a new entry, the arrays describing the structure can grow at their respective ends and can be conveniently updated via push!. No copying of existing data is necessary.

  • m::Integer: Number of rows
  • n::Integer: Number of columns
  • nnz::Integer: Number of nonzeros
  • nentries::Integer: Length of arrays

  • colptr::Vector{Ti} where Ti<:Integer: Linked list of column entries. Initial length is n, it grows with each new entry.

    colptr[index] contains the next index in the list or zero, in the later case terminating the list which starts at index 1<=j<=n for each column j.

  • rowval::Vector{Ti} where Ti<:Integer: Row numbers. For each index it contains the zero (initial state) or the row numbers corresponding to the column entry list in colptr.

    Initial length is n, it grows with each new entry.

  • nzval::Vector: Nonzero entry values correspondin to each pair (colptr[index],rowval[index])

    Initial length is n, it grows with each new entry.

Base.:\ — Method
 lufact
hs

Solve LU factorization problem.

Base.:\ — Method
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.

Base.:\ — Method
\(symm_ext, b)

\ for Hermitian{ExtendableSparse}

Base.:\ — Method
\(symm_ext, b)

\ for Symmetric{ExtendableSparse}

Base.getindex — Method
getindex(lnk, i, j)

Return value stored for entry or zero if not found

Base.getindex — Method
getindex(ext, i, j)

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

Base.getindex — Method
getindex(lnk, i, j)

Return value stored for entry or zero if not found

Base.setindex! — Method
setindex!(lnk, v, i, j)

Update value of existing entry, otherwise extend matrix if v is nonzero.

Base.setindex! — Method
setindex!(lnk, v, i, j)

Update value of existing entry, otherwise extend matrix if v is nonzero.

Base.setindex! — Method
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.

Base.similar — Method
similar(m)

Create similar but emtpy extendableSparseMatrix

Base.size — Method
size(lnk)

Return tuple containing size of the matrix.

Base.size — Method
size(lnk)

Return tuple containing size of the matrix.

ExtendableSparse.MKLPardisoLU — Function
MKLPardisoLU(;iparm::Vector, mtype::Int)

MKLPardisoLU(matrix; iparm, mtype)

LU factorization based on pardiso. For using this, you need to issue using Pardiso. This version uses the early 2000's fork in Intel's MKL library.

The optional keyword arguments mtype and iparm are Pardiso internal parameters.

For setting them you can also access the PardisoSolver e.g. like

using Pardiso
plu=MKLPardisoLU()
Pardiso.set_iparm!(plu.ps,5,13.0)
ExtendableSparse.PardisoLU — Function
PardisoLU(;iparm::Vector, 
           dparm::Vector, 
           mtype::Int)

PardisoLU(matrix; iparm,dparm,mtype)

LU factorization based on pardiso. For using this, you need to issue using Pardiso and have the pardiso library from pardiso-project.org installed.

The optional keyword arguments mtype, iparm and dparm are Pardiso internal parameters.

Forsetting them, one can also access the PardisoSolver e.g. like

using Pardiso
plu=PardisoLU()
Pardiso.set_iparm!(plu.ps,5,13.0)
ExtendableSparse.allow_views — Method
allow_views(::preconditioner_type)

Factorizations on matrix partitions within a block preconditioner may or may not work with array views. E.g. the umfpack factorization cannot work with views, while ILUZeroPreconditioner can. Implementing a method for allow_views returning false resp. true allows to dispatch to the proper case.

ExtendableSparse.eliminate_dirichlet! — Method
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.

ExtendableSparse.eliminate_dirichlet — Method
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.

ExtendableSparse.factorize! — Method
factorize!(factorization, matrix)

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

ExtendableSparse.fdrand! — Method
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

ExtendableSparse.fdrand — Method
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:

Parameter + default valeDescription
nxNumber of unknowns in x direction
nyNumber of unknowns in y direction
nzNumber of unknowns in z direction
matrixtype = SparseMatrixCSCMatrix type
update = (A,v,i,j)-> A[i,j]+=vElement update function
rand =()-> rand()Random number generator
symmetric=trueWhether to create symmetric matrix or not

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

ExtendableSparse.findindex — Method
findindex(csc, i, j)

Return index corresponding to entry [i,j] in the array of nonzeros, if the entry exists, otherwise, return 0.

ExtendableSparse.flush! — Method
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.

ExtendableSparse.flush! — Method
flush!(csc)

Trival flush! method for allowing to run the code with both ExtendableSparseMatrix and SparseMatrixCSC.

ExtendableSparse.lu — Function
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.

ExtendableSparse.pattern_equal — Method
pattern_equal(a::SparseMatrixCSC,b::SparseMatrixCSC)

Check if sparsity patterns of two SparseMatrixCSC objects are equal. This is generally faster than comparing hashes.

ExtendableSparse.rawupdateindex! — Method
rawupdateindex!(lnk, op, v, i, j)

Update element of the matrix with operation op. It assumes that op(0,0)==0. If v is zero a new entry is created nevertheless.

ExtendableSparse.rawupdateindex! — Method
rawupdateindex!(lnk, op, v, i, j)

Update element of the matrix with operation op. It assumes that op(0,0)==0. If v is zero a new entry is created nevertheless.

ExtendableSparse.simple! — Method
simple!(u,A,b;
                 abstol::Real = zero(real(eltype(b))),
                 reltol::Real = sqrt(eps(real(eltype(b)))),
                 log=false,
                 maxiter=100,
                 P=nothing
                 ) -> solution, [history]

Simple iteration scheme $u_{i+1}= u_i - P^{-1} (A u_i -b)$ with similar API as the methods in IterativeSolvers.jl.

ExtendableSparse.sprand! — Method
sprand!(A, xnnz)

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

ExtendableSparse.sprand_sdd! — Method
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.

ExtendableSparse.updateindex! — Method
updateindex!(lnk, op, v, i, j)

Update element of the matrix with operation op. It assumes that op(0,0)==0. If v is zero, no new entry is created.

ExtendableSparse.updateindex! — Method
updateindex!(csc, op, v, i, j)

Update element of the matrix with operation op whithout introducing new nonzero elements.

This can replace the following code and save one index search during acces:

using ExtendableSparse # hide
using SparseArrays # hide
A=spzeros(3,3)
A[1,2]+=0.1
A
using ExtendableSparse # hide
using SparseArrays # hide
A=spzeros(3,3)
updateindex!(A,+,0.1,1,2)
A
ExtendableSparse.updateindex! — Method
updateindex!(lnk, op, v, i, j)

Update element of the matrix with operation op. It assumes that op(0,0)==0. If v is zero, no new entry is created.

LinearAlgebra.ldiv! — Method
ldiv!(u,factorization,v)
ldiv!(factorization,v)

Solve factorization.

LinearAlgebra.lu! — Method
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.

ExtendableSparse.@makefrommatrix — Macro

" @makefrommatrix(fact)

For an AbstractFactorization MyFact, provide methods

    MyFact(A::ExtendableSparseMatrix; kwargs...)
    MyFact(A::SparseMatrixCSC; kwargs...)