27

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 ...


17

I originally didn't want to give an answer because this deserves a very long treatment, and hopefully someone else will still give it. However, I can certainly give a very brief overview of the recommended approach: Perform a thorough literature search. If that fails, try every preconditioner that makes sense that you can get your hands on. MATLAB, PETSc, ...


16

Multigrid and multilevel domain decomposition methods have so much in common that each can usually be written as a special case of the other. The analysis frameworks are somewhat different, as a consequence of the different philosophies of each field. Generally speaking, multigrid methods use moderate coarsening rates and simple smoothers while domain ...


16

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 ...


15

This question is very closely related to (and possibly a duplicate of) Is variable scaling essential when solving some PDE problems numerically?. There are still good practical reasons to nondimensionalize equations, if possible: It reduces the number of independent parameters for parametric studies (which was one of the original reasons for ...


15

Have a look at Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods (Barrett et al.). You can find it here. Here's why I'm recommending this over other references: the "flowchart of iterative methods" in appendix D (last page) covers both linear solvers and preconditioners, it is short (100 pages or so), does not go into too ...


14

Others have already commented on the issue of preconditioning what I will call "monolithic" matrices, i.e. for example the discretized form of a scalar equation such as the Laplace equation, the Helmholtz equation or, if you want to generalize it, the vector-valued elasticity equation. For these things, it is clear that multigrid (either algebraic or ...


14

Interesting that this question came yesterday, since I just finished an implementation yesterday that does this. My Background Just to start of, let me know that while my education background is from scientific computing, all work I have done since graduating, including my current Ph.D. work, has been in computational electromagnetics. So, I guess our ...


13

Jack has given a good procedure for finding a preconditioner. I will try an address the question, "What makes a good preconditioner?". The operational definition is: A Preconditioner M accelerates the iterative solution of $A x = b$, and $M^{-1}$ can be applied cheaply compared to $A^{-1}$. however this does not give us any insight into designing a ...


13

Warning Solving saddle point problems involves a lot more choices than definite problems, and there are a lot more things that can go wrong. Use monitors for all levels to debug convergence, to sure that null spaces are handled correctly when auxiliary operators are singular (usually just a constant null space), and to ensure that preconditioners are ...


13

You can use additive $$ P_a^{-1} x = (B^T B)^{-1} x + (C^T C)^{-1} x, $$ multiplicative $$ P_m^{-1} x = (B^T B)^{-1} x + (C^T C)^{-1} \Big(x - A (B^T B)^{-1} x \Big), $$ or symmetric multiplicative. Methods of this class are available in PETSc using PCCOMPOSITE in PETSc. For example, petsc/src/ksp/ksp/examples/tutorials$ ./ex2 -m 100 -n 100 -...


12

To perform reasonably, polynomial preconditioners need fairly accurate spectral estimates. For ill-conditioned elliptic problems the smallest eigenvalues are usually separated such that methods like Chebyshev are far from optimal. The most interesting property of polynomial methods is that they do not require any inner products. It's actually quite popular ...


12

Though it is a relatively rare situation when you actually have to calculate an inverse of the matrix, not all techniques were created equally. I would use the term badly-conditioned instead of ill-conditioned. For badly conditioned matrices, you might opt in the SVD-route to calculate the inverse: $$ A=U\Sigma V^H \implies A^{-1}=V\Sigma^{-1}U^H. $$ If ...


11

Practical experience shows that trying to get good initial iterates has little value. For example, in the context of solving partial differential equations, if you take the solution from one mesh, interpolate it onto a finer mesh, and use that as the starting guess for something like a CG iteration to solve the same problem on the finer mesh, then it turns ...


10

It is critical to know more about the structure. It matters whether the random entries are uniformly or normally distributed and whether there is a shift or not. If there is no structure at all, then you cannot asymptotically beat a direct solve. Some comments on your proposed approaches Incomplete LU is complete LU when applied to a dense matrix. You could ...


9

I'm going to answer my own question since the following method seems to be very effective. I'm making it an answer so people can upvote or downvote it independently of the question if they think it is good or bad. Answer: use randomized matrix probing applied to the diagonal of the matrix. Let $A$ be the operator we wish to find the diagonal of, and let $\...


9

A preconditioner, say M, is an approximation on the system matrix, say A that changes the problem into another problem with improved eigenvalue spectrum. A perfect preconditioner would be inverse of A i.e inv(M) = A. Unfortunately, this inverse is normally not avaiable, too complicated to compute, requires more space to store because of the fill-in's ...


9

This is an excellent writeup but I think saying that (multilevel) DD and MG have a lot in common is not accurate, or at least not useful. The methods are very different and I don't think that expertise in one is very useful in the other. First, the two communities use different definitions of complexity: DD optimizes the condition number of the ...


9

The preconditioner for the FDM method that corresponds to the one you outline for the FEM (i.e., the Sylvester-Wathen approach) will still contain the Schur complement of the FDM matrix. The Schur complement of the FDM matrix will, in general, have the same kind of structure it has for the FEM, i.e., it will be spectrally equivalent to the identity operator ...


8

Such a nested Krylov subspace method may work quite well in practice. It may be of interest for non symmetric linear systems for which restarted GMRES stagnates and unrestarted GMRES is too expensive or uses too much memory. Some literature: GMRESR: A family of nested GMRES methods, van der Vorst, Vuik Flexible inner-outer Krylov subspace methods, Simoncini,...


7

Sure, any approximate nonlinear solver can be used as a nonlinear preconditioner. For example, if your system has directional stiffness, you could use a 1D shooting method as a preconditioner. Cai and Keyes' (2002) Nonlinearly preconditioned inexact Newton algorithms provide another useful example. In PETSc, nonlinear methods such as quasi-Newton, nonlinear ...


6

I think if there were a preconditioner which was general in the sense of working for an arbitrary dense matrix without defined structure or contents accelerated your solution rate without other drawback then you would always use it in your solution algorithim, and thus it would become part of your solver.


6

Depending on the problem structure, you can solve the ill-conditioned Augmented Lagrangian system directly. For example, BDDC/FETI-DP can solve almost-incompressible elasticity in primal form with a convergence rate independent of the Poisson ratio (piecewise constant on subdomains, but with arbitrary jumps). Similarly, multigrid methods that exactly ...


6

This is something like a truncated SVD or eigenvector expansion of your solution. If you take $$x_m = \sum_{j=1}^m \frac{q_j\cdot b}{\lambda_j}q_j$$ with $m=n$, this is the exact solution to $Ax=b$. If you take instead the $m<n$ eigenvalues in the sum instead, you get an approximation. The problem is that it's a pretty poor approximation; since your ...


6

I don't know for sure whether there is an existing name for this method, but @jessechan's suggestion of "truncated eigenvector expansion" sounds perfectly fine to me (and most people would understand it). Your question had a second part embedded, namely why is this not what everyone does to get a good initial guess for iterative solvers? The answer to this ...


6

I don't particularly care for the notation $M^{-1}$ precisely because of the confusion you find yourself in. I (and others) simply call the preconditioner $P$. The point, however, is that for iterative solvers, you generally never need a matrix or its inverse explicitly. All you need to be able to do is apply it to a vector. So, let's say you want to solve $...


5

No, Jacobi only ever corrects relative scales. It does nothing for "smooth" ill-conditioning, such as the $\kappa(A) \in O(h^{-2})$ asymptotics for second order elliptic problems. If you are using a Krylov method, the global scale is automatically corrected, but with a stationary iteration, the (constant) scaling is needed somehow (could just be in the ...


5

Meijerink and van der Vorst showed that incomplete Cholesky dose not break down for $M$-matrices. As for finite element mass matrices, you will have to specify a basis to have any hope of making that claim. Given any basis $\Phi = [\phi_0 | \phi_1 | \phi_2 | \phi_3]$ such that the mass matrix $A = \Phi^T \Phi$ is SPD, we can construct a new basis $\hat \Phi ...


5

For dense matrices without structure, polynomial preconditioning is probably the only viable method, though it has limitations (see, e.g., http://amath.colorado.edu/pub/iterative/psi-phi.ps.Z). If your matrix has entries of different magnitudes, it may be necessary that you first scale your matrix using a matching http://www.cerfacs.fr/algor/reports/1997/...


5

I have not kept up with this literature for a while, but I see these methods as not quite a genuine "another track of development". My recollection is the theory was that it invoked classic linear algebra at the end of the day and end up looking like incomplete factorization methods. That said if you have some ILU solver that you motivate with graph theory ...


Only top voted, non community-wiki answers of a minimum length are eligible