This is a difficult question to answer objectively for the reasons outlined in Matt’s talk that you’ve cited. I’ll write my opinion here based on my past work as a consumer of computational linear algebra research and not a contributor.
I’ll consider two specific methods here: stationary (direct factorisation) and iterative (preconditioned Krylov subspace search), each in the context of the d- dimensional finite element discretisation of the Laplace operator defined on a mesh yielding n degrees of freedom:
-
Direct solution (e.g. mumps, superlu_dist):
- Complexity \mathcal{O}(n^{\frac{3}{2}}) in 2D and \mathcal{O}(n^{2}) in 3D.
- Robust. Unless your system is singular, you will arrive at the solution.
- Precise. The solution of the linear system is computed to machine pecision.
- Fast (for small problems). The growth constant is typically smaller than iterative solvers. If you have a small problem, you will likely solve that problem faster with a direct solver than an iterative solver.
- “Black box”. Very few parameters need to be tuned to arrive at a solution.
-
Algebraic multigrid preconditioned Krylov iterative method:
- Complexity \mathcal{O}(n).
- The correct algorithm to employ depends on your discretised operator: Is it hermitian? Is it positive-definite / semi-definite / indefinite?
- Designing a preconditioner can be difficult. Are there orders of magnitude difference in the problems’ material coefficients? How does one capture geometric effects arising in algebraic systems? How best to coarsen and refine the problem for a hierarchy of grids?
- Convergence rate depends on mesh quality.
- Many parameters to consider: rate of coarsening, paritioning, which solvers to employ on which levels, the cycle pattern to be employed…
- The computed solution is an approximation of the true linear algebra solution. You are introducing a new numerical error, in addition to finite element approximation error, by imprecisely solving the linear system. Typically more computational effort is required for a more precise approximation of the linear system’s solution.
With the pros and cons of the above problems in place, what should we use in modern finite element analysis? Clearly direct solvers have earned their popularity. They are fast and existing packages are easy to use. But there’s a problem…
Although a solver may perform exceptionally well for a problem with n \backsim 10^5 degrees of freedom, what if our scientific application demands the fidelity only offered by a problem with n \geq 10^9 degrees of freedom? For these large problems, even if the growth constant is small, the growth rates of \mathcal{O}(n^\frac{3}{2}) and \mathcal{O}(n^2) cannot be considered merely a disadvantage. They ensure that computing a solution of a large system with a direct method is impossible in reasonable time.
Therefore I’m not solely concerned with marginal performance improvements offered by spending hours tweaking iterative/direct solver schemes by trial and error. What’s crucially important in modern computational simulations is scaling. One must be very careful that the relationship between problem size and time-to-solution is optimal (i.e. linear). We therefore rely on iterative methods, and are humbled by their development by very talented numericists in the literature.