Preconditioned Krylov Methods Research Articles

Rates of robust superlinear convergence of preconditioned Krylov methods for elliptic FEM problems

This paper considers the iterative solution of finite element discretizations of second-order elliptic boundary value problems. Mesh independent estimations are given for the rate of superlinear convergence of preconditioned Krylov methods, involving the connection between the convergence rate and the Lebesgue exponent of the data. Numerical examples demonstrate the theoretical results.

Preconditioners for fractional diffusion equations based on the spectral symbol

AbstractIt is well known that the discretization of fractional diffusion equations with fractional derivatives , using the so‐called weighted and shifted Grünwald formula, leads to linear systems whose coefficient matrices show a Toeplitz‐like structure. More precisely, in the case of variable coefficients, the related matrix sequences belong to the so‐called generalized locally Toeplitz class. Conversely, when the given FDE has constant coefficients, using a suitable discretization, we encounter a Toeplitz structure associated to a nonnegative function , called the spectral symbol, having a unique zero at zero of real positive order between one and two. For the fast solution of such systems by preconditioned Krylov methods, several preconditioning techniques have been proposed in both the one‐ and two‐dimensional cases. In this article we propose a new preconditioner denoted by which belongs to the ‐algebra and it is based on the spectral symbol . Comparing with some of the previously proposed preconditioners, we show that although the low band structure preserving preconditioners are more effective in the one‐dimensional case, the new preconditioner performs better in the more challenging multi‐dimensional setting.

Numerical solutions to linear transfer problems of polarized radiation

Context. Numerical solutions to transfer problems of polarized radiation in solar and stellar atmospheres commonly rely on stationary iterative methods, which often perform poorly when applied to large problems. In recent times, stationary iterative methods have been replaced by state-of-the-art preconditioned Krylov iterative methods for many applications. However, a general description and a convergence analysis of Krylov methods in the polarized radiative transfer context are still lacking. Aims. We describe the practical application of preconditioned Krylov methods to linear transfer problems of polarized radiation, possibly in a matrix-free context. The main aim is to clarify the advantages and drawbacks of various Krylov accelerators with respect to stationary iterative methods and direct solution strategies. Methods. After a brief introduction to the concept of Krylov methods, we report the convergence rate and the run time of various Krylov-accelerated techniques combined with different formal solvers when applied to a 1D benchmark transfer problem of polarized radiation. In particular, we analyze the GMRES, BICGSTAB, and CGS Krylov methods, preconditioned with Jacobi, (S)SOR, or an incomplete LU factorization. Furthermore, specific numerical tests were performed to study the robustness of the various methods as the problem size grew. Results. Krylov methods accelerate the convergence, reduce the run time, and improve the robustness (with respect to the problem size) of standard stationary iterative methods. Jacobi-preconditioned Krylov methods outperform SOR-preconditioned stationary iterations in all respects. In particular, the Jacobi-GMRES method offers the best overall performance for the problem setting in use. Conclusions. Krylov methods can be more challenging to implement than stationary iterative methods. However, an algebraic formulation of the radiative transfer problem allows one to apply and study Krylov acceleration strategies with little effort. Furthermore, many available numerical libraries implement matrix-free Krylov routines, enabling an almost effortless transition to Krylov methods.

Length scale control for high-resolution three-dimensional level set–based topology optimization

This work considers an efficient approach for length scale control in high-resolution three-dimensional level set–based topology optimization. A contrast parameter–based cut element method is employed to ensure a crisp interface representation while working on fixed background hexahedral meshes. This enables the use of efficient multigrid preconditioned Krylov methods in massively parallelized computations. The minimum length scale is controlled using a projection filter without the need for β-continuation. The capabilities of the proposed approach is demonstrated on several numerical examples using more than 62 million hexahedral elements.

A short survey on preconditioners and Korovkin-type theorems

