Factorizations & Preconditioners

This functionality is deprecated as of version 1.6 an will be removed with version 2.0. Use the integration with LinearSolve.jl.

Factorizations

In this package, preconditioners and LU factorizations are both seen as complete or approximate factorizations. Correspondingly we provide a common API for them.

ExtendableSparse.AbstractLUFactorization

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

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

source
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.

source

LU Factorizations

Handling of the LU factorizations is meant to support a workflow where sequences of problems are solved based on the same matrix, where one possibly wants to re-use existing symbolic factorization data.

The support comes in two flavors.

  • Using factorize! which can work as a drop-in replacement for lu!:
using ExtendableSparse, LinearAlgebra
A = fdrand(20, 20, 1; matrixtype = ExtendableSparseMatrix)
n = size(A, 1)
b = rand(n)
factorization = SparspakLU()
factorize!(factorization, A)
nm1 = norm(factorization \ b)

# mock update from Newton etc.
for i = 4:(n - 3)
    A[i, i + 3] -= 1.0e-4
end
factorize!(factorization, A)
nm2 = norm(factorization \ b)
nm1, nm2
(1746.0085061899467, 1779.2833284680926)
  • Using update!, where the matrix only needs to be given at construction time:
using ExtendableSparse, LinearAlgebra
A = fdrand(20, 20, 1; matrixtype = ExtendableSparseMatrix)
n = size(A, 1)
b = rand(n)
factorization = CholeskyFactorization(A)
nm1 = norm(factorization \ b)

# mock update from Newton etc.
for i = 4:(n - 3)
    A[i, i + 3] -= 1.0e-4
    A[i - 3, i] -= 1.0e-4
end
update!(factorization)
nm2 = norm(factorization \ b)
nm1, nm2
(1841.2789601350842, 1989.290578304452)
Base.:\ — Function
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.

source
\(symm_ext, b)

\ for Symmetric{ExtendableSparse}

source
\(symm_ext, b)

\ for Hermitian{ExtendableSparse}

source
 lufact
hs

Solve LU factorization problem.

source

Preconditioners

The API is similar to that for LU factorizations.

The support comes in two flavors.

using ExtendableSparse, LinearAlgebra
using IterativeSolvers, IncompleteLU
A = fdrand(20, 20, 1; matrixtype = ExtendableSparseMatrix)
n = size(A, 1)
b = rand(n)
preconditioner = ILUTPreconditioner(; droptol = 1.0e-2)
factorize!(preconditioner, A)

# mock update from Newton etc.
nm1 = norm(bicgstabl(A, b, 1; Pl = preconditioner))
for i = 4:(n - 3)
    A[i, i + 3] -= 1.0e-4
end
factorize!(preconditioner, A)
nm2 = norm(bicgstabl(A, b, 1; Pl = preconditioner))
nm1, nm2
(1962.2503660911302, 2002.358258851302)
using ExtendableSparse, LinearAlgebra
using IterativeSolvers
A = fdrand(20, 20, 1; matrixtype = ExtendableSparseMatrix)
n = size(A, 1)
b = rand(n)
preconditioner = ILU0Preconditioner(A)
nm1 = norm(cg(A, b; Pl = preconditioner))

# mock update from Newton etc.
for i = 4:(n - 3)
    A[i, i + 3] -= 1.0e-4
    A[i - 3, i] -= 1.0e-4
end
update!(preconditioner)
nm2 = norm(cg(A, b; Pl = preconditioner))
nm1, nm2
(2124.851317084639, 2216.6042573587065)
ExtendableSparse.BlockPreconBuilder — Type
 BlockPreconBuilder(;precs=UMFPACKPreconBuilder(),  
                     partitioning = A -> [1:size(A,1)]

Return callable object constructing a left block Jacobi preconditioner from partition of unknowns.

  • partitioning(A) shall return 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].

  • precs(A,p) shall return a left precondioner for a matrix block.

source
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.

source
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.

source

Experimental/deprecated

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)
source
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)
source

Iteration schemes

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

source