|
1 | | -# Integration with LinearSolve.jl |
| 1 | +# Linear System solution |
| 2 | + |
| 3 | +## The `\` operator |
| 4 | +The packages overloads `\` for the ExtendableSparseMatrix type. |
| 5 | +The following example uses [`fdrand`](@ref) to create a test matrix and solve |
| 6 | +the corresponding linear system of equations. |
| 7 | + |
| 8 | +```@example |
| 9 | +using ExtendableSparse |
| 10 | +A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix) |
| 11 | +x = ones(1000) |
| 12 | +b = A * x |
| 13 | +y = A \ b |
| 14 | +sum(y) |
| 15 | +``` |
| 16 | + |
| 17 | +This works as well for number types besides `Float64` and related, in this case, |
| 18 | +by default a LU factorization based on Sparspak ist used. |
| 19 | + |
| 20 | +```@example |
| 21 | +using ExtendableSparse |
| 22 | +using MultiFloats |
| 23 | +A = fdrand(Float64x2, 10, 10, 10; matrixtype = ExtendableSparseMatrix) |
| 24 | +x = ones(Float64x2,1000) |
| 25 | +b = A * x |
| 26 | +y = A \ b |
| 27 | +sum(y) |
| 28 | +``` |
| 29 | + |
| 30 | +## Solving with LinearSolve.jl |
2 | 31 |
|
3 | 32 | Starting with version 0.9.6, ExtendableSparse is compatible |
4 | 33 | with [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl). |
5 | 34 | Since version 0.9.7, this is facilitated via the |
6 | | -AbstractSparseMatrixCSC interface. |
| 35 | +AbstractSparseMatrixCSC interface. |
7 | 36 |
|
8 | | -```@autodocs |
9 | | -Modules = [ExtendableSparse] |
10 | | -Pages = ["linearsolve.jl"] |
11 | | -``` |
12 | 37 |
|
13 | | -We can create a test problem and solve it with the `\` operator. |
| 38 | +The same problem can be solved via `LinearSolve.jl`: |
14 | 39 |
|
15 | 40 | ```@example |
16 | | -using ExtendableSparse # hide |
| 41 | +using ExtendableSparse |
| 42 | +using LinearSolve |
17 | 43 | A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix) |
18 | 44 | x = ones(1000) |
19 | 45 | b = A * x |
20 | | -y = A \ b |
| 46 | +y = solve(LinearProblem(A, b)).u |
21 | 47 | sum(y) |
22 | 48 | ``` |
23 | 49 |
|
24 | | -The same problem can be solved by the tools available via `LinearSolve.jl`: |
| 50 | +```@example |
| 51 | +using ExtendableSparse |
| 52 | +using LinearSolve |
| 53 | +using MultiFloats |
| 54 | +A = fdrand(Float64x2, 10, 10, 10; matrixtype = ExtendableSparseMatrix) |
| 55 | +x = ones(Float64x2,1000) |
| 56 | +b = A * x |
| 57 | +y = solve(LinearProblem(A, b), SparspakFactorization()).u |
| 58 | +sum(y) |
| 59 | +``` |
| 60 | + |
| 61 | +## Preconditioned Krylov solvers with LinearSolve.jl |
| 62 | + |
| 63 | +Since version 1.6, preconditioner constructors can be passed to iterative solvers via the [`precs` |
| 64 | +keyword argument](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/#prec) |
| 65 | +to the iterative solver specification. |
25 | 66 |
|
26 | 67 | ```@example |
27 | | -using ExtendableSparse # hide |
28 | | -using LinearSolve # hide |
| 68 | +using ExtendableSparse |
| 69 | +using LinearSolve |
| 70 | +using ExtendableSparse: ILUZeroPreconBuilder |
29 | 71 | A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix) |
30 | 72 | x = ones(1000) |
31 | 73 | b = A * x |
32 | | -y = solve(LinearProblem(A, b), SparspakFactorization()).u |
| 74 | +y = LinearSolve.solve(LinearProblem(A, b), |
| 75 | + KrylovJL_CG(; precs=ILUZeroPreconBuilder())).u |
33 | 76 | sum(y) |
34 | 77 | ``` |
35 | 78 |
|
36 | | -Also, the iterative method interface works with ExtendableSparse. |
| 79 | +## Available preconditioners |
| 80 | +ExtendableSparse provides constructors for preconditioners wich can be used as `precs`. |
| 81 | +These generally return a tuple `(Pl,I)` of a left preconditioner and a trivial right preconditioner. |
| 82 | + |
| 83 | +ExtendableSparse has a number of package extensions which construct preconditioners |
| 84 | +from some other packages. In the future, these packages may provide this functionality on their own. |
| 85 | + |
| 86 | +```@docs |
| 87 | +ExtendableSparse.ILUZeroPreconBuilder |
| 88 | +ExtendableSparse.ILUTPreconBuilder |
| 89 | +ExtendableSparse.SmoothedAggregationPreconBuilder |
| 90 | +ExtendableSparse.RugeStubenPreconBuilder |
| 91 | +``` |
| 92 | + |
| 93 | +In addition, ExtendableSparse implements some preconditioners: |
| 94 | + |
| 95 | +```@docs |
| 96 | +ExtendableSparse.JacobiPreconBuilder |
| 97 | +``` |
| 98 | + |
| 99 | +LU factorizations of matrices from previous iteration steps may be good |
| 100 | +preconditioners for Krylov solvers called during a nonlinear solve via |
| 101 | +Newton's method. For this purpose, ExtendableSparse provides a preconditioner constructor |
| 102 | +which wraps sparse LU factorizations supported by LinearSolve.jl |
| 103 | +```@docs |
| 104 | +ExtendableSparse.LinearSolvePreconBuilder |
| 105 | +``` |
| 106 | + |
| 107 | + |
| 108 | +Block preconditioner constructors are provided as well |
| 109 | +```@docs; canonical=false |
| 110 | +ExtendableSparse.BlockPreconBuilder |
| 111 | +``` |
| 112 | + |
37 | 113 |
|
| 114 | +The example beloww shows how to create a block Jacobi preconditioner where the blocks are defined by even and odd |
| 115 | +degrees of freedom, and the diagonal blocks are solved using UMFPACK. |
38 | 116 | ```@example |
39 | | -using ExtendableSparse # hide |
40 | | -using LinearSolve # hide |
41 | | -using SparseArrays # hide |
42 | | -using ILUZero # hide |
| 117 | +using ExtendableSparse |
| 118 | +using LinearSolve |
| 119 | +using ExtendableSparse: LinearSolvePreconBuilder, BlockPreconBuilder |
43 | 120 | A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix) |
44 | 121 | x = ones(1000) |
45 | 122 | b = A * x |
46 | | -y = LinearSolve.solve(LinearProblem(A, b), KrylovJL_CG(); |
47 | | - Pl = ILUZero.ilu0(SparseMatrixCSC(A))).u |
| 123 | +partitioning=A->[1:2:size(A,1), 2:2:size(A,1)] |
| 124 | +umfpackprecs=LinearSolvePreconBuilder(UMFPACKFactorization()) |
| 125 | +blockprecs=BlockPreconBuilder(;precs=umfpackprecs, partitioning) |
| 126 | +y = LinearSolve.solve(LinearProblem(A, b), KrylovJL_CG(; precs=blockprecs)).u |
48 | 127 | sum(y) |
49 | 128 | ``` |
| 129 | +`umpfackpreconbuilder` e.g. could be replaced by `SmoothedAggregationPreconBuilder()`. Moreover, this approach |
| 130 | +works for any `AbstractSparseMatrixCSC`. |
| 131 | + |
| 132 | + |
| 133 | +## Deprecated API |
| 134 | +Passing a preconditioner via the `Pl` or `Pr` keyword arguments |
| 135 | +will be deprecated in LinearSolve. ExtendableSparse used to |
| 136 | +export a number of wrappers for preconditioners from other packages |
| 137 | +for this purpose. This approach is deprecated as of v1.6 and will be removed |
| 138 | +with v2.0. |
50 | 139 |
|
51 | | -However, ExtendableSparse provides a number of wrappers around preconditioners |
52 | | -from various Julia packages. |
53 | 140 | ```@example |
54 | | -using ExtendableSparse # hide |
55 | | -using LinearSolve # hide |
56 | | -using ILUZero # hide |
| 141 | +using ExtendableSparse |
| 142 | +using LinearSolve |
| 143 | +using SparseArray |
| 144 | +using ILUZero |
57 | 145 | A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix) |
58 | 146 | x = ones(1000) |
59 | 147 | b = A * x |
60 | 148 | y = LinearSolve.solve(LinearProblem(A, b), KrylovJL_CG(); |
61 | | - Pl = ILUZeroPreconditioner(A)).u |
| 149 | + Pl = ILUZero.ilu0(SparseMatrixCSC(A))).u |
62 | 150 | sum(y) |
63 | 151 | ``` |
| 152 | + |
0 commit comments