Korovkin-type theorems in the commutative as well as noncommutative set up have many applications in various fields. In this survey, we discuss the applications of noncommutative Korovkin-type theorems in numerical linear algebra. Consider the equation $$Ax=b$$ where A is an $$n \times n$$ matrix with Toeplitz structure. In 1988, Tony Chan invented the optimal circulant preconditioner C for the system and observed that they are useful in getting faster convergence for iteration processes associated with classical splittings and with (preconditioned) Krylov methods. In the same spirit, several researchers have considered this problem to get efficient preconditioners such as Hartley, $$\epsilon$$ -circulant etc. In 1999, Stefano Serra-Capizzano generalized this technique by considering an arbitrary sequence of matrix algebras $$\{M_{U_n}\}$$ that includes algebras associated with fast transforms like Fourier, Trigonometric, Hartley, Wavelet etc. He showed that in this general situation, a sequence of matrices $$\{P_{U_n}(A_n)\}$$ in $$M_{U_n}$$ is a preconditioner to the Toeplitz matrix sequence $$\{A_n\}$$ if $$P_{U_n}(A_n)-A_n$$ converges to the 0 matrix in the sense of the singular value clustering or eigenvalue clustering, if all the involved matrices are Hermitian. In 2013, these notions were generalized into the setting of operators acting on infinite dimensional Hilbert spaces. The connection of preconditioners with the noncommutative Korovkin-type theorems is the most interesting aspect of this study. The classical theorem of Korovkin unified several approximation processes such as Bernstein, Weierstrass, Fejer etc. Similarly, the Korovkin-type theorems obtained in the context of matrix/operator sequences here unify several preconditioning techniques such as Fourier, Hartley etc. These results find applications in spectral approximation problem for bounded self-adjoint operators. The most recent interests include Korovkin-type theorems for nonnormal Toeplitz operators on various function spaces. A short survey of all these developments are presented in this article.

Alternating Anderson–Richardson method: An efficient alternative to preconditioned Krylov methods for large, sparse linear systems

We present the Alternating Anderson–Richardson (AAR) method: an efficient and scalable alternative to preconditioned Krylov solvers for the solution of large, sparse linear systems on high performance computing platforms. Specifically, we generalize the recently proposed Alternating Anderson–Jacobi (AAJ) method (Pratapa et al., 2016) to include preconditioning, discuss efficient parallel implementation, and provide serial MATLAB and parallel C/C++ implementations. In serial applications to nonsymmetric systems, we find that AAR is comparably robust to GMRES, using the same preconditioning, while often outperforming it in time to solution; and find AAR to be more robust than Bi-CGSTAB for the problems considered. In parallel applications to the Helmholtz and Poisson equations, we find that AAR shows superior strong and weak scaling to GMRES, Bi-CGSTAB, and Conjugate Gradient (CG) methods, using the same preconditioning, with consistently shorter times to solution at larger processor counts. Finally, in massively parallel applications to the Poisson equation, on up to 110,592 processors, we find that AAR shows superior strong and weak scaling to CG, with shorter minimum time to solution. We thus find that AAR offers a robust and efficient alternative to current state-of-the-art solvers, with increasing advantages as the number of processors grows. Program summaryProgram Title: AARProgram Files doi:http://dx.doi.org/10.17632/669xrdcwry.1Licensing provisions: GNU General Public License 3 (GPL)Programming language: MATLAB for the MATLAB version. C/C++ for the PETSc version. C/C++ for the standalone version.External routines/libraries: MATLAB 2014a or later for the MATLAB version.PETSc 3.5.3 (http://www.mcs.anl.gov/petsc) or later andMVAPICH2 2.1 (http://mvapich.cse.ohio-state.edu/) or later for the PETSc version.MVAPICH2 2.1 or later for the standalone version.Nature of problem: Linear systems of equations.Solution method: Anderson extrapolation at periodic intervals within a preconditioned Richardson iteration.Restrictions: Jacobi preconditioning for the standalone version.

Accelerated multigrid for graph Laplacian operators

We consider multigrid type techniques for the numerical solution of large linear systems, whose coefficient matrices show the structure of (weighted) graph Laplacian operators. We combine ad hoc coarser-grid operators with iterative techniques used as smoothers. Empirical tests suggest that the most effective smoothers have to be of Krylov type with subgraph preconditioners, while the projectors, which define the coarser-grid operators, have to be designed for maintaining as much as possible the graph structure of the projected matrix at the inner levels. The main theoretical contribution of the paper is the characterization of necessary and sufficient conditions for preserving the graph structure. In this framework it is possible to explain why the classical projectors inherited from differential equations are good in the differential context and why they may behave unsatisfactorily for unstructured graphs. Furthermore, we report and discuss several numerical experiments, showing that our approach is effective even in very difficult cases where the known approaches are rather slow. As a conclusion, the main advantage of the proposed approach is the robustness, since our multigrid type technique behaves uniformly well in all cases, without requiring either the setting or the knowledge of critical parameters, as it happens when using the best known preconditioned Krylov methods.

