Xavier Vasseur

Flexible GMRES with Deflated Restarting

In many situations, it has been observed that significant convergence improvements can be achieved in preconditioned Krylov subspace methods by enriching them with some spectral information. On the other hand, effective preconditioning strategies are often designed where the preconditioner varies from one step to the next so that a flexible Krylov solver is required. In this paper, we present a new numerical technique for nonsymmetric problems that combines these two features. We illustrate the numerical behavior of the new solver both on a set of small academic test examples as well as on large industrial simulation arising in wave propagation simulations.

Flexible Variants of Block Restarted GMRES Methods with Application to Geophysics

In a wide number of applications in computational science and engineering the solution of large linear systems of equations with several right-hand sides given at once is required. Direct methods based on Gaussian elimination are known to be especially appealing in that setting. Nevertheless, when the dimension of the problem is very large, preconditioned block Krylov space solvers are often considered as the method of choice. The purpose of this paper is thus to present iterative methods based on block restarted GMRES that allow variable preconditioning for the solution of linear systems with multiple right-hand sides. The use of flexible methods is especially of interest when approximate possibly iterative solvers are considered in the preconditioning phase. First we introduce a new variant of block flexible restarted GMRES that includes a strategy for detecting when a linear combination of the systems has approximately converged. This explicit block size reduction is often called deflation. We analyze the main properties of this flexible method based on deflation and notably prove that the Frobenius norm of the block residual is always nonincreasing. We also present a flexible variant based on both deflation and truncation to especially be used in case of limited memory. Finally we illustrate the numerical behavior of these flexible block methods for large industrial simulations arising in geophysics, where indefinite linear systems of size up to 1 billion unknowns with multiple right-hand sides have been successfully solved in a parallel distributed memory environment.

Constrained Regional Recovery of Continental Water Mass Time-variations from GRACE-based Geopotential Anomalies over South America

We propose a “constrained” least-squares approach to estimate regional maps of equivalent-water heights by inverting GRACE-based potential anomalies at satellite altitude. According to the energy integral method, the anomalies of difference of geopotential between the two GRACE vehicles are derived from along-track K-Band Range-Rate (KBRR) residuals that correspond mainly to the continental water storage changes, once a priori known accelerations (i.e. static field, polar movements, atmosphere and ocean masses including tides) are removed during the orbit adjustment process. Newtons first law merely enables the Difference of Potential Anomalies from accurate KBRR data and the equivalent-water heights to be recovered. Spatial constraints versus spherical distance between elementary surface tiles are introduced to stabilize the linear system to cancel the effects of the north-south striping. Unlike the “mascons” approach, no basis of orthogonal functions (e.g., spherical harmonics) is used, so that the proposed regional method does not suffer from drawbacks related to any spectrum truncation. Time series of 10-day regional maps over South America for 2006–2009 also prove to be consistent with independent data sets, namely the outputs of hydrological models, “mascons” and global GRACE solutions.

An improved two-grid preconditioner for the solution of three-dimensional Helmholtz problems in heterogeneous media‡

SUMMARY In this paper, we address the solution of three-dimensional heterogeneous Helmholtz problems discretized with second-order finite difference methods with application to acoustic waveform inversion in geophysics. In this setting, the numerical simulation of wave propagation phenomena requires the approximate solution of possibly very large indefinite linear systems of equations. For that purpose, we propose and analyse an iterative two-grid method acting on the Helmholtz operator where the coarse grid problem is solved inaccurately. A cycle of a multigrid method applied to a complex shifted Laplacian operator is used as a preconditioner for the approximate solution of this coarse problem. A single cycle of the new method is then used as a variable preconditioner of a flexible Krylov subspace method. We analyse the properties of the resulting preconditioned operator by Fourier analysis. Numerical results demonstrate the effectiveness of the algorithm on three-dimensional applications. The proposed numerical method allows us to solve three-dimensional wave propagation problems even at high frequencies on a reasonable number of cores of a distributed memory computer. Copyright

A Modified Block Flexible GMRES Method with Deflation at Each Iteration for the Solution of Non-Hermitian Linear Systems with Multiple Right-Hand Sides

We propose a variant of the block GMRES method for the solution of linear systems of equations with multiple right-hand sides. We investigate a deflation strategy to detect when a linear combination of approximate solutions is already known that avoids performing expensive computational operations with the system matrix. This is especially useful when the cost of the preconditioner is supposed to be larger than the cost of orthogonalization in the block Arnoldi procedure. We specifically focus on the block GMRES method incorporating deflation at the end of each iteration proposed by Robbe and Sadkane [M. Robbe and M. Sadkane, Linear Algebra Appl., 419 (2006), pp. 265--285]. We extend their contribution by proposing that deflation be performed also at the beginning of each cycle. This change leads to a modified least-squares problem to be solved at each iteration and gives rise to a different behavior especially when multiple restarts are required to reach convergence. Additionally we investigate truncatio...

