Why is my iterative linear solver not converging?

Question

What can go wrong when using preconditoned Krylov methods from KSP (PETSc's linear solver package) to solve a sparse linear system such as those obtained by discretizing and linearizing partial differential equations?

What steps can I take to determine what is going wrong for my problem?

What changes can I make to successfully and efficiently solve my linear system?

Answer 1

Initial advice

• Always run with -ksp_converged_reason -ksp_monitor_true_residual when trying to learn why a method is not converging.
• Make the problem size and number of processes as small as possible to demonstrate the failure. You often gain insight by determining what small problems exhibit the behavior that is causing your method to break down and the turn-around time is reduced. Additionally, there are some investigation techniques that can only be used for small systems.
• If the issue only arises after a large number of time steps, continuation steps, or nonlinear solve steps, consider writing the model state out when failure occurs so that you can experiment quickly.
• Alternatively, especially if your software does not have checkpoint capability, use -ksp_view_binary or MatView() to save the linear system, then use the code at $PETSC_DIR/src/ksp/ksp/examples/tutorials/ex10.c to read in the matrix and solve it (possibly with a different number of processes). This requires an assembled matrix, so it's usefulness can be somewhat limited.
• There are many possible solver choices (e.g. an infinite number available at the command line in PETSc due to an arbitrary number of levels of composition), see this question for general advice on choosing linear solvers.

Common reasons for KSP not converging

• The equations are singular by accident (e.g. forgot to impose boundary conditions). Check this for a small problem using -pc_type svd -pc_svd_monitor. Also try a direct solver with -pc_type lu (via a third-party package in parallel, e.g. -pc_type lu -pc_factor_mat_solver_package superlu_dist).
• The equations are intentionally singular (e.g. constant null space), but the Krylov method was not informed, see KSPSetNullSpace().
• The equations are intentionally singular and KSPSetNullSpace() was used, but the right hand side is not consistent. You may have to call MatNullSpaceRemove() on the right hand side before calling KSPSolve().
• The equations are indefinite so that standard preconditioners don't work. Usually you will know this from the physics, but you can check with -ksp_compute_eigenvalues -ksp_gmres_restart 1000 -pc_type none. For simple saddle point problems, try -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_detect_saddle_point. See the User's Manual and PCFIELDSPLIT man page for more details. For more difficult problems, read the literature to find robust methods and ask here (or [email protected] or [email protected]) if you want advice about how to implement them. For example, see this question for high frequency Helmholtz. For modest problem sizes, see if you can live with just using a direct solver.
• If the method converges in preconditioned residual, but not in true residual, the preconditioner is likely singular or nearly so. This is common for saddle point problems (e.g. incompressible flow) or strongly nonsymmetric operators (e.g. low-Mach hyperbolic problems with large time steps).
• The preconditioner is too weak or is unstable. See if -pc_type asm -sub_pc_type lu improves the convergence rate. If GMRES is losing too much progress in the restart, see if longer restarts help -ksp_gmres_restart 300. If a transpose is available, try -ksp_type bcgs or other methods that do not require a restart. (Note that convergence with these methods is frequently erratic.)
• The preconditioning matrix may not be close to the (possibly unassembled) operator. Try solving with a direct solver, either in serial with -pc_type lu or in parallel using a third-party package (e.g. -pc_type lu -pc_factor_mat_solver_package superlu_dist, or mumps). The method should converge in one iteration if the matrices are the same, and in a "small" number of iterations otherwise. Try -snes_type test to check the matrices if solving a nonlinear problem.
• The preconditioner is nonlinear (e.g. a nested iterative solve), try -ksp_type fgmres or -ksp_type gcr.
• You are using geometric multigrid, but some equations (often boundary conditions) are not scaled compatibly between levels. Try -pc_mg_galerkin to algebraically construct a correctly scaled coarse operator or make sure that all the equations are scaled in the same way if you want to use rediscretized coarse levels.
• The matrix is very ill-conditioned. Check the condition number using the methods described here. Try to improve it by choosing the relative scaling of components/boundary conditions. Try -ksp_diagonal_scale -ksp_diagonal_scale_fix. Perhaps change the formulation of the problem to produce more friendly algebraic equations. If you cannot correct the scaling, you may need to use a direct solver.
• The matrix is nonlinear (e.g. evaluated using finite differencing of a nonlinear function). Try different differencing parameters (e.g. -mat_mffd_type ds). Try using higher precision to make the differencing more accurate, ./configure --with-precision=__float128 --download-f2cblaslapack. Check if it converges in "easier" parameter regimes.
• A symmetric method is being used for a non-symmetric problem.
• Classical Gram-Schmidt is becoming unstable, try -ksp_gmres_modifiedgramschmidt or use a method that orthogonalizes differently, e.g. -ksp_type gcr.

Answer 2

My advice to students is to try a direct solver in these cases. The reason is that there are two classes of reasons why a solver may not converge: (i) the matrix is wrong, or (ii) there is a problem with the solver/preconditioner. Direct solvers almost always yield something that you can compare with the solution you expect, so if the answer of the direct solver looks correct, then you know the problem is with the iterative solver/precondition. On the other hand, if the answer looks wrong, the problem is with assembling the matrix and right hand side.

I typically just use UMFPACK as the direct solver. I'm sure it's simple to try something similar with PETSC.