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?

  • Do you intend this question to be a question about iterative linear solvers specifically in PETSc (which is what I would've gathered from the question body text), or to be a question about potential algorithmic failures of iterative linear solvers in a mostly software-agnostic context (which is what I would've gathered from looking at the headline alone)? – Geoff Oxberry Dec 23 '11 at 8:34
  • 4
    It has the petsc tag. The methodology is general, but I think the answer would be less useful if each "try this" did not also include the "how". Alternatively, the "how" would need to be much longer (and more error-prone for the viewer) if it needed to be explained in a software-agnostic way. If someone wants to explain how to do all these things using a different package, I will happily make the question software-agnostic and change my answer to state that it describes what to do in PETSc. Note: I added this, which is an enhanced version of an FAQ, so I could like people to this site. – Jed Brown Dec 23 '11 at 15:46
up vote 23 down vote accepted

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 petsc-users@mcs.anl.gov or petsc-maint@mcs.anl.gov) 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.

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.

  • 3
    -pc_type lu -pc_factor_mat_solver_type umfpack to use UMFPACK (or -pc_type cholesky -pc_factor_mat_solver_package cholmod for SPD problems) through PETSc, but note that UMFPACK and CHOLMOD are serial. For parallel, use -pc_factor_mat_solver_package superlu_dist or mumps, pastix, spooles. – Jed Brown Dec 23 '11 at 18:34
  • 1
    Just to be clear, the complete option set for using (e.g.) superlu_dist would be -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package superlu_dist. Is that right? – Leon Avery Aug 30 '17 at 13:27
  • I don't know. What happens if you do this? – Wolfgang Bangerth Aug 30 '17 at 23:20

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.