# Defining a preconditioner for IterativeSolvers
**URL:** https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977
**Category:** Numerics
**Tags:** iterative-solvers
**Created:** [September 9, 2022, 9:33am UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977 "2022-09-09T09:33:04Z")
**Posts on this page:** 19
**Page:** 1
### Author:  [@Gravlax](https://discourse.julialang.org/u/Gravlax)
#### Post date: [September 9, 2022, 9:33am UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/1 "2022-09-09T09:33:04Z")
Dear all,
I would like to solve a linear(ised) system of equations using IterativeSolvers, however I am not sure how to define the preconditioner.
The preconditonner is a symmetric positive-definite sparse matrix and the action of the preconditonner is defined in 2 steps: (1) compute Cholesky factorisation of the preconditionner (done before calling the iterative solver), (2) apply Cholesky factors during the iterative solver (backsubstitutions).
I understand that the action of the preconditionner is dealt internally with in-place backsubstitutions, which is great.
But how should I input the Cholesky factors (SuiteSparse.CHOLMOD.Factor{Float64}) to say, `gmres!` or `bicgstabl!`? Can I simply do this?:
```julia
PC_fact = cholesky(PC) # Factor preconditionner
gmres!(x, A, b; Pl=PC_fact)
```
So far I use a self programmed `gcr!` translated from an earlier MATLAB code (where I actually hard-coded the action of the preconditionner) but I am confident that I’d be better off using solvers from IterativeSolvers.jl 😃
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [September 9, 2022, 11:10am UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/2 "2022-09-09T11:10:27Z")
> [@Gravlax](#):
>
> Can I simply do this?
Yes, I think so. According to [the docs](https://iterativesolvers.julialinearalgebra.org/dev/preconditioning/), the preconditioner `P` is used by IterativeSolvers.jl via `P \ x` etcetera, and [`cholesky(A)` in Julia](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.cholesky) returns a factorization object that supports `\` and `ldiv!` (out-of-place and in-place left division), implemented via back/forward-substitution.
---
### Author:  [@j-fu](https://discourse.julialang.org/u/j-fu)
#### Post date: [September 9, 2022, 3:18pm UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/3 "2022-09-09T15:18:21Z")
Hi, this probably would work, but `cholesky(A)` is a full factorization. For a Preconditioner see e.g. IncompleteLU.jl or AlgebraicMultigrid.jl.
The point is that a preconditioner M for a matrix A is much simpler to factorize than the full matrix A. In practice one often doesn’t represent M, but just the factorization.
Also, cholesky is for symmetric matrices, in that case you can use cg instead of gmres.
I have some material on this:
[https://www.wias-berlin.de/people/fuhrmann/SciComp-WS2122/week4/](https://www.wias-berlin.de/people/fuhrmann/SciComp-WS2122/week4/)
For a new take on sparse system solution which incorporates all these concepts and packages see LinearSolve.jl.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [September 9, 2022, 4:08pm UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/4 "2022-09-09T16:08:23Z")
> [@j-fu](#):
>
> Hi, this probably would work, but `cholesky(A)` is a full factorization. For a Preconditioner see e.g. IncompleteLU.jl or AlgebraicMultigrid.jl.
My understanding is that the original poster is factorizing some matrix P _different_ from the original matrix A. It can be very sensible to use the _full_ factorization of P in this case if P is much cheaper to factorize than A (e.g. P is banded).
---
### Author:  [@Gravlax](https://discourse.julialang.org/u/Gravlax)
#### Post date: [September 15, 2022, 5:00pm UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/5 "2022-09-15T17:00:28Z")
Thanks for both of your inputs!
@steveng : exact, I’m solving for a Newton step but the Jacobian is not symmetric. Hence, I use the Cholesky factorisation of the system resulting from Picard linearisation (symmetric positive definite). To my knowledge, it’s an efficient way to solve 2D non-linear problems (FDM/FEM). Not a candidate for 3D though…
@j-fu : thanks for the suggestion, I’ll definitely try the available ILU and AMG once I’m back at work.
---
### Author:  [@amontoison](https://discourse.julialang.org/u/amontoison)
#### Post date: [September 15, 2022, 9:04pm UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/6 "2022-09-15T21:04:36Z")
@Gravlax `ldiv!(y, F, x)` is not implemented for `F::SuiteSparse.CHOLMOD.Factor{Float64}`.
You can use [LDLFactorizations.jl](https://github.com/JuliaSmoothOptimizers/LDLFactorizations.jl), if you want to use `ldiv!` with the factorization of a symmetric (positive definite) matrix.
I give some examples with preconditioners in my package [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) if it can help you:
[https://juliasmoothoptimizers.github.io/Krylov.jl/stable/preconditioners/#Examples](https://juliasmoothoptimizers.github.io/Krylov.jl/stable/preconditioners/#Examples)
The examples should also work with IterativeSolvers.jl.
---
### Author:  [@Gravlax](https://discourse.julialang.org/u/Gravlax)
#### Post date: [September 16, 2022, 2:10pm UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/7 "2022-09-16T14:10:32Z")
Oh nice, that’s indeed very useful. If I understand well, you manage to overcome the `ldiv!`/cholesky factor issue by “overloading” or adding a method to ‘ldiv!’ with the following line, right?
```julia
ldiv!(y::Vector{T}, F::SuiteSparse.CHOLMOD.Factor{T}, x::Vector{T}) where T = (y .= F \ x)
```
---
### Author:  [@amontoison](https://discourse.julialang.org/u/amontoison)
#### Post date: [September 16, 2022, 6:09pm UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/8 "2022-09-16T18:09:26Z")
Yes, you’re right.
Note that we allocate a new vector at each call of `ldiv!` for `F \ x`.
It’s not a “real” ldiv! routine.
We can’t do better because the in-place backward and forward sweeps with permutations are not implemented in `CHOLMOD`…
---
### Author:  [@vlc](https://discourse.julialang.org/u/vlc)
#### Post date: [September 16, 2022, 6:33pm UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/9 "2022-09-16T18:33:52Z")
Hi there, I’ve ran into a similar issue recently, but even factorisations that do provide in-place `ldiv!` method (such as UMFPACK’s LU decomposition) turn out to allocate quite a lot. However looking into UMFPACK’s documentation, it looks like there it defines a set of functions (umfpack\_\*\_wsolve) that takes a pointer to a preallocated buffer from the user (and therefore does not allocate itself).
I’m curious as to how hard it would be to expose this feature in the Julia wrapper, much like what [FastLapackInterface](https://github.com/DynareJulia/FastLapackInterface.jl) does for LAPACK (only problem is that the matrix formats LAPACK supports is quite limited…).
---
### Author:  [@amontoison](https://discourse.julialang.org/u/amontoison)
#### Post date: [September 16, 2022, 8:15pm UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/10 "2022-09-16T20:15:02Z")
@vlc
>
>
> The method \`ldiv!\` has very large allocations happening for sparse LU factorizat…ions. See the example below:
> \`\`\`julia
> using SparseArrays
> using LinearAlgebra
>
> N = 1000
>
> A = sparse(I, N, N) + sprand(N, N, .01)
> fact = lu(A)
> b = rand(N)
> x = zeros(N)
>
> @time ldiv!(x, fact, b)
>
> \#Convert to a full matrix and recompute the factorization
> A = Matrix(A)
> fact = lu(A)
>
> @time ldiv!(x, fact, b);
> \`\`\`
> A characteristic test runs gives
> \`\`\`
> 0.001000 seconds (2 allocations: 46.906 KiB)
> 0.000669 seconds (2 allocations: 96 bytes)
> \`\`\`
>
> The case of the sparse LU decomposition causes a large amount of memory allocation in \`ldiv!\` when I would expect it to no different from the case of a dense matrix.
>
> This is on Julia Version 1.6.1.
I remember that someone fixed this problem recently.
>
>
> closes JuliaSparse/SparseArrays.jl#112, passes all tests (and some new). The onl…y possible painpoint is that it makes \`ldvi!(x, l, b)\` edit l and makes it thread unsafe. I could add a a lock...
We will have this problem solved with the future releases of Julia.
---
### Author:  [@vlc](https://discourse.julialang.org/u/vlc)
#### Post date: [September 18, 2022, 5:14pm UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/11 "2022-09-18T17:14:26Z")
Thanks @amontoison that’s exactly what I was looking for!
---
### Author:  [@amontoison](https://discourse.julialang.org/u/amontoison)
#### Post date: [September 20, 2022, 11:01pm UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/12 "2022-09-20T23:01:40Z")
I opened an issue about the problem with CHOLMOD in [DrTimothyAldenDavis/SuiteSparse](https://github.com/DrTimothyAldenDavis/SuiteSparse/issues/136).
---
### Author:  [@Gravlax](https://discourse.julialang.org/u/Gravlax)
#### Post date: [September 21, 2022, 6:06pm UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/13 "2022-09-21T18:06:14Z")
Great, thanks a lot for this. I did not understand that the problem was inherent to CHOLMOD. I thought it was related the way it was wrapped in Julia.
That said, it is true that when using the library in C:
` s = cholmod_solve (CHOLMOD_A, Lfact, f, c);`
`s` is internally allocated. So, it’s probably a similar issue.
---
### Author:  [@tduretz](https://discourse.julialang.org/u/tduretz)
#### Post date: [December 20, 2024, 7:04pm UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/14 "2024-12-20T19:04:20Z")
Hi there,
I am also wondering if there is an easy way to use a precomputed Cholesky factorisation as preconditioner to `gcr!`.
It seems that the `ldiv!` issue is still going on, is there a simple workaround to enable the use of a Cholesky factorisation as preconditioner?
thanks!
---
### Author:  [@tduretz](https://discourse.julialang.org/u/tduretz)
#### Post date: [November 18, 2025, 9:42am UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/15 "2025-11-18T09:42:49Z")
This is indeed working in 1.12.1. Below is a MWE, in case it can be useful to anyone.
This should now allow to implement solvers like those used to make Fig 4 in this [study](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2016GC006727), cool !
```julia
using IterativeSolvers, SparseArrays, LinearAlgebra
let
# System matrix
# Non symmetric for whatever reason
M = [1.0 0.0 0.0 0.0 0.0;
0.0 2.6 -1.0 0.0 0.0;
0.0 -0.5 6.1 -1.9 0.0;
0.0 0.0 -1.3 7.0 0.0;
0.0 0.0 0.0 0.0 1.0]
# Preconditioner
# Symmetric postive definite (e.g.: Picard preconditioner for a Newton solve)
N = [1.0 0.0 0.0 0.0 0.0;
0.0 2.0 -1.0 0.0 0.0;
0.0 -1.0 2.0 -1.0 0.0;
0.0 0.0 -1.0 2.0 0.0;
0.0 0.0 0.0 0.0 1.0]
# Right-hand side
b = [1.0; 2.0; 2.5; 1.5; 4.0]
# Cholesky factors
N_chol = cholesky(N)
# N_chol = lu(M) # This leads to the correct solution in 1 GMRES iteration but it's more expensive than cholesky(N)
# Allocate solution arrays
x = zeros(size(b))
# Call GMRES with cholesky factors as preconditionner
gmres!(x, M, b; Pl=N_chol, verbose=true) # internally uses ldiv!(M_chol, b)
# Check that's it's similar to M⁻¹b
@show x .- M\b
end
```
---
### Author:  [@tduretz](https://discourse.julialang.org/u/tduretz)
#### Post date: [March 18, 2026, 11:38am UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/16 "2026-03-18T11:38:38Z")
In fact this only works if matrices are dense.
So the right preconditioner can be a Cholesky factorisation if:
`typeof(N_chol) = Cholesky{Float64, Matrix{Float64}} `
but not if
`typeof(N_chol) = SparseMatrixCSC{Float64, Int64}`
In the MWE above adding the line `N = sparse(N)` is sufficient to make it break, with the message that `ldiv!()` does not support this type for the factorization.
Is there a why to get it working, else than using dense matrices?
---
### Author:  [@ziolai](https://discourse.julialang.org/u/ziolai)
#### Post date: [March 19, 2026, 11:21am UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/17 "2026-03-19T11:21:53Z")
Replace `ldiv!()` by merely `\` (backslash), assuming here that backslash recognizes the (lower or upper) triangular nature of the matrix?
Alternatively, writing the lower or upper triangular solve as an explicit `for` loop.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 19, 2026, 11:25am UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/18 "2026-03-19T11:25:41Z")
> [@amontoison](#):
>
> `ldiv!(y, F, x)` is not implemented for `F::SuiteSparse.CHOLMOD.Factor{Float64}`.
Update: `ldiv!` was implemented for `CHOLMOD` (sparse Cholesky) factors in 2024 (released in Julia 1.12): [implement in-place `ldiv!` for Cholesky factorization by ranocha · Pull Request #547 · JuliaSparse/SparseArrays.jl · GitHub](https://github.com/JuliaSparse/SparseArrays.jl/pull/547) and corresponding updates were recently incorporated into LinearSolve.jl ([Clean up ldiv!() overloads by JamesWrigley · Pull Request #910 · SciML/LinearSolve.jl · GitHub](https://github.com/SciML/LinearSolve.jl/pull/910))
---
### Author:  [@tduretz](https://discourse.julialang.org/u/tduretz)
#### Post date: [March 20, 2026, 3:10pm UTC](https://discourse.julialang.org/t/defining-a-preconditioner-for-iterativesolvers/86977/19 "2026-03-20T15:10:30Z")
Yes, in-place `ldiv` for sparse Cholesky is the reason why I updated to 1.12. I thought this would work out of the box, and why I got confused. Now I understand this functionality will come with the next release of LinearSolve, thanks for the hint !