Polynomial Preconditioned GMRES and GMRES-DR

We look at solving large nonsymmetric systems of linear equations using polynomial preconditioned Krylov methods. We give a simple way to find the polynomial. It is shown that polynomial preconditioning can significantly improve restarted GMRES for difficult problems, and the reasons for this are examined. Stability is discussed, and algorithms are given for increased stability. Next, we apply polynomial preconditioning to GMRES with deflated restarting. It is shown that this is worthwhile for sparse matrices and for problems with many small eigenvalues. Multiple right-hand sides are also considered.

Auto-tuned Krylov methods on cluster of graphics processing unit

Exascale computers are expected to have highly hierarchical architectures with nodes composed by multiple core processors (CPU; central processing unit) and accelerators (GPU; graphics processing unit). The different programming levels generate new difficult algorithm issues. In particular when solving extremely large linear systems, new programming paradigms of Krylov methods should be defined and evaluated with respect to modern state of the art of scientific methods. Iterative Krylov methods involve linear algebra operations such as dot product, norm, addition of vectors and sparse matrix–vector multiplication. These operations are computationally expensive for large size matrices. In this paper, we aim to focus on the best way to perform effectively these operations, in double precision, on GPU in order to make iterative Krylov methods more robust and therefore reduce the computing time. The performance of our algorithms is evaluated on several matrices arising from engineering problems. Numerical experiments illustrate the robustness and accuracy of our implementation compared to the existing libraries. We deal with different preconditioned Krylov methods: Conjugate Gradient for symmetric positive-definite matrices, and Generalized Conjugate Residual, Bi-Conjugate Gradient Conjugate Residual, transpose-free Quasi Minimal Residual, Stabilized BiConjugate Gradient and Stabilized BiConjugate Gradient (L) for the solution of sparse linear systems with non symmetric matrices. We consider and compare several sparse compressed formats, and propose a way to implement effectively Krylov methods on GPU and on multicore CPU. Finally, we give strategies to faster algorithms by auto-tuning the threading design, upon the problem characteristics and the hardware changes. As a conclusion, we propose and analyse hybrid sub-structuring methods that should pave the way to exascale hybrid methods.

Hierarchical Krylov and nested Krylov methods for extreme-scale computing

The solution of large, sparse linear systems is often a dominant phase of computation for simulations based on partial differential equations, which are ubiquitous in scientific and engineering applications. While preconditioned Krylov methods are widely used and offer many advantages for solving sparse linear systems that do not have highly convergent, geometric multigrid solvers or specialized fast solvers, Krylov methods encounter well-known scaling difficulties for over 10,000 processor cores because each iteration requires at least one vector inner product, which in turn requires a global synchronization that scales poorly because of internode latency. To help overcome these difficulties, we have developed hierarchical Krylov methods and nested Krylov methods in the PETSc library that reduce the number of global inner products required across the entire system (where they are expensive), though freely allow vector inner products across smaller subsets of the entire system (where they are inexpensive) or use inner iterations that do not invoke vector inner products at all. Nested Krylov methods are a generalization of inner-outer iterative methods with two or more layers. Hierarchical Krylov methods are a generalization of block Jacobi and overlapping additive Schwarz methods, where each block itself is solved by Krylov methods on smaller blocks. Conceptually, the hierarchy can continue recursively to an arbitrary number of levels of smaller and smaller blocks. As a specific case, we introduce the hierarchical FGMRES method, or h-FGMRES, and we demonstrate the impact of two-level h-FGMRES with a variable preconditioner on the PFLOTRAN subsurface flow application. We also demonstrate the impact of nested FGMRES, BiCGStab and Chebyshev methods. These hierarchical Krylov methods and nested Krylov methods significantly reduced overall PFLOTRAN simulation time on the Cray XK6 when using 10,000 through 224,000 cores through the combined effects of reduced global synchronization due to fewer global inner products and stronger inner hierarchical or nested preconditioners.

Semi-automatic sparse preconditioners for high-order finite element methods on non-uniform meshes

