Preconditioner for scalar laplacian system

Suppose that I have a large (on the order of 10^6 unknowns) 3D scalar Poisson system which I would like to precondition. The boundary conditions have been treated so that the system is SPD. I.e.,

$$\nabla\cdot(\kappa(x)\nabla u(x)) = f$$

the matrices tend to have the following structure (this is a very small problem, the larger problems are indeed sparse):

Is there a consensus about how to precondition such systems (ideally in a way that can be implemented in parallel)?

The linear system is SPD, so it seems natural to choose a PCG scheme, but I haven't been able to find a consensus about how to choose the preconditioner, since many papers deal with more complicated Stokes-like systems instead.

• Choosing the right preconditioner is more of an art than an actual science. If you are dealing with a specific set of equations, then I recommend that you provide the matrices. Otherwise, I am afraid, that it will be hard to give you more than general pointers, say, try an incomplete Cholesky factorization or attempt to reorder the system and extract a narrow banded preconditioner, etc. Moreover, it is entirely possible that a direct solver is applicable. Jun 12, 2017 at 19:31
• I have updated the question. Is it the case that computational scientists typically develop their own preconditioners for such problems? I simply don't want to skip over a resource that might provide a satisfactory answer to the problem if it's well understood. Jun 12, 2017 at 19:47
• What are the properties of the coefficient $\kappa$? Does it vary strongly?
– cfh
Jun 13, 2017 at 6:35
• $\kappa$ varies in space with small differences compared to that of the solution. Jun 13, 2017 at 12:35

2 Answers

Preconditioning and iterative solvers are cool, but did you try to solve your problem with some kind of sparse direct solver? If not, try it first.

State–of–the–art preconditioning techniques for this kind of problems are multigrid techniques. In many cases multigrid allows you to solve your (elliptic) problem in $O(n)$ time.

1. If you don’t want to go deeper, a good start is probably algebraic multigrid. You may try it both as a stand–alone solver and a preconditioner for CG. It is implemented e.g. in PETSc.

2. If you do want to go deeper, a good start is this paper. It has a nice collection of examples revealing multigrid’s behavior depending on problem’s geometry, parameters $f$ and $\kappa$, and so forth.

With SPD matrix, call it $A$, PCG is a good method.

The kind of preconditioner depend a lot by the matrix structure (read kind of matrix) that you have got and there is not a unique answer.

The basic criteria to choose the preconditioner, $P$, are more than one. Considering the PCG the three most important are:

• $AP^{-1} \approx I$
• $P$ is cheap to built
• calculated $P^{-1}r_{k+1}$ is cheap

In general for problem like laplacian you can try the two follow.

Diagonal $D^{-1}$ The diagonal preconditioner is very easy to implement and cheap to apply, but is not very efficient. The matrix contains the verse of $A$ diagonal terms. Better if it is used with matrix $A$ diagonally dominant, but when the extra diagonal terms are big you can find problem.

Incomplete Cholesky $\tilde{L}\tilde{L^T}$ Whit $A$ SPD is defined the Cholesky decomposition. The idea of this preconditioner is to use the decomposition but to conserve the sparsity pattern of $A$, i.e you use the term calculated by the decomposition only for the terms of $A$ different from zero. The existence of this incomplete decomposition depend by the structure of $A$, for example if $A$ is a $M$ matrix it is ok. This kind of preconditioner can improve if you reorder the nodes of $A$ so that the mustache of $A$ are shorter (here you can see an example that the structure of $A$ is important). This preconditioned is apply with the method forward and back (as for the $LU$ decomposition)

When you talk about parallel execution here is important to know the kind of hardware. Anyway you must parallelize the preconditioner application, for example with $D^{-1}$ you parallelize the $\frac{v_i}{d_i}$ operation for the whole vector.

Since you mentioned the Stokes's problem I add some rows about parallelition. If your goal is to use heavy parallelization as gpu is important to use preconditioner with an easy application for gpu, as the product matrix vector sparse (SpMV). In the article of Li and Saad [1] there is a comparison of various preconditioners for gpu.

In the case of gpu polynomial preconditioner is good. I personally use this kind with good performance respect other traditional gpu preconditioner, see this link to CUSP group where I give some detail, with reference, and also the code with some result test. This is a least square preconditioner, it uses Chebyshev polynomials so explicit quadrature formula are not involved, detail are in the article by J. Erhel, F. Guyomarc, and Y. Saad [2].

[1] Ruipeng Li and Yousef Saad. Gpu-accelerated preconditioned iterative linear solvers. The Journal of Supercomputing, 63(2):443–466, 2013.

[2] J. Erhel, F. Guyomarc, and Y. Saad, Least-squares polynomial filters for ill-conditioned lin- ear systems, Tech. Report umsi-2001-32, Minnesota Supercomputer Institute, University of Minnesota, Minneapolis, MN, 2001.