Hello,

I am interested in solving the linear system M x = b with an iterative method, where M is positive definite (PD) with the following block structure:

M = \begin{bmatrix}
\Sigma_a & 0\\
0 & \Sigma_b
\end{bmatrix} + \begin{bmatrix}A \Sigma_c A^\top & A \Sigma_c B^\top \\
B \Sigma_c A^\top & B \Sigma_c B^\top \end{bmatrix} \\ = \begin{bmatrix}
\Sigma_a & 0\\
0 & \Sigma_b
\end{bmatrix} + \begin{bmatrix} A \\ B \end{bmatrix} \Sigma_c \begin{bmatrix}A \\ B \end{bmatrix}^\top

where

- A and B are sparse linear operators (e.g. discrete differential operators).
- \Sigma_a and \Sigma_b are PD matrices, such that u \mapsto \Sigma_{a}^{-1} u and v \mapsto \Sigma_{b}^{-1} v are easy to perform.
- \Sigma_c is only positive semi-definite for which only w \mapsto \Sigma_{c} w can be performed.

I am leaning towards a conjugate gradient iterative solver.

I am wondering if people have experience with this kind of matrix, in particular in the design of preconditioners.

Thank you for your help,

You can perhaps try the preconditionner P = diag(Sigma_a^-1, Sigma_b^1)

As @rveltz suggested, the block Jacobi preconditioner `diag(Sigma_a^{-1}, Sigma_b^{-1})`

is a good idea.

I don’t see many alternatives because other blocks are matrix-free operators.

If you already have `M`

as a matrix-free operator, I suggest to use `BlockDiagonalOperator`

from LinearOperators.jl:

You can easily combine your operators to build a preconditioner `P = BlockDiagonalOperator(P1, P2)`

.

To build P1 and P2, you can use `cholesky`

from SuiteSparse or LDLFactorization.jl.

I give some examples of them in the documentation of Krylov.jl:

Your operators P1 and P2 can just perform triangular solves with the factors of `Sigma_a`

and `Sigma_b`

.

If.you use Krylov.jl fo the Krylov method, I suggest the following methods: CG, CR, CAR or MINRES.

```
using Krylov
x, stats = cg(M, b, M=P)
```

1 Like