High-order finite elements often have a higher accuracy per degree of freedom than the classical low-order finite elements. However, in the context of implicit time-stepping methods, high-order finite elements present challenges to the construction of efficient simulations due to the high cost of inverting the denser finite element matrix. There are many cases where simulations are limited by the memory required to store the matrix and/or the algorithmic components of the linear solver. We are particularly interested in preconditioned Krylov methods for linear systems generated by discretization of elliptic partial differential equations with high-order finite elements. Using a preconditioner like Algebraic Multigrid can be costly in terms of memory due to the need to store matrix information at the various levels. We present a novel method for defining a preconditioner for systems generated by high-order finite elements that is based on a much sparser system than the original high-order finite element system. We investigate the performance for non-uniform meshes on a cube and a cubed sphere mesh, showing that the sparser preconditioner is more efficient and uses significantly less memory. Finally, we explore new methods to construct the sparse preconditioner and examine their effectiveness for non-uniform meshes. We compare results to a direct use of Algebraic Multigrid as a preconditioner and to a two-level additive Schwarz method.

Analysis of Multigrid Preconditioning for Implicit PDE Solvers for Degenerate Parabolic Equations

In this paper an implicit numerical method designed for nonlinear degenerate parabolic equations is proposed. A convergence analysis and the study of the related computational cost are provided. In fact, due to the nonlinear nature of the underlying mathematical model, the use of a fixed point scheme is required. The chosen scheme is the Newton method and its convergence is proven under mild assumptions. Every step of the Newton method implies the solution of large, locally structured, linear systems. A special effort is devoted to the spectral analysis of the relevant matrices and to the design of appropriate multigrid preconditioned Krylov methods. Numerical experiments for the validation of our analysis complement this contribution.

ADI preconditioned Krylov methods for large Lyapunov matrix equations

In the present paper, we propose preconditioned Krylov methods for solving large Lyapunov matrix equations AX + XA T + BB T = 0 . Such problems appear in control theory, model reduction, circuit simulation and others. Using the Alternating Direction Implicit (ADI) iteration method, we transform the original Lyapunov equation to an equivalent symmetric Stein equation depending on some ADI parameters. We then define the Smith and the low rank ADI preconditioners. To solve the obtained Stein matrix equation, we apply the global Arnoldi method and get low rank approximate solutions. We give some theoretical results and report numerical tests to show the effectiveness of the proposed approaches.

Preconditioned Implicit Solvers for Nonlinear PDEs in Monument Conservation

Mathematical models for the description, in a quantitative way, of the damages induced on monuments by the action of specific pollutants are often systems of nonlinear, possibly degenerate, parabolic equations. Although some of the asymptotic properties of the solutions are known, one needs a numerical approximation scheme in order to have a quantitative forecast at any time of interest. In this paper an implicit numerical method is proposed, analyzed, and numerically tested for parabolic equations of porous media type and on a systems of two PDEs that models the sulfation of marble in monuments. Due to the nonlinearity inside the differential operator of the underlying mathematical model, the use of an iterative solver for a large system of nonlinear equations is required, and every step implies the solution of large, locally structured, linear systems. A special effort is devoted to the spectral analysis of the relevant matrices and to the design of appropriate iterative solvers, with special attention to preconditioned Krylov methods and to multigrid procedures. Numerical experiments for the validation of the analysis complement this contribution.

3D phase-field simulations of interfacial dynamics in Newtonian and viscoelastic fluids

This work presents a three-dimensional finite-element algorithm, based on the phase-field model, for computing interfacial flows of Newtonian and complex fluids. A 3D adaptive meshing scheme produces fine grid covering the interface and coarse mesh in the bulk. It is key to accurate resolution of the interface at manageable computational costs. The coupled Navier–Stokes and Cahn–Hilliard equations, plus the constitutive equation for non-Newtonian fluids, are solved using second-order implicit time stepping. Within each time step, Newton iteration is used to handle the nonlinearity, and the linear algebraic system is solved by preconditioned Krylov methods. The phase-field model, with a physically diffuse interface, affords the method several advantages in computing interfacial dynamics. One is the ease in simulating topological changes such as interfacial rupture and coalescence. Another is the capability of computing contact line motion without invoking ad hoc slip conditions. As validation of the 3D numerical scheme, we have computed drop deformation in an elongational flow, relaxation of a deformed drop to the spherical shape, and drop spreading on a partially wetting substrate. The results are compared with numerical and experimental results in the literature as well as our own axisymmetric computations where appropriate. Excellent agreement is achieved provided that the 3D interface is adequately resolved by using a sufficiently thin diffuse interface and refined grid. Since our model involves several coupled partial differential equations and we use a fully implicit scheme, the matrix inversion requires a large memory. This puts a limit on the scale of problems that can be simulated in 3D, especially for viscoelastic fluids.

