From bf6221ece445373b3abb9ad779aaa769d513e945 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Mon, 21 Oct 2024 21:59:49 +0200 Subject: [PATCH 1/6] Document `precs` --- docs/src/basics/Preconditioners.md | 50 ++++++++++++++++++++++++++---- 1 file changed, 44 insertions(+), 6 deletions(-) diff --git a/docs/src/basics/Preconditioners.md b/docs/src/basics/Preconditioners.md index 05faedb06..d4f333548 100644 --- a/docs/src/basics/Preconditioners.md +++ b/docs/src/basics/Preconditioners.md @@ -4,15 +4,13 @@ Many linear solvers can be accelerated by using what is known as a **preconditio an approximation to the matrix inverse action which is cheap to evaluate. These can improve the numerical conditioning of the solver process and in turn improve the performance. LinearSolve.jl provides an interface for the definition of -preconditioners which works with the wrapped packages. +preconditioners which works with the wrapped iterative solver packages. ## Using Preconditioners ### Mathematical Definition -Preconditioners are specified in the keyword arguments to `init` or `solve`: `Pl` for left -and `Pr` for right preconditioner, respectively. -The right preconditioner, ``P_r`` transforms the linear system ``Au = b`` into the form: +A right preconditioner, ``P_r`` transforms the linear system ``Au = b`` into the form: ```math AP_r^{-1}(P_r u) = AP_r^{-1}y = b @@ -31,10 +29,12 @@ A two-sided preconditioned system is of the form: P_l^{-1}A P_r^{-1} (P_r u) = P_l^{-1}b ``` -By default, if no preconditioner is given, the preconditioner is assumed to be +### Specifying Preconditioners + +One way to specify preconditioners uses the `Pl` and `Pr` keyword arguments to `init` or `solve`: `Pl` for left +and `Pr` for right preconditioner, respectively. By default, if no preconditioner is given, the preconditioner is assumed to be the identity ``I``. -### Using Preconditioners In the following, we will use the `DiagonalPreconditioner` to define a two-sided preconditioned system which first divides by some random numbers and then @@ -55,6 +55,44 @@ sol = solve(prob, KrylovJL_GMRES(), Pl = Pl) sol.u ``` +Alternatively, preconditioners can be specified via the `precs` argument to the constructor of +an iterative solver specification. This argument shall deliver a function mapping `A` and a +parameter `p` to a tuple `(Pl,Pr)` consisting a left and a right preconditioner. + + +```@example precon2 +using LinearSolve, LinearAlgebra +n = 4 +s = rand(n) + +A = rand(n, n) +b = rand(n) + +prob = LinearProblem(A, b) +sol = solve(prob, KrylovJL_GMRES(precs = (A,p)->(Diagonal(A),I)) ) +sol.u +``` +This approach has the advantage that the specification of the preconditioner is possible without +the knowledge of a concrete matrix `A`. It also allows to specifiy the preconditioner via a functor struct: + +```@example precon2 +using LinearSolve, LinearAlgebra + +struct DiagonalPrecs end + +(::DiagonalPrecs)(A,p) = (Diagonal(A),I) + +n = 4 +s = rand(n) + +A = rand(n, n) +b = rand(n) + +prob = LinearProblem(A, b) +sol = solve(prob, KrylovJL_GMRES(precs = DiagonalPrecs()) ) +sol.u +``` + ## Preconditioner Interface To define a new preconditioner you define a Julia type which satisfies the From 89ed45b74c1e523ced1794f6a5faee36cfb8400d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 22 Oct 2024 09:19:15 +0200 Subject: [PATCH 2/6] fix typo; replace "functor" by "callable object" Googling "Callable object" leads to the right place in the julia docs and is easier to understand --- docs/src/basics/Preconditioners.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/basics/Preconditioners.md b/docs/src/basics/Preconditioners.md index d4f333548..9fa7e4c02 100644 --- a/docs/src/basics/Preconditioners.md +++ b/docs/src/basics/Preconditioners.md @@ -73,7 +73,7 @@ sol = solve(prob, KrylovJL_GMRES(precs = (A,p)->(Diagonal(A),I)) ) sol.u ``` This approach has the advantage that the specification of the preconditioner is possible without -the knowledge of a concrete matrix `A`. It also allows to specifiy the preconditioner via a functor struct: +the knowledge of a concrete matrix `A`. It also allows to specify the preconditioner via a callable object: ```@example precon2 using LinearSolve, LinearAlgebra From 8b8b1714b08f3ebccc76e1085aeac91f01a66843 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 22 Oct 2024 23:52:08 +0200 Subject: [PATCH 3/6] Provide a more thorough explanation of precs, including reinit! with reuse_precs=true. --- docs/Project.toml | 1 + docs/src/basics/Preconditioners.md | 39 ++++++++++++++++-------------- 2 files changed, 22 insertions(+), 18 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 5934dd13c..4900c58df 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,5 +1,6 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" [compat] diff --git a/docs/src/basics/Preconditioners.md b/docs/src/basics/Preconditioners.md index 9fa7e4c02..7b2919243 100644 --- a/docs/src/basics/Preconditioners.md +++ b/docs/src/basics/Preconditioners.md @@ -36,34 +36,30 @@ and `Pr` for right preconditioner, respectively. By default, if no preconditione the identity ``I``. -In the following, we will use the `DiagonalPreconditioner` to define a two-sided -preconditioned system which first divides by some random numbers and then -multiplies by the same values. This is commonly used in the case where if, instead -of random, `s` is an approximation to the eigenvalues of a system. +In the following, we will use a left sided diagonal (Jacobi) preconditioner. -```@example precon +```@example precon1 using LinearSolve, LinearAlgebra n = 4 -s = rand(n) -Pl = Diagonal(s) A = rand(n, n) b = rand(n) +Pl=Diagonal(A) + prob = LinearProblem(A, b) sol = solve(prob, KrylovJL_GMRES(), Pl = Pl) sol.u ``` Alternatively, preconditioners can be specified via the `precs` argument to the constructor of -an iterative solver specification. This argument shall deliver a function mapping `A` and a +an iterative solver specification. This argument shall deliver a factory method mapping `A` and a parameter `p` to a tuple `(Pl,Pr)` consisting a left and a right preconditioner. ```@example precon2 using LinearSolve, LinearAlgebra n = 4 -s = rand(n) A = rand(n, n) b = rand(n) @@ -73,26 +69,33 @@ sol = solve(prob, KrylovJL_GMRES(precs = (A,p)->(Diagonal(A),I)) ) sol.u ``` This approach has the advantage that the specification of the preconditioner is possible without -the knowledge of a concrete matrix `A`. It also allows to specify the preconditioner via a callable object: +the knowledge of a concrete matrix `A`. It also allows to specify the preconditioner via a callable object +and to pass parameters to the constructor of the preconditioner instances. The example below also shows how +to reuse the preconditioner once constructed for the subsequent solution of a modified problem. -```@example precon2 +```@example precon3 using LinearSolve, LinearAlgebra -struct DiagonalPrecs end +Base.@kwdef struct WeightedDiagonalBuilder + w::Float64 +end -(::DiagonalPrecs)(A,p) = (Diagonal(A),I) +(builder::WeightedDiagonalBuilder)(A,p) = (builder.w*Diagonal(A),I) n = 4 -s = rand(n) - -A = rand(n, n) +A = n*I-rand(n, n) b = rand(n) prob = LinearProblem(A, b) -sol = solve(prob, KrylovJL_GMRES(precs = DiagonalPrecs()) ) +sol = solve(prob, KrylovJL_GMRES(precs = WeightedDiagonalBuilder(w=0.9)) ) sol.u -``` +B=A.+0.1 +cache=sol.cache +reinit!(cache,A=B, reuse_precs=true) +sol = solve!(cache, KrylovJL_GMRES(precs = WeightedDiagonalBuilder(w=0.9)) ) +sol.u +``` ## Preconditioner Interface To define a new preconditioner you define a Julia type which satisfies the From 3cb97cc05fd3fbde607a4a5136cd919695e5f340 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Wed, 23 Oct 2024 23:06:28 +0200 Subject: [PATCH 4/6] *Builder -> *PreconBuilder in proposed pattern --- docs/src/basics/Preconditioners.md | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/src/basics/Preconditioners.md b/docs/src/basics/Preconditioners.md index 7b2919243..28304731f 100644 --- a/docs/src/basics/Preconditioners.md +++ b/docs/src/basics/Preconditioners.md @@ -76,24 +76,24 @@ to reuse the preconditioner once constructed for the subsequent solution of a mo ```@example precon3 using LinearSolve, LinearAlgebra -Base.@kwdef struct WeightedDiagonalBuilder +Base.@kwdef struct WeightedDiagonalPreconBuilder w::Float64 end -(builder::WeightedDiagonalBuilder)(A,p) = (builder.w*Diagonal(A),I) +(builder::WeightedDiagonalPreconBuilder)(A,p) = (builder.w*Diagonal(A),I) n = 4 A = n*I-rand(n, n) b = rand(n) prob = LinearProblem(A, b) -sol = solve(prob, KrylovJL_GMRES(precs = WeightedDiagonalBuilder(w=0.9)) ) +sol = solve(prob, KrylovJL_GMRES(precs = WeightedDiagonalPreconBuilder(w=0.9)) ) sol.u B=A.+0.1 cache=sol.cache reinit!(cache,A=B, reuse_precs=true) -sol = solve!(cache, KrylovJL_GMRES(precs = WeightedDiagonalBuilder(w=0.9)) ) +sol = solve!(cache, KrylovJL_GMRES(precs = WeightedDiagonalPreconBuilder(w=0.9)) ) sol.u ``` ## Preconditioner Interface @@ -128,14 +128,14 @@ The following preconditioners match the interface of LinearSolve.jl. Implementations of the algebraic multigrid method. Must be converted to a preconditioner via `AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.precmethod(A))`. Requires `A` as a `AbstractMatrix`. Provides the following methods: - + + `AlgebraicMultigrid.ruge_stuben(A)` + `AlgebraicMultigrid.smoothed_aggregation(A)` - [PyAMG](https://github.com/cortner/PyAMG.jl): Implementations of the algebraic multigrid method. Must be converted to a preconditioner via `PyAMG.aspreconditioner(PyAMG.precmethod(A))`. Requires `A` as a `AbstractMatrix`. Provides the following methods: - + + `PyAMG.RugeStubenSolver(A)` + `PyAMG.SmoothedAggregationSolver(A)` - [ILUZero.ILU0Precon(A::SparseMatrixCSC{T,N}, b_type = T)](https://github.com/mcovalt/ILUZero.jl): @@ -154,7 +154,7 @@ The following preconditioners match the interface of LinearSolve.jl. and `HYPRE.BoomerAMG`. - [KrylovPreconditioners.jl](https://github.com/JuliaSmoothOptimizers/KrylovPreconditioners.jl/): Provides GPU-ready preconditioners via KernelAbstractions.jl. At the time of writing the package provides the following methods: - + + Incomplete Cholesky decomposition `KrylovPreconditioners.kp_ic0(A)` + Incomplete LU decomposition `KrylovPreconditioners.kp_ilu0(A)` + Block Jacobi `KrylovPreconditioners.BlockJacobiPreconditioner(A, nblocks, device)` From 350f26b57058a4963ca2946268e0973830c69042 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Wed, 30 Oct 2024 21:06:43 +0100 Subject: [PATCH 5/6] remove IncompleteLU from docs project --- docs/Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index 4900c58df..5934dd13c 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,5 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" [compat] From bd0d087c4c2338a40be65811b675cb64251e6839 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Mon, 4 Nov 2024 20:02:49 +0100 Subject: [PATCH 6/6] ... format --- docs/src/basics/Preconditioners.md | 28 ++++++++++++++-------------- test/resolve.jl | 2 +- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/docs/src/basics/Preconditioners.md b/docs/src/basics/Preconditioners.md index 28304731f..d2148de54 100644 --- a/docs/src/basics/Preconditioners.md +++ b/docs/src/basics/Preconditioners.md @@ -35,7 +35,6 @@ One way to specify preconditioners uses the `Pl` and `Pr` keyword arguments to and `Pr` for right preconditioner, respectively. By default, if no preconditioner is given, the preconditioner is assumed to be the identity ``I``. - In the following, we will use a left sided diagonal (Jacobi) preconditioner. ```@example precon1 @@ -45,7 +44,7 @@ n = 4 A = rand(n, n) b = rand(n) -Pl=Diagonal(A) +Pl = Diagonal(A) prob = LinearProblem(A, b) sol = solve(prob, KrylovJL_GMRES(), Pl = Pl) @@ -56,7 +55,6 @@ Alternatively, preconditioners can be specified via the `precs` argument to th an iterative solver specification. This argument shall deliver a factory method mapping `A` and a parameter `p` to a tuple `(Pl,Pr)` consisting a left and a right preconditioner. - ```@example precon2 using LinearSolve, LinearAlgebra n = 4 @@ -65,9 +63,10 @@ A = rand(n, n) b = rand(n) prob = LinearProblem(A, b) -sol = solve(prob, KrylovJL_GMRES(precs = (A,p)->(Diagonal(A),I)) ) +sol = solve(prob, KrylovJL_GMRES(precs = (A, p) -> (Diagonal(A), I))) sol.u ``` + This approach has the advantage that the specification of the preconditioner is possible without the knowledge of a concrete matrix `A`. It also allows to specify the preconditioner via a callable object and to pass parameters to the constructor of the preconditioner instances. The example below also shows how @@ -80,22 +79,23 @@ Base.@kwdef struct WeightedDiagonalPreconBuilder w::Float64 end -(builder::WeightedDiagonalPreconBuilder)(A,p) = (builder.w*Diagonal(A),I) +(builder::WeightedDiagonalPreconBuilder)(A, p) = (builder.w * Diagonal(A), I) n = 4 -A = n*I-rand(n, n) +A = n * I - rand(n, n) b = rand(n) prob = LinearProblem(A, b) -sol = solve(prob, KrylovJL_GMRES(precs = WeightedDiagonalPreconBuilder(w=0.9)) ) +sol = solve(prob, KrylovJL_GMRES(precs = WeightedDiagonalPreconBuilder(w = 0.9))) sol.u -B=A.+0.1 -cache=sol.cache -reinit!(cache,A=B, reuse_precs=true) -sol = solve!(cache, KrylovJL_GMRES(precs = WeightedDiagonalPreconBuilder(w=0.9)) ) +B = A .+ 0.1 +cache = sol.cache +reinit!(cache, A = B, reuse_precs = true) +sol = solve!(cache, KrylovJL_GMRES(precs = WeightedDiagonalPreconBuilder(w = 0.9))) sol.u ``` + ## Preconditioner Interface To define a new preconditioner you define a Julia type which satisfies the @@ -128,14 +128,14 @@ The following preconditioners match the interface of LinearSolve.jl. Implementations of the algebraic multigrid method. Must be converted to a preconditioner via `AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.precmethod(A))`. Requires `A` as a `AbstractMatrix`. Provides the following methods: - + + `AlgebraicMultigrid.ruge_stuben(A)` + `AlgebraicMultigrid.smoothed_aggregation(A)` - [PyAMG](https://github.com/cortner/PyAMG.jl): Implementations of the algebraic multigrid method. Must be converted to a preconditioner via `PyAMG.aspreconditioner(PyAMG.precmethod(A))`. Requires `A` as a `AbstractMatrix`. Provides the following methods: - + + `PyAMG.RugeStubenSolver(A)` + `PyAMG.SmoothedAggregationSolver(A)` - [ILUZero.ILU0Precon(A::SparseMatrixCSC{T,N}, b_type = T)](https://github.com/mcovalt/ILUZero.jl): @@ -154,7 +154,7 @@ The following preconditioners match the interface of LinearSolve.jl. and `HYPRE.BoomerAMG`. - [KrylovPreconditioners.jl](https://github.com/JuliaSmoothOptimizers/KrylovPreconditioners.jl/): Provides GPU-ready preconditioners via KernelAbstractions.jl. At the time of writing the package provides the following methods: - + + Incomplete Cholesky decomposition `KrylovPreconditioners.kp_ic0(A)` + Incomplete LU decomposition `KrylovPreconditioners.kp_ilu0(A)` + Block Jacobi `KrylovPreconditioners.BlockJacobiPreconditioner(A, nblocks, device)` diff --git a/test/resolve.jl b/test/resolve.jl index 5d0beb343..892401087 100644 --- a/test/resolve.jl +++ b/test/resolve.jl @@ -2,7 +2,7 @@ using LinearSolve, LinearAlgebra, SparseArrays, InteractiveUtils, Test using LinearSolve: AbstractDenseFactorization, AbstractSparseFactorization for alg in vcat(InteractiveUtils.subtypes(AbstractDenseFactorization), - InteractiveUtils.subtypes(AbstractSparseFactorization)) + InteractiveUtils.subtypes(AbstractSparseFactorization)) if alg in [PardisoJL] ## Pardiso has extra tests in test/pardiso/pardiso.jl continue