Two-Level preconditioned Krylov subspace methods for the solution of three-dimensional heterogeneous Helmholtz problems in seismics

In this paper we address the solution of three-dimensional heterogeneous Helmholtz problems discretized with compact fourth-order finite difference methods with application to acoustic waveform inversion in geophysics. In this setting, the numerical simulation of wave propagation phenomena requires the approximate solution of possibly very large linear systems of equations. We propose an iterative two-grid method where the coarse grid problem is solved inexactly. A single cycle of this method is used as a variable preconditioner for a flexible Krylov subspace method. Numerical results demonstrate the usefulness of the algorithm on a realistic three-dimensional application. The proposed numerical method allows us to solve wave propagation problems with single or multiple sources even at high frequencies on a reasonable number of cores of a distributed memory cluster.

A Flexible Generalized Conjugate Residual Method with Inner Orthogonalization and Deflated Restarting

This work is concerned with the development and study of a minimum residual norm subspace method based on the generalized conjugate residual method with inner orthogonalization (GCRO) method that allows flexible preconditioning and deflated restarting for the solution of nonsymmetric or non-Hermitian linear systems. First we recall the main features of flexible generalized minimum residual with deflated restarting (FGMRES-DR), a recently proposed algorithm of the same family but based on the GMRES method. Next we introduce the new inner-outer subspace method named FGCRO-DR. A theoretical comparison of both algorithms is then made in the case of flexible preconditioning. It is proved that FGCRO-DR and FGMRES-DR are algebraically equivalent if a collinearity condition is satisfied. While being nearly as expensive as FGMRES-DR in terms of computational operations per cycle, FGCRO-DR offers the additional advantage to be suitable for the solution of sequences of slowly changing linear systems (where both the matrix and right-hand side can change) through subspace recycling. Numerical experiments on the solution of multidimensional elliptic partial differential equations show the efficiency of FGCRO-DR when solving sequences of linear systems.

Multigrid based preconditioners for the numerical solution of two-dimensional heterogeneous problems in geophysics

We study methods for the numerical solution of the Helmholtz equation for two-dimensional applications in geophysics. The common framework of the iterative methods in our study is a combination of an inner iteration with a geometric multigrid method used as a preconditioner and an outer iteration with a Krylov subspace method. The preconditioning system is based on either a pure or shifted Helmholtz operator. A multigrid iteration is used to approximate the inverse of this operator. The proposed solution methods are evaluated on a complex benchmark in geophysics involving highly variable coefficients and high wavenumbers. We compare this preconditioned iterative method with a direct method and a hybrid method that combines our iterative approach with a direct method on a reduced problem. We see that the hybrid method outperforms both the iterative and the direct approach.

Low rank updates in preconditioning the saddle point systems arising from data assimilation problems

The numerical solution of saddle point systems has received a lot of attention over the past few years in a wide variety of applications such as constrained optimization, computational fluid dynamics and optimal control, to name a few. In this paper, we focus on the saddle point formulation of a large-scale variational data assimilation problem, where the computations involving the constraint blocks are supposed to be much more expensive than those related to the (1, 1) block of the saddle point matrix. New low-rank limited memory preconditioners exploiting the particular structure of the problem are proposed and analysed theoretically. Numerical experiments performed within the Object-Oriented Prediction System are presented to highlight the relevance of the proposed preconditioners.

Limited memory preconditioners for symmetric indefinite problems with application to structural mechanics

Summary This paper presents a class of limited memory preconditioners (LMP) for solving linear systems of equations with symmetric indefinite matrices and multiple right-hand sides. These preconditioners based on limited memory quasi-Newton formulas require a small number k of linearly independent vectors and may be used to improve an existing first-level preconditioner. The contributions of the paper are threefold. First, we derive a formula to characterize the spectrum of the preconditioned operator. A spectral analysis of the preconditioned matrix shows that the eigenvalues are all real and that the LMP class is able to cluster at least k eigenvalues at 1. Secondly, we show that the eigenvalues of the preconditioned matrix enjoy interlacing properties with respect to the eigenvalues of the original matrix provided that the k linearly independent vectors have been prior projected onto the invariant subspaces associated with the eigenvalues of the original matrix in the open right and left half-plane, respectively. Third, we focus on theoretical properties of the Ritz-LMP variant, where Ritz information is used to determine the k vectors. Finally, we illustrate the numerical behaviour of the Ritz limited memory preconditioners on realistic applications in structural mechanics that require the solution of sequences of large-scale symmetric saddle-point systems. Numerical experiments show the relevance of the proposed preconditioner leading to a significant decrease in terms of computational operations when solving such sequences of linear systems. A saving of up to 43% in terms of computational effort is obtained on one of these applications. Copyright