Linear System solution

The \ operator

The packages overloads \ for the ExtendableSparseMatrix type. The following example uses fdrand to create a test matrix and solve the corresponding linear system of equations.

using ExtendableSparse
A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix)
x = ones(1000)
b = A * x
y = A \ b
sum(y)
1000.0000000000001

This works as well for number types besides Float64 and related, in this case, by default a LU factorization based on Sparspak ist used.

using ExtendableSparse
using MultiFloats
A = fdrand(Float64x2, 10, 10, 10; matrixtype = ExtendableSparseMatrix)
x = ones(Float64x2,1000)
b = A * x
y = A \ b
sum(y)
999.999999999999999999999999999862690973379727207

Solving with LinearSolve.jl

Starting with version 0.9.6, ExtendableSparse is compatible with LinearSolve.jl. Since version 0.9.7, this is facilitated via the AbstractSparseMatrixCSC interface.

The same problem can be solved via LinearSolve.jl:

using ExtendableSparse
using LinearSolve
A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix)
x = ones(1000)
b = A * x
y = solve(LinearProblem(A, b)).u
sum(y)
1000.0
using ExtendableSparse
using LinearSolve
using MultiFloats
A = fdrand(Float64x2, 10, 10, 10; matrixtype = ExtendableSparseMatrix)
x = ones(Float64x2,1000)
b = A * x
y = solve(LinearProblem(A, b), SparspakFactorization()).u
sum(y)
999.999999999999999999999999999838163512840608

Preconditioned Krylov solvers with LinearSolve.jl

Since version 1.6, preconditioner constructors can be passed to iterative solvers via the precs keyword argument to the iterative solver specification.

using ExtendableSparse
using LinearSolve
using ExtendableSparse: ILUZeroPreconBuilder
A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix)
x = ones(1000)
b = A * x
y = LinearSolve.solve(LinearProblem(A, b),
                      KrylovJL_CG(; precs=ILUZeroPreconBuilder())).u
sum(y)
999.9999999378895

Available preconditioners

ExtendableSparse provides constructors for preconditioners wich can be used as precs. These generally return a tuple (Pl,I) of a left preconditioner and a trivial right preconditioner.

ExtendableSparse has a number of package extensions which construct preconditioners from some other packages. In the future, these packages may provide this functionality on their own.

In addition, ExtendableSparse implements some preconditioners:

ExtendableSparse.JacobiPreconBuilder — Type
JacobiPreconBuilder()

Return callable object constructing a left Jacobi preconditioner to be passed as the precs parameter to iterative methods wrapped by LinearSolve.jl.

source

LU factorizations of matrices from previous iteration steps may be good preconditioners for Krylov solvers called during a nonlinear solve via Newton's method. For this purpose, ExtendableSparse provides a preconditioner constructor which wraps sparse LU factorizations supported by LinearSolve.jl

Block preconditioner constructors are provided as well

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

The example beloww shows how to create a block Jacobi preconditioner where the blocks are defined by even and odd degrees of freedom, and the diagonal blocks are solved using UMFPACK.

using ExtendableSparse
using LinearSolve
using ExtendableSparse: LinearSolvePreconBuilder, BlockPreconBuilder
A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix)
x = ones(1000)
b = A * x
partitioning=A->[1:2:size(A,1), 2:2:size(A,1)]
umfpackprecs=LinearSolvePreconBuilder(UMFPACKFactorization())
blockprecs=BlockPreconBuilder(;precs=umfpackprecs, partitioning)
y = LinearSolve.solve(LinearProblem(A, b), KrylovJL_CG(; precs=blockprecs)).u
sum(y)
999.9999999554125

umpfackpreconbuilder e.g. could be replaced by SmoothedAggregationPreconBuilder(). Moreover, this approach works for any AbstractSparseMatrixCSC.

Deprecated API

Passing a preconditioner via the Pl or Pr keyword arguments will be deprecated in LinearSolve. ExtendableSparse used to export a number of wrappers for preconditioners from other packages for this purpose. This approach is deprecated as of v1.6 and will be removed with v2.0.

using ExtendableSparse
using LinearSolve
using SparseArray
using ILUZero
A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix)
x = ones(1000)
b = A * x
y = LinearSolve.solve(LinearProblem(A, b), KrylovJL_CG();
                      Pl = ILUZero.ilu0(SparseMatrixCSC(A))).u
sum(y)