Robin–Robin preconditioned Krylov methods for fluid–structure interaction problems

In this work, we propose a Robin–Robin preconditioner combined with Krylov iterations for the solution of the interface system arising in fluid–structure interaction (FSI) problems. It can be seen as a partitioned FSI procedure and in this respect it generalizes the ideas introduced in [S. Badia, F. Nobile, C. Vergara, J. Comput. Phys. 227 (2008) 7027–7051]. We analyze the convergence of GMRES iterations with the Robin–Robin preconditioner on a model problem and compare its efficiency with some existing algorithms. The method is shown to be very efficient for many challenging fluid–structure interaction problems, such as those characterized by a large added-mass effect or by enclosed fluids. In particular, the possibility to solve balloon-type problems without any special treatment makes this algorithm very appealing compared to the computationally intensive existing approaches.

Adaptive numerical components for PDE‐based simulations

AbstractNumerical simulations based on nonlinear partial differential equations (PDEs) using Newton‐based methods require the solution of large, sparse linear systems of equations at each nonlinear iteration. Typically in large‐scale parallel simulations such linear systems are solved by using preconditioned Krylov methods. In many cases, especially in time‐dependent problems, the attributes of the linear systems can change throughout the stimulation, potentially leading to varying times for solving the linear systems during different nonlinear iterations. We present an approach to characterizing the nonlinear and linear system solution and using the resulting application performance information to dynamically select linear solver methods, with the goal of reducing the total time to solution. We discuss the effect of these adaptive heuristics on fluid dynamics and radiation transport codes. We also introduce general component infrastructure to support dynamic algorithm selection and adaptation in applications involving the solution of nonlinear PDEs. (© 2008 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)

A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier–Stokes equations

In recent years, considerable effort has been placed on developing efficient and robust solution algorithms for the incompressible Navier–Stokes equations based on preconditioned Krylov methods. These include physics-based methods, such as SIMPLE, and purely algebraic preconditioners based on the approximation of the Schur complement. All these techniques can be represented as approximate block factorization (ABF) type preconditioners. The goal is to decompose the application of the preconditioner into simplified sub-systems in which scalable multi-level type solvers can be applied. In this paper we develop a taxonomy of these ideas based on an adaptation of a generalized approximate factorization of the Navier–Stokes system first presented in [A. Quarteroni, F. Saleri, A. Veneziani, Factorization methods for the numerical approximation of Navier–Stokes equations, Computational Methods in Applied Mechanical Engineering 188 (2000) 505–526]. This taxonomy illuminates the similarities and differences among these preconditioners and the central role played by efficient approximation of certain Schur complement operators. We then present a parallel computational study that examines the performance of these methods and compares them to an additive Schwarz domain decomposition (DD) algorithm. Results are presented for two and three-dimensional steady state problems for enclosed domains and inflow/outflow systems on both structured and unstructured meshes. The numerical experiments are performed using MPSalsa, a stabilized finite element code.

Adaptive preconditioners for nonlinear systems of equations

The use of preconditioned Krylov methods is in many applications mandatory for computing efficiently the solution of large sparse nonlinear systems of equations. However, the available preconditioners are often sub-optimal, due to the changing nature of the linearized operator. In this work we introduce and analyse an adaptive preconditioning technique based on the Krylov subspace information generated at previous steps in the nonlinear iteration. In particular, we use an adaptive technique suggested in [J. Baglama, D. Calvetti, G.H. Golub, L. Reichel, Adaptively preconditioned GMRES algorithms, SIAM J. Sci. Comput. 20(1) (1998) 243–269] for restarted GMRES to enhance existing preconditioners with information available from previous stages in the nonlinear iteration. Numerical experiments drawn from domain decomposition techniques and fluid flow applications are used to validate the increased efficiency of our approach.

Wavelet Domain Image Reconstruction by Compactly-Supported Radial Basis Functions

In this paper we propose the use of wavelets to accelerate the solution of the System of Linear Algebraic Equations that arise from the formulation of the problem of image interpolation from scattered data by means of Compactly-Supported Radial Basis Functions. Examples demonstrate the superiority of the solution in the wavelet domain using preconditioned iterative Krylov methods.