Desire Nuentsa
email: de-sire.nuentsa_wakam@inria.fr

Jocelyne Erhel
email: joce-lyne.erhel@inria.fr

PARALLELISM AND ROBUSTNESS IN GMRES WITH THE NEWTON BASIS AND THE DEFLATED RESTARTING

Keywords: Augmented Krylov subspaces, Adaptive Deated GMRES, Newton basis, Hybrid linear solvers AMS subject classications. 65F10, 65F15, 65F22

The GMRES iterative method is widely used as Krylov subspace technique for solving sparse linear systems when the coecient matrix is nonsymmetric and indenite. The Newton basis implementation has been proposed on distributed memory computers as an alternative to the classical approach with the Arnoldi process. The aim of our work here is to introduce a modication based on deation techniques. This approach builds an augmented subspace in an adaptive way to accelerate the convergence of the restarted formulation. In our numerical experiments, we show the benets of using this implementation with hybrid direct/iterative methods to solve large linear systems.

1. Introduction. In this paper, we are interested in solving large systems of linear algebraic equations Ax = b, (1.1) where A is a n × n real nonsingular matrix, b and x are n-dimensional real vectors. Practical algorithms transform the original problem (1.1) to the following

M -1 L AM -1 R x = M -1 L b, x = M R x (1.2)
where M -1 L and M -1 R are the action of preconditioning the system at left (M R = I), at right (M L = I) or both. On parallel computers, we assume that these preconditioners are formulated from some algebraic decomposition of the input matrix. However, they can be any approximation of the inverse of the matrix A and we refer the reader to the survey on the preconditioning techniques [START_REF] Benzi | Preconditioning techniques for large linear systems: a survey[END_REF]. These preconditioners are generally combined with Krylov subspace methods as accelerators. The GMRES method [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] is widely used in this context. From this method, many improvements have been proposed to enhance its robustness and parallel eciency; see for instance [START_REF] Bai | A Newton basis GMRES implementation[END_REF][START_REF] Burrage | On the performance of various adaptive preconditioned GMRES strategies[END_REF][START_REF] Chapman | Deated and augmented Krylov subspace techniques[END_REF][START_REF] Erhel | A parallel GMRES version for general sparse matrices[END_REF][START_REF] Erhel | Restarted GMRES preconditioned by deation[END_REF][START_REF]A parallel preconditioned GMRES algorithm for sparse matrices[END_REF][START_REF] Morgan | A restarted GMRES method augmented with eigenvectors[END_REF][START_REF] Sidje | Alternatives for parallel Krylov subspace basis computation[END_REF][START_REF] Hoemmen | Communication-avoiding Krylov subspace methods[END_REF][START_REF] Philippe | On the generation of Krylov subspace bases[END_REF]. In this work, we propose a new formulation of the method which combines two main approaches, namely the Newton basis GMRES [START_REF] Bai | A Newton basis GMRES implementation[END_REF] and the augmented basis for the restarted GMRES [START_REF] Morgan | A restarted GMRES method augmented with eigenvectors[END_REF]. Our approach benets from the enhanced parallelism in the former and the robustness in the latter. For the sake of clarity, we give here the formulation of the GMRES algorithm as rst proposed by Saad and Schultz [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF].

We consider in this paper the right preconditioned matrix B ≡ AM -1 . The proposed algorithms can be derived with less eort for the left preconditioned matrix. Given an initial guess x 0 , the GMRES method nds the j-th approximate solution x j of the form x j ∈ x 0 + M -1 K j (B, r 0 ), (1.3) where r 0 = b -Ax 0 is the initial residual vector and K j (B, r 0 ) is the j-th Krylov subspace dened as K j (B, r 0 ) = span{r 0 , Br 0 , . . . , B j-1 r 0 }. (1.4) The goal behind GMRES is to minimize at each step the Euclidian norm of the residual, i.e b -Ax j = min (1.5)

An orthonormal basis V j+1 = [v 0 , . . . , v j ] 1 of K j+1 (B, r 0 ) is generated such that

v 0 = r 0 /β, β = r 0 , BV j = V j+1 H j+1,j = V j H j + h {j+1,j} v j e T j , (1.6)
It is therefore proved [START_REF] Saad | Iterative methods for sparse linear systems[END_REF] that (1.5) reduces to βe 1 -H j+1,j y j = min y∈R j βe 1 -H j+1,j y

(1.7)

and the approximate solution x j can be written as

x j = x 0 + M -1 V j y j .
(1.8)

Our work combines two improvements of this method. In GMRES(m), the method restarts at some step m to save the storage and the computational requirements as the iterations proceed. The deated and augmented approaches [START_REF] Baglama | Adaptively preconditioned GM-RES algorithms[END_REF][START_REF] Burrage | On the performance of various adaptive preconditioned GMRES strategies[END_REF][START_REF] Erhel | Restarted GMRES preconditioned by deation[END_REF][START_REF] Kharchenko | Eigenvalue translation based preconditioners for the GMRES(k) method[END_REF][START_REF] Morgan | A restarted GMRES method augmented with eigenvectors[END_REF] keep some useful information at the time of the restart to enhance robustness. We briey review these methods in section 2. The second improvement builds the orthonormal basis with a parallel algorithm that reduces the number of exchanged MPI messages on distributed-memory computers. Indeed, the time to exchange a MPI message of size n words is given by γ + δ • n where γ is the network latency (in seconds) and δ is the inverse bandwith (in seconds per word). Hence, reducing the number of MPI messages is equivalent to minimize the overhead due to this latency. From the original formulation of GMRES, the basis V is built and orthogonalized by the modied Gram-Schmidt implementation (MGS) of the Arnoldi process (V is referred to as Arnoldi basis). This process induces a high communication overhead due to the numerous inner products : for instance, a GMRES cycle of m iterations requires approximately 1 2 (m 2 + 3m) global communications for the inner products. On high latency networks, the start-up time due to these collective communication can easily dominate. Moreover, the kernel operations in MGS have a very low granularity which does not fully benet from the computer architecture. In the classical Gram-Schmidt implementation (CGS) of the Arnoldi process, the communication time can be reduced by accumulating and broadcasting multiple inner products together. However the low granularity of the kernel operations in the orthogonalization procedure remains because of the sequential form of the Arnoldi process. Moreover, CGS is more sensitive to rounding errors than MGS [START_REF] Jalby | Stability analysis and improvement of the block Gram-Schmidt algorithm[END_REF]. Alternative implementations [START_REF] Bai | A Newton basis GMRES implementation[END_REF][START_REF] Sturler | Iterative Methods on Distributed Memory Computers[END_REF][START_REF] Erhel | A parallel GMRES version for general sparse matrices[END_REF][START_REF]A parallel preconditioned GMRES algorithm for sparse matrices[END_REF][START_REF] Hoemmen | Communication-avoiding Krylov subspace methods[END_REF][START_REF] Joubert | Parallelizable restarted iterative methods for nonsymmetric linear systems. part II: Parallel implementation[END_REF][START_REF] Sidje | Alternatives for parallel Krylov subspace basis computation[END_REF][START_REF] Philippe | On the generation of Krylov subspace bases[END_REF] have been proposed. They divide the process into two main phases : rst, a nonorthogonal basis of the Krylov subspace is generated and orthogonalized as a group of vectors in the second phase. As rst proposed by Bai, Hu and Reichel [START_REF] Bai | A Newton basis GMRES implementation[END_REF], the a priori basis is built with the aid of Newton polynomials.

We will refer to this as the Newton-basis GMRES. Later on, the orthogonalization is done by replacing the vector-vector operations of the MGS method by the task of computing a QR factorization of a dense method. De Sturler [START_REF] Sturler | Iterative Methods on Distributed Memory Computers[END_REF] analyzes the parallel implementation of the second phase and suggests a distributed MGS factorization to overlap communication with computation. Sidje and Philippe [START_REF] Sidje | parallel krylov subspace basis computation[END_REF], Erhel [START_REF] Erhel | A parallel GMRES version for general sparse matrices[END_REF] and Sidje [START_REF] Sidje | Alternatives for parallel Krylov subspace basis computation[END_REF] use a dierent orthogonalization strategy called RODDEC which combines the Householder factorization with Givens rotations and requires only point-to-point communication. Demmel et al [START_REF] Demmel | Communication-optimal parallel and sequential QR and LU factorizations[END_REF] propose a dierent QR factorization called Tall Skinny QR (TSQR) which reorganizes the computation to reduce the memory access and exploit the data locality.

Our proposal in this work is to combine the Newton-basis GMRES with the augmented and deated GMRES. The new approach is simple and can be used together with any of the previous orthogonalization strategies once the augmented a priori basis is built. The motivation of our work is two-fold: previous studies [START_REF] Philippe | On the generation of Krylov subspace bases[END_REF] have shown that when the size of the Newton basis grows, the vectors become increasingly dependent. As a result, the method may experience a slow convergence rate. With the new approach, the basis is kept small and augmented with some useful approximate eigenvectors. The second motivation is related to GMRES(m) preconditioned by domain decomposition methods. Indeed, with Schwarz-based preconditioners, when the number of subdomains increases, the preconditioner becomes less and less robust and the method requires more iterations to converge. In this situation, the basis size is usually increased to prevent the stagnation. In the proposed approach, we show that by adding adaptively more approximate eigenvectors, the convergence rate is improved. Recently, Mohiyuddin et al [START_REF] Mohiyuddin | Minimizing communication in sparse matrix solvers[END_REF] and Hoemmen [START_REF] Hoemmen | Communication-avoiding Krylov subspace methods[END_REF] propose a new formulation in their communication avoiding GMRES which does not require the Krylov basis size to be equal to the number of vectors generated a priori. Their formulation builds the Krylov basis with several steps of the Arnoldi process where each step builds a set of vectors with the Newton polynomials. Our proposed approach can be used as well with their formulation for the problems that are very sensitive to the restarting procedure in GMRES. We show indeed in section 3 how the augmented basis can be formulated in their proposed approach. The remaining part of this paper is organized as follows : in section 2, we review briey how the deation of eigenvalues is used in the restarted GMRES. In section 3, we derive the new approach combining deation to the Newton basis GMRES and we discuss on the parallel implementation. section 4 is focused on numerical experiments to show the benets of the proposed approach.

2. Restarted GMRES accelerated by deation. A practical implementation of GMRES is based on restarting a minimum residual iteration when the correction space reaches a given dimension m. At the time of restart, information from the previous Krylov subspace is discarded and the orthogonality between successive Krylov subspaces is not preserved. The worst case is when the successive generated Krylov subspaces are very close. As a result, there is no signicant reduction in the residual norm and the iterative process stagnates. Deation techniques are a class of acceleration strategies that collects useful information at the time of restart mainly to avoid this stagnation and improve the convergence rate. The main idea behind these methods is to remove the smallest eigencomponents from the residual vector as they are known to slow down the convergence of GMRES [START_REF] Van Der Vorst | The superlinear convergence behaviour of GMRES[END_REF]. For a general analysis of acceleration strategies in the minimal residual methods, we refer the reader to Eiermann, Ernst, and Schneider [START_REF] Eiermann | Analysis of acceleration strategies for restarted minimal residual methods[END_REF]. For the general Krylov subspace methods, the recent reviews in [START_REF] Erhel | Some properties of Krylov projection methods for large linear systems[END_REF][START_REF] Simoncini | Recent computational developments in Krylov subspace methods for linear systems[END_REF] are also of great interest.

In deation techniques, the Krylov subspaces are enriched by some approximation of invariant subspaces associated to a selected group of eigenvalues (generally the smallest ones). Two strategies are often used, namely by preconditioning the linear system [START_REF] Burrage | On the performance of various adaptive preconditioned GMRES strategies[END_REF][START_REF] Erhel | A parallel GMRES version for general sparse matrices[END_REF][START_REF] Kharchenko | Eigenvalue translation based preconditioners for the GMRES(k) method[END_REF] or by augmenting the Krylov subspace [START_REF] Morgan | A restarted GMRES method augmented with eigenvectors[END_REF][START_REF]GMRES with deated restarting[END_REF].

Let λ 1 , λ 2 , . . . , λ n be the eigenvalues of B. Consider an exact invariant subspace of dimension r corresponding to the r smallest eigenvalues λ 1 , λ 2 , . . . , λ r and an orthonormal basis of this subspace, written U in the matrix form. Using the preconditioner M -1 dened as:

M -1 = I n + U (|λ n |T -1 -I r )U T , , T = U T BU, (2.1)
the eigenvalues of B M -1 are λ r+1 , . . . , λ n , |λ n | with a multiplicity at least r [START_REF] Erhel | Restarted GMRES preconditioned by deation[END_REF]. Then, GMRES(m) applied to B M -1 converges faster since the smallest eigencomponents that slow down the convergence are deated. In practice, U is replaced by a basis of an approximate invariant subspace.

The actual implementation in [START_REF] Erhel | Restarted GMRES preconditioned by deation[END_REF] and the improvements in [START_REF] Burrage | On the performance of various adaptive preconditioned GMRES strategies[END_REF] rely on the approximation of U which is updated at each restart by computing the Ritz values from the Arnoldi relation in equation (1.6) to yield a more accurate basis. In [START_REF] Baglama | Adaptively preconditioned GM-RES algorithms[END_REF], U is computed by the Implicit-Restarted Arnoldi (IRA) process and the result is used to form a left preconditioner. The adaptive preconditioner by Kharchenko and Yeremin [START_REF] Kharchenko | Eigenvalue translation based preconditioners for the GMRES(k) method[END_REF] is built such that the Ritz values which approximate the largest eigenvalues of B are translated to a cluster around one.

Complexity overhead is similar in the preconditioning and the augmented subspace approaches. In both cases, computing a new basis at each restart involves a small eigenvalue problem. Also, it is easy to update the preconditioning matrix with a new basis U [START_REF] Burrage | On the performance of various adaptive preconditioned GMRES strategies[END_REF]. However, applying the preconditioner M -1 in equation (2.1) may induce overhead with a large number of processors, because of global MPI reduction operations. On the other hand, in the augmented approach presented below [START_REF] Morgan | A restarted GMRES method augmented with eigenvectors[END_REF][START_REF]GMRES with deated restarting[END_REF], the communication is reduced since there is no global reduction operation.

The augmented approaches form the new approximation with a projection onto a subspace C = K m (B, r 0 ) + W, where W = span{u 0 , . . . , u r-1 }. Minimal changes are required to the existing kernel operations as the vectors are directly added to the existing Krylov basis. Moreover, when the vectors u 0 , . . . , u r-1 are the harmonic Ritz vectors, Morgan [START_REF]GMRES with deated restarting[END_REF] shows that the augmented subspace C is itself a Krylov subspace and writes

C = K m (B, r 0 ) + W = K m+r (B, q m (B)r 0 ), (2.2)
where q m (B) is a polynomial of degree at most m.

3. Deated GMRES in the Newton basis. In this section, We derive the new implementation of the GMRES algorithm where the Krylov subspace is spanned by the Newton polynomials and augmented with eigenvectors. Compared to the previous Newton basis implementations [START_REF] Bai | A Newton basis GMRES implementation[END_REF][START_REF] Erhel | A parallel GMRES version for general sparse matrices[END_REF][START_REF] Sidje | Alternatives for parallel Krylov subspace basis computation[END_REF][START_REF] Hoemmen | Communication-avoiding Krylov subspace methods[END_REF], the main dierence is that the new approach uses the deation strategies to recover the information that is lost at the time of restart. Hence for the problems that are sensitive to the restarting procedure, our implementation should converge faster than the previous approaches for the same basis size. Compared to the GMRES-E by Morgan [START_REF] Morgan | A restarted GMRES method augmented with eigenvectors[END_REF], our approach communicates less and should produce kernels that are better suited for parallel computations. Compared to GMRESDR of Morgan [START_REF]GMRES with deated restarting[END_REF] however, we have not investigated whether the proposed augmented basis is itself a Krylov basis and we left it as future work.

3.1. Augmenting the Newton basis. We now derive the proposed approach.

Our motivation is to get a Arnoldi-like relation for the augmented basis when the eigenvectors are added at the end of the Newton basis. Let B be the preconditioned matrix, x 0 an initial guess and r 0 the initial residual vector. A m-dimensional Krylov subspace is spanned by the Newton polynomials applied to r 0 of the form P j+1 (B) := σ j+1 (B -λ j+1 I)P j (B), j = 0, 1, 2, . . . , (3.1) where P 0 (B) := r 0 , σ j and λ j ∈ R (See [START_REF] Bai | A Newton basis GMRES implementation[END_REF][START_REF] Philippe | On the generation of Krylov subspace bases[END_REF]). We discuss later on the choice of these scalars.

Let W = span{u 0 , u 1 , . . . , u r-1 } be a coarse subspace of dimension r. We discuss about the computation of these vectors in section 3.2.5. Let s = m + r and C s the augmented subspace dened by

C s = K m (B, r 0 ) + W (3.2)
as previously. We write the basis of W in the matrix form as

U r = u 0 . . . u r-1 .
We also dene the augmented subspace C s+1 by

Ĉs+1 = K m+1 (B, r 0 ) + BW (3.3) Proposition 3.1. Suppose that dim(C s+1 ) = s + 1, there exists an orthonormal basis V s+1 = V m+1 Vr ∈ R n×(s+1) of the subspace Ĉs+1 such that V m+1 is an orthonormal basis of K m+1 (B, r 0 ) and W s = V m U r (3.4)
is a basis of C s . These basis are related by an Arnoldi-like relation written

BW s = V s+1 Hs = V s H s + h s+1,s v s+1 e T s , (3.5)
where Hs ∈ R (s+1)×s is an upper Hessenberg matrix. Moreover, the vector x s ∈ R n given by

x s = x 0 + M -1 W s y s , (3.6)
where y s ∈ R n solves the least-square problem J s (y) dened by

J s (y) = βe 1 -Hs y 2 , β = r 0 2 (3.7) minimizes the residual norm b -Ax s over x 0 + M -1 C s .
Proof. From k 0 = r 0 / r 0 2 , a set of vectors k j can be generated such that

σ j+1 k j+1 =    (B -λ j+1 I)k j if 0 ≤ j ≤ m -1 Bu j-m if m ≤ j ≤ s -1.
(3.8)

where λ j and σ j , (j = 1, 2, . . .) are user-specied real scalars. We discuss in sections 3.2.1 and 3.2.2 on their optimal choice. In matrix form, the relation (3.8) writes

BK m = K m+1 Tm BU r = Kr D r (3.9) with K m+1 = k 0 k 1 . . . k m , Kr = k m+1 . . . k s , Tm =            λ 1 σ 1 λ 2 σ 2 λ 3 . . . . . . λ m-1 σ m-1 λ m σ m            ∈ R (m+1)×m , (3.10)
and D r = diag{σ m+1 , . . . , σ s } ∈ R r×r . A QR factorization of K m+1 Kr yields

K m+1 K r = V s+1 R s+1 = V m+1 Vr R m+1 R m+1,r 0 R r .
(3.11)

Thus we get BV m = BK m R -1 m = V m+1 R m+1 Tm R -1 m and BU r = Kr D r = (V m+1 R m+1,r + Vr R r )D r .
Using equations (3.9) and (3.11), we get

BW s = V s+1 Hs , (3.12)
where

Hs = R s+1 Tm 0 0 D r R -1 m 0 0 I r = H s h s+1,s e T s . (3.13)
The rst part is thus proved (equation 3.5). The second part of the proposition is similar to the optimality property in the augmented GMRES [START_REF] Chapman | Deated and augmented Krylov subspace techniques[END_REF]Algorithm 2.1]. Consider an arbitrary vector x s = x 0 + M -1 W s y in the ane space x 0 + M -1 C s , the corresponding residual vector can be expressed as:

b -Ax s = b -A(x 0 + M -1 W s y) (3.14) = r 0 -V s+1 Hs y = βk 0 -V s+1 Hs y = V s+1 (βe 1 -Hs y) (3.15)
where β = r 0 and e 1 = [1, 0, . . . , 0] T . If we denote by J s (y) the function

J s (y) = b -A[x 0 + M -1 W s y] 2 ,
it comes from equation 3.15 and the fact that V s+1 is orthogonal that

J s (y) = βe 1 -Hs y 2 .
(3.16) Thus by taking the vector y s ∈ R m which minimizes the function J s (y), the approximate solution

x s = x 0 + M -1 W s y s will have the smallest residual in x 0 + M -1 C s .
We refer to the matrix W s as the augmented Newton basis of the subspace C s , and to the matrix V s+1 as the augmented orthonormal basis of the subspace Ĉs+1 .

The induced GMRES is the augmented Newton basis GMRES, which we denote by AGMRES.

It is also possible to choose K m U r as a basis of C s , avoiding the computation of R -1 m . However, the vectors U r are computed at each restart by using a Rayleigh-Ritz procedure (see section 3.2.5), which requires the storage of V s+1 and the computation of V T s+1 W s . Thus, it saves memory and computation to choose

W s = V m U r .
This proof assumes that the basis vectors k j , j = 0, . . . , s are generated through one pass in the kernel computation of equation 3.8. There are some situations where m is too large to guarantee a good robustness (ill-conditioned basis) or good performance (best value for data locality in multicore nodes). In a recent work, Hoemmen [START_REF] Hoemmen | Communication-avoiding Krylov subspace methods[END_REF] uses the µ-step Arnoldi of Kim and Chronopoulos [START_REF] Kim | An ecient parallel algorithm for extreme eigenvalues of sparse nonsymmetric matrices[END_REF] 2 in his Arnoldi(µ, t) to build the basis vectors through multiple passes of the kernel computation in equation 3.8. Hence the process of computing s basis vectors is divided into t steps where each step generates µ basis vectors with the Newton polynomials. The restart length is thus s = µ • t. We show in the following that the proposition 3.1 holds in the case of a µ-step basis. We explain the basis idea with t = 2.

Let k 0 be the starting vector and s = µ • t + r = m + r. As from the rst part of equation 3.8 we generate the sequence of vectors

σ j+1 k j+1 = (B -λ j+1 I)k j , 0 ≤ j ≤ µ -1.
(3.17)

It comes that

BK (0) µ = K (0) µ+1 T (0) µ (3.18)
where

K (0) µ+1 = [k 0 , k 1 , . . . , k µ ] ∈ R n×(µ+1) , T (0) µ ∈ R (µ+1)×µ is a bidiagonal matrix. A QR factorization of K (0) µ+1 gives K (0) µ+1 = V (0) µ+1 R (0) µ+1 (3.19)
and thus

BK (0) µ = V (0) µ+1 H(0) µ = V (0) µ H (0) µ + h µ+1,µ V (0) µ+1 e µ e T µ , (3.20)
where

H(0) µ = R (0) µ+1 T (0)
µ and e µ is the µ th unit vector. This rst step is just the derivation of the Arnoldi-like relation for the (non-augmented) Newton basis. Note that we don't have a mathematically equivalent Arnoldi relation as in equation 1.6.

H(0)

µ+1 is not equal in exact arithmetic to the Hessenberg matrix H of that equation as we avoid dealing with the term (R (0) µ+1 ) -1 . However, the columns of V (0) µ+1 form an orthogonal basis of a µ-dimensional Krylov basis. From its last column as a starting vector, we can thus build the second µ-step basis. At this step, we add the eigenvectors in the subspace by augmenting the µ-step basis.

Let k µ = V (0) µ+1 e µ , a µ-step augmented basis is generated as follows:

σ j+1 k j+1 =    (B -λ µ-j+1 I)k j if µ ≤ j ≤ m -1 Bu j-m if m ≤ j ≤ s -1. (3.21)
which is a Hessenberg matrix. The rst part of the proposition 3.1 is thus proved and the second part is similar to the previous proof. GMRES with the µ-step Newton basis is useful to control the conditioning of the basis generated with the Newton polynomials by choosing a suitable value of µ. On multicore nodes, a well chosen value of µ will also improve the data locality during the computation of the kernel computations (Generation of the basis and orthogonalization) [START_REF] Demmel | Communication-optimal parallel and sequential QR and LU factorizations[END_REF][START_REF] Hoemmen | Communication-avoiding Krylov subspace methods[END_REF]. The drawback with this formulation is that when the new set of basis vectors is orthogonalized against all the previous vectors already computed, it is important to perform a good QR factorization update. Sometimes when a block Gram-Schmidt process is used, a reorthogonalization strategy should be performed to avoid loss of orthogonality, see for instance [START_REF] Jalby | Stability analysis and improvement of the block Gram-Schmidt algorithm[END_REF]. This process induces more computational cost as the number t of steps increases. As for the scalar formulation, the augmented basis will thus help to reduce this cost by reducing the number of steps t. We do not further investigate reorthogonalization and we focus in this paper to the basic implementation of a (µ + r)-step augmented Newton basis.

3.2. AGMRES : Augmented Newton-basis GMRES. This section discusses the parallel implementation of the GMRES method where the Newton basis is augmented with a few selection of approximate eigenvectors. The main steps are outlined in Algorithm 3.2.

If we compare AGMRES with the previous related implementations of the GM-RES method, we can make the following observations :

• Compared to the standard GMRES method, AGMRES produces kernel that are better suited for parallel computations during the generation of the orthogonal Krylov basis in steps 6 and 8. However, in addition to the basis V m+1 , it keeps 2r vectors U r and Vr . It is thus worth mentioning that AGMRES(m, r) requires as much memory as GMRES(m + 2r). Nevertheless, we show in the provided numerical experiments that, on most test cases, AGMRES(m, r) produces a better convergence acceleration than GMRES(m+ 2r) and even GMRES(2m). • The GMRES-E of Morgan [START_REF] Morgan | A restarted GMRES method augmented with eigenvectors[END_REF] keeps a second basis as well. However its implementation is based on the Arnoldi process. It will thus communicate more for the same convergence behavior. Our implementation includes an adaptive strategy that will allow to increase the number of extracted eigvenvectors if necessary. • Compared to CA-GMRES of Hoemmen [START_REF] Hoemmen | Communication-avoiding Krylov subspace methods[END_REF], our implementation is limited to one µ-step Newton basis. However, we show in the previous section how an augmented basis can be dened for more than one µ-step basis. For the same restart length, CA-GMRES(µ, t) and GMRES(µ • t) produce the same convergence behavior. AGMRES(µ • t, r) is more likely to produce a faster convergence than these two approaches when the convergence rate is aected by the restarting procedure. So far, the algorithm starts with an initial approximation of the solution vector x 0 (in practice, we use a zero vector), the size m of the Krylov basis, the maximum number of iterations itmax allowed and the desired accuracy . The remaining input values are used for the augmented basis: the number of eigenvectors r that are added at each step; the parameters l, rmax, smv and bgv for the adaptive strategy, see section 3.2.6. The main steps of the algorithm are the computation of the shifts (steps 1 and 3), the generation of the augmented Newton basis at step 6 and its orthogonalization in section 8. The approximate solution is updated at step 12. At step 21, we update the Algorithm 1 AGMRES(m, r) : Augmented Newton-basis GMRES Require: x 0 , m, itmax, , r, l, rmax, smv, bgv; Compute the basis vectors of equation 3.8 [START_REF] Sidje | Alternatives for parallel Krylov subspace basis computation[END_REF][START_REF] Erhel | A parallel GMRES version for general sparse matrices[END_REF] to get K m+1 and Tm .

7:

Compute Kr and D r using the second part of equation 3.8.

8:

Compute the QR factorization K m+1 Kr = V s+1 R s+1 . 9: Compute end if 24: end while eigenvectors to be added in the Newton basis. The adaptive strategy is implemented in steps 15-23 All these steps are explained in the next sections.

Derive W s = V m U r .
x s = x 0 + M -1 W s y s r s = b -Ax s , it ← it + s 13: if( r s < or it > itmax)

Computation of the shifts. The generation of the Krylov subspace

with the Newton polynomials uses the scalars λ j , j = 1, . . . , m to produce a stable basis. Bai et al [START_REF] Bai | A Newton basis GMRES implementation[END_REF] show that a good choice would be to use the eigenvalues of B numbered according to the following modied Leja order (see [START_REF] Reichel | Newton interpolation at Leja points[END_REF]):

|λ 1 | = max j=1,...,m |λ j | j k=1 |λ j+1 -λ k | = max l=1,...,m j k=1 |λ l -λ k |, j = 1, . . . , m -1 (3.29)
In practice, the spectrum of B is not available and very expensive to compute. In this situation, the Ritz values of B which are the eigenvalues of the Hessenberg matrix H m in equation 1.6 are used. This implies that m steps of the Arnoldi process should be performed to get these values. At step 1, we perform one cycle of the Arnoldi-GMRES method. From this step, we get an approximation of the solution x m and the associated residual r m . This vector is used as the initial search direction for Newtonbasis GMRES from step 5. At step 3, each process computes the eigenvalues of its own copy of the Hessenberg matrix H m and order them with the Leja ordering. This step uses so far the parallelism inside the matrix-vector product and the preconditioning operation. But it requires global communication as pointed out in section 1. Note that when m gets large, it may be expensive to perform this step of Arnoldi-GMRES. The cost here is compared to one step of the Newton basis GMRES mainly in terms of granularity and the volume of MPI messages. In practice, we use a small value of m to show the benets of augmenting the basis. Nevertheless, if a large basis should be used, a solution could be to use a µ-step basis as explained in the previous section. Another solution, as advised by Philippe and Reichel [START_REF] Philippe | On the generation of Krylov subspace bases[END_REF] is to perform one cycle of Arnoldi-GMRES with a smaller basis to get a subset of these values. From this subset, a convex hull is dened and continuously updated with new values collected during the Newton-basis GMRES iterations.

Computation of the Newton basis with scaling factors. The rst

m + 1 vectors of K s+1 can be generated using Algorithm 1.1 in [START_REF] Sidje | Alternatives for parallel Krylov subspace basis computation[END_REF]; then it is easy to generate the last r vectors from U r . Note that when a particular λ j+1 is complex, and assuming that Im(λ j+1 ) > 0 (This is always possible from the complex conjugate pairs and the modied Leja ordering), the complex arithmetic is avoided by writing the rst part of equation 3.8 as

σ j+1 k j+1 = (B -Re(λ j+1 )Ik j σ j+2 k j+2 = (B -λ j+1 I)(B -λj+1 I)k j = σ j+1 (B -Re(λ j+1 )I)k j+1 + Im(λ j+1 ) 2 k j .
In this case, the matrix T ∈ R s+1×s of equations 3.5 and 3.22 is tridiagonal. So far, the scalars σ j , j = 1, . . . , m in equation 3.5 are used to control the growth of the vectors {k j } j=m j=1 . The common choice is to take σ j = k j . The parallelism inside this step is through the preconditioning and the parallel matrix-vector operations (AM -1 -λI)k ≡ A(M -1 k) -λk. When σ j = k j , then there are (m + r) global communications, which are far less than the 1 2 (m 2 +3m) global communications in the Arnoldi process. With some particular cases, this norm can be computed distributively. When using for instance the explicit formulation of multiplicative Schwarz, the basis vectors are computed in a pipeline accross all the subdomains. Each process is thus able to compute its own contribution to the norm and the basis vectors are normalized a posteriori [START_REF] Atenekeng-Kahou | Parallélisation de GMRES préconditionné par une itération de Schwarz multiplicatif[END_REF][START_REF] Wakam | Parallel GMRES with a multiplicative Schwarz preconditioner[END_REF]. When the size of the basis is small enough, the rows and columns of the matrix can be equilibrated and no scaling, thus no global communication, should be needed during the computation of the basis vectors [START_REF] Hoemmen | Communication-avoiding Krylov subspace methods[END_REF].

3.2.3. Orthogonalization of the basis. After the basis vectors are computed, they should be orthogonalized between each other at step 8 of Algorithm 3.2 to produce the orthogonal system V s+1 . At the end of the step 6, the vectors K s+1 are distributed on all processors as a contiguous blocks of rows which is equivalent to the classical 1D rowwise partitioning for the matrix-vector products. Any algorithm for the parallel dense QR factorization can now be used to orthogonalize the system K s+1 . In our implementation, we use the RODDEC algorithm described in [40, section 4.2]. This method performs rst a Householder orthogonalization on each block of rows. This is done in a perfect parallel phase by all the processes having the rows. After that, the Givens rotations are used to annihilate the blocks below the rst one. During this second step, the processors are placed on a ring topology and each process sends the required data on this ring. This step requires O(m 2 ) point-to-point messages on a ring topology and the average message length is (m + 1)/2 double precision number. The approach advocated in [START_REF] Da Cunha | New parallel (rank-revealing) QR factorization algorithms[END_REF] maps the processors on a binary tree to eliminate the extra triangular factors of RODDEC. The TSQR algorithm of Demmel et al [START_REF] Demmel | Communication-optimal parallel and sequential QR and LU factorizations[END_REF] which gives a more general divide-and-conquer algorithm can be used as well at this step. It requires O(log(P )) MPI messages where P is the total number of MPI processes sharing the system K s+1 .

3.2.4. Updating the current approximation. At the end of the QR factorization, the triangular matrix R s of equation 3.11 is usually available on one process. In the RODDEC algorithm, it is available in the last process. It can be broadcasted to all other processes such that the steps 10 and 11 are done by all the processes. When the number of MPI processes gets large, it is more ecient to perform these steps on the last process and to broadcast only the result of the least-squares problem at step 11. In our implementation, we choose to send a copy of the matrix since it is required by all processes to update the eigenvectors, see section 3.2.5. So far, the Hessenberg matrix Hs is assembled from R s and Ts using a modication of Algorithm 1.2 in [START_REF] Sidje | Alternatives for parallel Krylov subspace basis computation[END_REF]. The modication allows to take into account the scaling factors of the augmented vectors in the basis. A QR factorization is perfomed on the output Hessenberg matrix to solve the least-squares problem in the minimization step. The LAPACK routine dgeqrf is used for this purpose. The output solution is used to compute the new approximate solution at step 12. Note that since we are using right preconditioning, we can obtain an estimate of the true residual norm without explicitly computing the residual vector r s . Nevertheless, at the time of restart, we need r s for the new search direction.

3.2.5. Updating the eigenvectors. When the iterative process starts at line 16 of Algorithm 3.2, the eigenpairs (u j , λ j ) of B ≡ AM -1 are approximated from the rst GMRES(m) cycle with a standard projection technique as follows :

V T m (B -λ j I)V m g j = 0
(3.30) leading to the eigenvalue problem H m g j = λ j g j .

(3.31) The Ritz values λ j , j = 0, . . . , m -1 are used as shifts for the Newton basis and the vectors u j = V m g j , j = 0, . . . , r -1 corresponding to the r smallest eigenvalues are used to augment the Newton basis.

Then to change the vectors U r at step 21, we use a Rayleigh-Ritz procedure. Indeed, as advised by the previous studies [START_REF] Chapman | Deated and augmented Krylov subspace techniques[END_REF][START_REF] Morgan | A restarted GMRES method augmented with eigenvectors[END_REF], this procedure does better at nding eigenvalues nearest zero. Using the augmented subspace C s , each extracted approximate eigenvector u is expressed as u = W s g i . Using BW s and W s , the Galerkin condition writes : (3.33)

(BW s ) T (B -λ j I)W s g j = 0.
We thus obtain a dense generalized eigenvalue problem of size s × s where (λ j , W s g j ) gives a harmonic Ritz pair of B. Multiplying F s and G s by H -T s , we get

H -T s G s = H -T s H T s αe s H s αe T m = H s + h 2 s+1,s H -T s e s e T s (3.34) H -T s F s = I s h s+1,s H -T s e s V T s+1 W s , (3.35
)

with V T s+1 W s = V T s+1 V m V T s+1 U r and V T s+1 V m = I m 0 .
The overhead here is small since it induces only r dot products of size n and a small generalized eigenvalue problem.

If the basis of C s is equal to K m U r , then the overhead is larger since it involves s dot products. Nevertheless, the numerical experiments show in most test cases that when the convergence is accelerated by deation, the time to update the approximate eigenvectors is negligible compared to the total time saved without the deation. Moreover, the adaptive strategy proposed next sets o deation only if convergence is too slow.

3.2.6. Adaptive strategy. When the desired accuracy is not achieved, the method restarts and r new approximate eigenvectors (corresponding to the eigenvalues to deate) are extracted from the s-dimensional subspace C s . This process may become expensive and is not benecial if the convergence rate is not improved enough. We thus propose an adaptive strategy which detects if the deation process will be benecial to speedup the convergence or to avoid stagnation. This approach is based upon the work by Sosonkina et al. [START_REF] Sosonkina | A new adaptive GMRES algorithm for achieving high accuracy[END_REF] which has been used successfully in another formulation of deated GMRES [START_REF] Nuentsa Wakam | Parallel adaptive deated GMRES[END_REF]. At lines 16, based on the convergence rate already achieved, we estimate the remaining number of steps (Iter) needed to reach the desired accuracy . We use a small multiple (smv) of the remaining number of steps to detect some insucient reduction in the residual norm. If it is greater than a small multiple (smv) of the number of steps allowed (itmax), then we switch to the deation. We use a large multiple (bgv) of itmax to detect a near-stagnation in the iterative process. In this case, the number of eigenvectors to augment is increased by a xed (small) value. Clearly with the parameters r, l, Iter, smv, bgv, the adaptive strategy can be sketched as follows :

• If Iter ≤ smv * itmax, the convergence rate is good enough and no more update should be done on the eigenvectors already computed. • If smv * itmax < Iter ≤ bgv * itmax, there is an insucient reduction in the residual norm and the r eigenvectors are updated for the next cycles of AGMRES. • If Iter > bgv * itmax, a stagnation may have occurred and we increase the number of eigenvalues to extract/update by a xed number l. This can be changed at runtime by the user. Otherwise stated, l = 1 in all our test cases. Note that there are more sophisticated methods to ensure that for some given values of m, GMRES(m) (and thus AGMRES(m)) will not stagnate; see for instance [START_REF] Simoncini | On a non-stagnation condition for GMRES and application to saddle point matrices[END_REF][START_REF]New conditions for non-stagnation of minimal residual methods[END_REF]. Our current stagnation test is computed a posteriori and should be mostly used to detect a very slow reduction in the residual norm. Although the proposed parameters are problem-dependent, they can be useful to avoid the stagnation if there are some previous knowledge in the convergence behavior for the problems under study. Some numerical results are given in this sense in the next section. [START_REF] Atenekeng-Kahou | Parallélisation de GMRES préconditionné par une itération de Schwarz multiplicatif[END_REF]. Numerical experiments. This section presents some numerical results to show the parallel eciency and the numerical robustness of the proposed approach. We rst present the template for all the numerical tests in section 4.1 and the test cases in section 4.2.

4.1. Test routines and implementation notes. Implementations are done using the PETSc routines and data structures [START_REF] Balay | PETSc users manual[END_REF][START_REF] Balay | PETSc Web page[END_REF]. The algorithm 3.2 has been implemented by a KSP module called AGMRES using a locally modied version of PETSc revision 3.1.p8. It uses the routines for matrix-vector product, the application of the preconditioner and the other parallel linear algebra functions. It can be used transparently with any preconditioner implemented in the package including the domain decomposition preconditioners. We use so far the Restricted Additive Schwarz (RAS) method [START_REF] Cai | A restricted additive Schwarz preconditioner for general sparse linear systems[END_REF] applied as a right preconditioner in all our tests. The main steps are outlined in Algorithm 2. Note from the step 7 of our test routine that we compare Algorithm 2 Test routine for the parallel computation of the system 1.2 using restricted additive Schwarz method and GMRES-based accelerator.

1: Read the matrix from a binary le and store it in a distributed CSR format. Read the right-hand side vector and store it accordingly. 2: Perform a parallel iterative row and column scaling on the matrix and the righthand side vector [START_REF] Amestoy | A parallel matrix scaling algorithm[END_REF]. 3: Partition the weighted graph of the matrix in parallel with PARMETIS. 4: Redistribute the matrix and right-hand-side according to the PARMETIS partitioning. 5: Dene the overlap between the submatrices for the additive Schwarz preconditioner. 6: Setup the submatrices (ILU or LU factorization using MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF]). 7: Solve iteratively the system using either the KSP AGMRES (Algorithm 3.2) or the PETSc built-in KSP GMRES [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]Algorithm 4]. 8: Write the solution vector to a binary le.

AGMRES with the classical implementation of GMRES. As stated earlier, either the classical Gram-Schmidt or the modied Gram-Schmidt can be used for the Arnoldi process. The main advantage of CGS over MGS is the number of MPI messages, the amount of MPI reductions and the granularity in the computational kernel. However, a practical implementation of CGS includes a possible renement strategy to be as stable as MGS. During our numerical experiments however, this renement was not used in CGS and we did not notice any dierence between GMRES-MGS and GMRES-CGS. We therefore give the results of GMRES with CGS. Unless stated, the stopping criterion of GMRES and AGMRES is b -Ax b < 10 -10 and the maximum number of iterations is 1,000. In AGMRES, the residual norm is computed only at each outer iteration. In GMRES, it is available during each inner iteration. Note that since we are using a right preconditioner, this residual norm is obtained cheaply from the Givens rotations that are used to transform the Hessenberg matrix in equation 1.7 into a triangular matrix.

In the following, since right preconditioning is used, the number of iterations is understood as the total number of matrix-vectors products and preconditioning steps. Hence in GMRES(m), it is equivalent to the counts of A(M -1 k). In AGMRES, it is equal to the size of the augmented basis times the restart cycles. So far, AGMRES(m) refers to the algorithm 3.2 without the deation (i.e r = 0, l = 0); In AGMRES(m, r), r vectors corresponding to the smallest harmonic Ritz values are added to the basis and replaced at each restart; With the adaptive strategy, when it is necessary, r is increased by l = 1 (otherwise stated).

In these numerical experiments, the basis of the augmented subspace C s is equal to K m U r .

Test problems. The matrices of tests arise from industrial applications in

uid dynamics and from convection-diusion problems. The main characteristics are listed in Table 4.1 The problems IM07R, VV11R and RM07R arise from design optimization in computational uid dynamics simulations. They are provided by the FLUOREM company, a CFD software editor 4 . Table 4.1 lists the coecient matrices with their main characteristics. The physical equations are the Reynolds-Averaged Navier-Stokes for compressible ows discretized using the nite volume methods as presented by [START_REF] Pacull | Study of ILU factorization for schwarz preconditioners with application to computational uid dynamics[END_REF]. The resulting matrix is formed of b×b blocks where b is the number of uid conservative variables (density, velocity, energy and turbulent variables). The matrix RM07R is available online in the University of Florida sparse matrix collection (see [START_REF] Davis | The University of Florida Sparse Matrix Collection[END_REF]) in the FLUOREM directory. The matrix is structurally symmetric in blocks. Regarding the values, the matrix is nonsymmetric and indenite. In [START_REF] Nuentsa Wakam | A comparative study of some distributed linear solvers on systems arising from uid dynamics simulations[END_REF][START_REF] Pacull | Study of ILU factorization for schwarz preconditioners with application to computational uid dynamics[END_REF], preliminary studies show that hybrid solvers based on GMRES and Schwarz-based preconditioners oer robust approaches to solve eciently these systems. As pointed in [START_REF] Pacull | Study of ILU factorization for schwarz preconditioners with application to computational uid dynamics[END_REF], we avoid the ILU factorization in the subdomain matrices because of its unpredictable behavior. We therefore rely on a direct solver (MUMPS) within each subdomain.

The test cases 3DCONSKY_121 and 3DCONSKY_161 correspond to the convective SkyScraper problem in [START_REF] Achdou | Low frequency tangential ltering decomposition[END_REF][START_REF] Kumar | Combinative preconditioning based on Relaxed Nested Factorization and Tangential Filtering preconditioner[END_REF]. The physical equation is given by the boundary value problem

div(a(x)u) -div(κ(x)∇u) = f in µ (4.1) u = 0 on ∂µ D (4.2) ∂u ∂n = 0 on ∂µ N (4.3)
where µ = [0, 1] 3 , ∂µ N = ∂µ\∂µ D . The tensor κ is isotropic and discontinuous. The domain contains many zones of high permeability which are isolated from each other. Let [x] denote the integer value of x then κ is given in 3D by

κ(x) = 10 3 * ([10 * x 2 ] + 1), if [10 * x i ] = 0 mod(2), i = 1, 2, 3, 1, otherwise (4.4)
The velocity eld a = (1000, 1000, 1000

) T and f = x 2 1 + x 2 2 + x 2
3 . The discretization is done using P2-type nite element methods in the Freefm++ 5 package. We consider 4 www.uorem.com/en/softwares/optimization/turb-opty-cfd 5 http://www.freefem.org/++/index.htm a uniform grid with n × n × n nodes and we choose = 121 and 161. During our numerical experiments, we rely on the ILU (1) factorization to approximate the solutions on the subdomains induced by the additive Schwarz method. 4.3. Platform of tests. Experiments are done on a distributed memory supercomputer Vargas6 which has 3,584 Power6 CPUs. Each Power6 CPU is a dual-core 2-way SMT with a peak frequency at 4.7 GHz. The computer is made of 112 nodes connected through an Inniband network. Each node has 32 Power6 CPUs that access 128GB of local memory in a non-uniform way (hardware NUMA nodes). The memory accessed by a single MPI process is limited to 3.2GB for the data and 0.5GB for the stack. [START_REF] Atenekeng-Kahou | Parallélisation de GMRES préconditionné par une itération de Schwarz multiplicatif[END_REF].4. Analysis of convergence with matrix RM07R. We rst consider the large test case RM07R from the FLUOREM collection and compare AGMRES(m,r) with GMRES(m). We rst use AGMRES without deation (r = 0). The goal is to conrm that the two methods have the same convergence behavior for a reasonable restart length. After that, we show the benets of using the deation when the restart length in AGMRES is small and when the number of subdomains increases. We nish this section by giving the benets of using an adaptive strategy.

In Figure 4.1, we plot the convergence of GMRES(m) and AGMRES(m) with three restart lengths, m = 32, 48 and 64. The number of subdomains is 32 and the LU factorization is used within the subdomains. The rst remark from Figure 4.1 is that there is no real dierence between the residual norm obtained from the two strategies. Secondly, the convergence curve of GMRES(m) indicates a periodic stagnation in the iterative process. These ticks occur at the time of restart and are more visible when m is small, hence the larger number of iterations. These ticks suggest that some information is lost at the time of restart and that the augmented basis could be benecial to improve the convergence rate on these cases. The other test cases give similar behaviors. Now we show the impact of deation by augmenting the basis. In Figure 4.2, we give the convergence history of GMRES(m) and AGMRES(m, r) with m = 32, 48 and r = 2, that is we compute two approximate eigenvectors at each restart and we use a basis of size s = m + 2. The number of subdomains is still 32. The adaptive strategy is not used at this point. It can be clearly noticed that adding only two eigenvectors in the basis is sucient to speedup the convergence in AGMRES. For instance, GMRES(32) requires 886 iterations while AGMRES(32,2) needs almost 272 iterations. When we increase the restart length to 48, GMRES benets greatly from that and requires almost 355 iterations to reach the desired accuracy while AGMRES(48,2) needs 250 iterations. The general notice here is that AGMRES(32,2) and AGMRES(48,2) have almost the same convergence rate while GMRES is more sensitive to the restart length. This is more visible when the number of subdomains vary.

The robustness of Schwarz preconditioners decreases as the number of subdomains increases. GMRES will thus require more and more iterations, particularly if the restart length is xed. We show this behavior in Figure 4.3, where the restart length is xed and the number of subdomains is increased. Clearly, as expected, the number of iterations in GMRES increases as we add more subdomains. For instance, GMRES(32) requires 886 iterations with 32 subdomains. With 64 subdomains, this number reaches 1000 iterations without reaching the prescribed tolerance of 10 -10 . In AGMRES(m, r), there is no such dierence. As we increase the number of subdomains, we observe that the convergence rates remain quite close. Indeed, AGM-RES(32,2) requires respectively 272 and 311 iterations for 32 and 64 subdomains. The fact that the number of iterations increases only slightly when increasing D has a great impact to the scalability of AGMRES. We give the timing results in the next section. In GMRES, a better convergence rate can be obtained if the restart length is increased as a function of the number of subdomains. We show in Table 4 and AGMRES(m, 4) as a function of the number of subdomains in the restricted additive Schwarz can be divided into three parts:

1. With 8 subdomains, GMRES needs less iterations than AGMRES for all values of the restart length. Note that this dierence is mainly due to the fact that the stopping test is computed only at each outer cycle in AGMRES. The accuracy achieved in AGMRES for these cases is always better. Typically, AGMRES gives an accuracy of 10 -14 in the computed residual while it is 10 -11 in GMRES.

2. For the same reasons, AGMRES needs more iterations than GMRES for 16 subdomains and large value of m. However for small values of m, AGMRES is cleary better than GMRES. 3. From 32 subdomains, AGMRES needs less iterations than GMRES for all restart lengths. The dash in GMRES [START_REF] Wakam | Parallel GMRES with a multiplicative Schwarz preconditioner[END_REF] for 64 subdomains denotes that the desired accuracy has not been reached within the 1,000 iterations allowed.

On the contrary, it requires almost 300 iterations with AGMRES(m, r) to converge. Thus, the main empirical conclusion from these experiments and others not reported here is that AGMRES is less sensitive to the restart length and the number of subdomains than GMRES. On the other hand, AGMRES is rather sensitive to the number of extracted eigenvectors. As for the basis size, it is dicult indeed to know how many vectors should be added to the basis to improve the convergence. If r is very large, the process of updating the eigenvectors could add more overhead. If r is small, the deation could not be benecial. The proposed adaptive strategy provides a trade-o between these two bounds. If some information about the convergence behavior has been collected before, then it can be used to dene the smv and bgv parameters in the adaptive strategy. Our goal is to show that this technique can be used to speedup the convergence by adaptively adjusting the frequency and the number of extracted eigenvalues. We take the smallest restart length m = 32, a large number of subdomains and the smallest number of harmonic Ritz vectors r = 1. Yet, we know from D = 16 that GMRES [START_REF] Wakam | Parallel GMRES with a multiplicative Schwarz preconditioner[END_REF] and thus AGMRES [START_REF] Wakam | Parallel GMRES with a multiplicative Schwarz preconditioner[END_REF] needs roughly 254 iterations. From the adaptive strategy, we still set the maximum number of iterations itmax = 1, 000 but now we set smv = 0.1 and bgv = 0.2. As explained in section 3.2.6, smv × itmax denes the lower bound below which it is not benecial to use an augmented basis, and bgv ×itmax denes the upper bound beyond which a slow convergence rate is expected and some action should be done. In this last case, we increase r by a xed value l. We take l = 2 in this case. It can be seen that when r = 1 and without adaptive strategy, the augmented basis does not contain enough spectral information to speed up the convergence. When r is adaptively increased, the basis recovers more and more spectral information and the convergence rate gets better. At the coonvergence for instance, r = 3 and r = 5 respectively for AGMRES [START_REF] Jalby | Stability analysis and improvement of the block Gram-Schmidt algorithm[END_REF][START_REF] Achdou | Low frequency tangential ltering decomposition[END_REF] and AGMRES [START_REF] Wakam | Parallel GMRES with a multiplicative Schwarz preconditioner[END_REF][START_REF] Achdou | Low frequency tangential ltering decomposition[END_REF] with the adaptive strategy while r = 1 for the non adaptive strategy. It should be noted however that the actual limitation of the proposed adaptive strategy is the choice of the right values of smv and bgv. It is heuristic and problem-dependent. Nevertheless, if there are some experimental knowledge about the convergence of GMRES on similar problems, a good interval can be set with smv and bgv around itmax to detect a near-stagnation and switch to the augmented basis. [START_REF] Atenekeng-Kahou | Parallélisation de GMRES préconditionné par une itération de Schwarz multiplicatif[END_REF].5. Convergence analysis for the other matrices. Now, we report the number of iterations in 4.3 for the two test cases IM07R and VV11R. AGMRES(m, r), we use r = 2 and adaptively increase it until rmax = 5. As noted before, we see that for a xed value of m, the number of iterations increases as D increases. It increases faster in GMRES than AGMRES. We note here again that deation is needed to reach a good accuracy for a large D. With IM07R test case for instance and for 32 sudomains in the additive Schwarz, neither GMRES [START_REF] Jalby | Stability analysis and improvement of the block Gram-Schmidt algorithm[END_REF] nor GMRES [START_REF] Wakam | Parallel GMRES with a multiplicative Schwarz preconditioner[END_REF] can produce the desired accuracy while AGMRES [START_REF] Wakam | Parallel GMRES with a multiplicative Schwarz preconditioner[END_REF] requires 724 iterations to converge. "AGMRES(24,1)-Adaptive" "AGMRES(24,1)-Non-Adaptive" "AGMRES(32,1)-Non-Adaptive" "AGMRES(32,1)-Adaptive" We end this section by considering the two matrices arising from the convectiondiusion problem. In Table 4.4, we give the number of matrix-vectors. Unlike the previous test cases, GMRES here is less sensitive to the restart parameter and the variation of subdomains. Hence the augmented basis is not as benecial to the convergence as in the previous cases. Nevertheless, AGMRES is still faster than GMRES if we consider the parallel eciency. This is the aim of the next sections.

Matrix

3DCONSKY_121 3DCONSKY_161 4.6. Analysis of the CPU time. To better show the benets of using an augmented subspace approach with the Newton basis, we analyze in this section the timing results. The paramount goal when showing these results is that, as we increase the number of subdomains, we should be able to get a decrease during the iterative time. In GMRES(m) and AGMRES(m), the best way is undoubtedly to increase the restart length as well. Even then, the time will not decrease eciently because of the negative eects of the restarting procedure and the weakness of one-level Schwarz preconditioner. In AGMRES(m, r), only a few extracted Ritz vectors are sucient to decrease the time and obtain a signicant eciency.

In Table 4.5, We compare GMRES(m), AGMRES(m) and AGMRES(m, r) on the test case RM07R by varying the number of subdomains D, the restart length m and the number of harmonic Ritz values r which increases adaptively until rmax = 5. The number of MPI processes is equal to the number of subdomains. The total time is the CPU time required to perform all the steps in Algorithm 2. The iterative time is the time spent in the step 7. The setup time is the dierence between the two times. It is independent of the method and of m. It decreases when D increases because the subdomains become smaller and the LU factorizations are faster. Thus, we concentrate from now on the iterative time. The time per iteration is the time of one cycle divided by the number of matrix-vectors products in the cycle, which is m or m + r. It includes the time to compute the orthonormal basis (with Arnoldi GMRES or the QR factorization for AGMRES) and the time to update the eigenvectors U for AGMRES(m, r). The iterative time is thus the product of the time per iteration by the number of iterations. The behaviors of both GMRES(m) and AGMRES(m) are similar. Increasing m has two opposite eects: it decreases the number of iterations (in some cases, the number of cycles remain the same for AGMRES) and increases the time per iteration, because of the orthogonalization steps. Thus, in most cases, there is an optimal value of m, which depends on D, with a minimal iterative time. Increasing D has also two opposite eects, but in the reverse way: it increases the number of iterations and decreases the time per iteration, thus there is in general an optimal value of D, which depends on m. Even though their behavior are similar, AGMRES(m) clearly performs faster than GMRES(m), for all but one congurations. This is mainly due to a faster time per iteration thanks to a more ecient parallel algorithm . This is explored in the next section by analyzing the communication volume.

The objective of deation in AGMRES(m, r) is two-fold: to get an algorithm less sensitive to m and to increase the number of subdomains (thus the number of MPI Timing statistics for RM07R; D: Number of subdomains and number of MPI processes. Total Time: CPU elapsed time in seconds, Iter. Time: CPU time in the iterative phase. Time/Iter : average time spent in each iteration (matrix-vector product and preconditioning step). MSG: MPI messages and reductions. With AGMRES(m, r), r is adaptively increased until rmax = 5. processes). For D xed, there is still an optimal value of m but it is smaller. The iterative time decreases from D = 8 until D = 64. Thus our method allows to choose a small value of m and to reduce the CPU time with a large number of subdomains. We get indeed a more ecient parallelism because the number of iterations does not inate. Clearly, AGMRES(m, r) gives the smallest CPU time. These results are conrmed with other test cases as shown in Table 4.6 and Figures 4. It is better for GMRES(m) to choose a small number of subdomains D and a large restart m. On the contrary, it is more ecient to choose a large number of Now between GMRES and AGMRES(m, r), the previous analysis holds as well but there are two situations: when the restart length is very close to the number of subdomains, then the communication for the computation of eigenvectors may dominate if there is no substantial acceleration in the convergence rate of AGMRES(m, r). This is observed in Table 4.6 for VV11R and IM07R. With a substantial gain in the convergence rate as in Table 4.5, AGMRES(m, r) benets from that and the communication volume decreases proportionally to the number of iterations. The second situation is when the number of subdomains is very large with respect to the basis size. Even if there is no substantial acceleration in AGMRES(m, r), the kernel computations of AGMRES will produce less communication volume than that in GMRES. This is observed in Figure 4.7. As the number of subdomains increase, the dierence between the two methods are more and more distincts. Between the two situations, a netuned adaptive strategy is still required to determine whether or not to augment the basis. 5. Concluding remarks. We have proposed the AGMRES(m, r) implementation, which combines the Newton basis GMRES implementation with the augmented subspace technique. This approach benets from the high level of parallelism during the kernel computation of the Krylov basis. The proposed augmented basis reduces the negative eects due to restarting and to a large number of subdomains.

The numerical results on the VARGAS supercomputer (IBM Power 6 processors) conrm that AGMRES communicates less than GMRES and produces a faster solution of large linear systems. Moreover, on the proposed test cases, AGMRES gives a fairly good convergence rate when few eigenvectors are added to the Krylov basis. The proposed implementation is done in the PETSc package. It thus benets from the optimized routines for the usual linear algebra operations on matrices and vectors. Its object-oriented interface allows to use transparently any parallel preconditioner implemented in the package, based on algebraic domain decomposition methods or multilevel methods. It can be used indeed as a smoother for algebraic multigrid methods [START_REF] Elman | A multigrid method enhanced by Krylov subspace iteration for discrete Helmhotz equations[END_REF].

Although the proposed augmented basis behaves well on the proposed test cases, there are some cases where it may not be useful and thus expensive to use. Hence a good analysis is still needed in the adaptive strategy to avoid the computation of eigenvectors on such cases.

u∈x0+M - 1

 1 Kj (B,r0) b -Au .

10 :

 10 Compute the (s + 1) × s Hessenberg matrix Hs from equation 3.13 11: Solve y s = min βe 1 -Hs y 2 12:

return 14 :

 14 Set x 0 ← x s and r 0 ← r s ; 15: if r > 0 then 16: Iter = s * log rs /log rs r0 17: if (Iter > smv * (itmax -it)) then 18: if ((Iter > bgv * (itmax -it)) and (r < rmax) and (l > 0)) then 19: r ← r + l /*Increase the number of eigenvalues to deate*the r approximate eigenvectors u 0 , u 1 , . . . , u r-1 by r harmonic Ritz vectors of B ≡ AM -1

Fig. 4 . 1 .Fig. 4 . 2 .

 4142 Fig. 4.1. RM07R : Inuence of the restart length in AGMRES and GMRES, 32 subdomains

Fig. 4 . 3 .

 43 Fig. 4.3. RM07R: Inuence of the number of subdomains in the convergence of GMRES. The restart length is xed and the benets of the augmented basis in AGMRES is given.

Figure 4 .

 4 [START_REF] Atenekeng-Kahou | Parallélisation de GMRES préconditionné par une itération de Schwarz multiplicatif[END_REF] gives the convergence history of AGMRES(m, 1) with m = 24 and m = 32.

5

 5

  and 4.6.

Fig. 4 . 5 .

 45 Fig. 4.5. CPU Time in the iterative phase for the 3D 121 × 121 × 121 convective SkyScraper problem (Matrix size 1,771,561; Nonzeros 50,178,241);16 to 256 subdomains, ILU(1) in subdomains; m = 16; r = 1

Fig. 4 . 7 .

 47 Fig. 4.7. Amount of MPI Messages in the iterative phase of GMRES and AGMRES for the convective skyscraper problems on the 3D 121 × 121 × 121 and 161 × 161 × 161 grids

  1: Perform one cycle of GMRES(m)[START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] Algorithm 4] to nd a new approximation x m , the residual r m and the matrices H m and V m satisfying equation 1.6. if( r m < ) return 2: Set x 0 ← x m and r 0 ← r m ; β = r 0 ; 3: Compute m Ritz values {λ j } j=m j=1 of AM -1 from H m and order them with the Leja ordering[START_REF] Bai | A Newton basis GMRES implementation[END_REF]. : If (r > 0) extract r Ritz vectors U r for the augmented basis. 5: while ( r 0 > ) do

	6:

4

Table 4

 4

		.1
	Characteristics of test matrices, N:Number of rows/columns, NNZ:Nonzero entries
	Matrix	N	N N Z	geometry
	IM07R	261,465	26,872,530	3D
	VV11R	277,095	30,000,952	3D
	RM07R	272,635	37,355,908	3D
	3DCONSKY_121 1,771,561 50,178,241	3D
	3DCONSKY_161 4,173,281 118,645,121	3D

  .2 that in such case, it is still benecial to have an augmented basis in AGMRES. These results

	KSP	GMRES(m)	AGMRES(m, 2) AGMRES(m, 4)
	H H D	H 8	H m 32 48 64 32 48 64 H 93 70 57 100 98 57 105 100 57 32 48 64
	16		254 169 123 169 148 130 177 153 132
	32		886 355 220 272 250 196 212 205 200
	64		-	702 445 311 303 265 287 258 270
				Table 4.2
	RM07R: Number of iterations in GMRES(m), AGMRES(m, 2)

  Number of iterations in GMRES(m), AGMRES(m, r) as a function of the number of subdomains in the restricted additive Schwarz, r is adaptively increased until rmax = 5.

	Fig. 4.4. RM07R : Benets of the adaptive deation strategy, Restart=24 and 32, 32 subdomains
	KSP	GMRES(m)	AGMRES(m, 2)
	H D H H	H m 24 32 48 24 32 48 H VV11R
	8		251 191 147 248 172 146
	16		499 458 288 492 304 207
	32		-	957 670 641 541 516
					IM07R
	8		240 235 189 249 203 195
	16		695 623 521 378 370 316
	24		927 913 759 492 444 408
	32		-	-	833 724 629 579
				Table 4.3
	VV11R & IM07R:

Table 4 . 5

 45

	D	Algo.	Total Time Iter. Time Time/Iter MSG (×10 4 )
	8	GMRES(32)	427.3	327.33	3.52	1.74
		GMRES(48)	386.1	291.64	4.166	1.41
		GMRES(64)	358.1	264.58	4.64	1.03
		AGMRES(32)	358.5	263.08	2.74	1.3
		AGMRES(48)	369.1	271.9	2.832	1.45
		AGMRES(64)	329.4	236.76	4.228	1.23
		AGMRES(32,2)	347.4	257.11	2.624	1.32
		AGMRES(48,2)	373.1	277.98	2.837	1.48
		AGMRES(64,2)	329.4	236.76	4.228	1.23
	16	GMRES(32)	379.3	349.97	1.378	13.1
		GMRES(48)	333.1	302.66	1.791	9.05
		GMRES(64)	286.8	257.03	2.09	6.88
		AGMRES(32)	305.8	276.1	1.079	8.1
		AGMRES(48)	263.0	230.5	1.201	6.78
		AGMRES(64)	256.8	227.82	1.78	5.56
		AGMRES(32,2)	224.1	193.39	1.316	9.84
		AGMRES(48,2)	240.9	210.56	1.376	10.01
		AGMRES(64,2)	231.4	201.05	1.547	5.66
	32	GMRES(32)	573.4	557.13	0.629	96.25
		GMRES(48)	239.5	223.54	0.63	39.74
		GMRES(64)	158.4	139.2	0.633	25.38
		AGMRES(32)	273.0	256.91	0.287	54.97
		AGMRES(48)	167.1	150.84	0.393	25.93
		AGMRES(64)	131.4	114.83	0.449	19.42
		AGMRES(32,4)	91.41	75.23	0.357	31.83
		AGMRES(48,4)	94.79	79.028	0.38	33.9
		AGMRES(64,4)	99.45	83.148	0.406	32.24
	64	GMRES(32)	-	-	-	-
		GMRES(48)	214.8	204.16	0.291	227.02
		GMRES(64)	165.6	156.44	0.352	145.69
		AGMRES(32)	-	-	-	-
		AGMRES(48)	167.0	157.72	0.219	132.42
		AGMRES(64)	97.87	86.066	0.192	88.67
		AGMRES(32,4)	62.39	52.839	0.202	101.53
		AGMRES(48,4)	67.0	57.733	0.22	110.99
		AGMRES(64,4)	63.15	53.788	0.203	116.08

Table 4 . 6

 46 Timing statistics in GMRES and AGMRES for test cases VV11R and IM07R. D: Number of subdomains and number of MPI processes. Iter. Time : time spent in the iterative phase. MSG: MPI messages and reductions. With AGMRES(m, r), r is adaptively increased until rmax = 5.

	H H D	H m H H	24 Iter. Time MSG Iter. Time MSG Iter. Time MSG 32 48
						GMRES(m)
	8	92.84	2.05		68.95	1.69	77.7	1.47
	16	101.1	12.27	89.37	11.47	63.2	7.66
	32 8 16	-52.8 51.8	-1.28 7.4	31.2 AGMRES(m, 2) 38.5 34.5	22.5 1.02 4.91	29.7 40.5 28.08	18.54 1.05 3.87	VV11R
	32	38.3	25.6		31.2	22.5	29.7	18.5
						GMRES(m)
	8	76.219	2.6		73.3	2.63	63.669	2.31
	16	111.74	20.06	96.246	18.25	83.583	15.76
	32 8 16	-45.781 36.492	-1.65 21.65 AGMRES(m, 2) -40.905 34.803	-5.48 24.12	77.066 40.85 33.65	59.87 1.52 23.64	IM07R
	32	33.262	94.54	27.837	93.27	27.109	105.35

Throughout this paper, we use a zero-based numbering for all the vectors in the basis

The original method is referred to as s-step Arnoldi instead of µ-step Arnoldi but we choose µ here to dierentiate with the size s of our augmented basis.

http://www.idris.fr/su/Scalaire/vargas/hw-vargas.html

Acknowledgments. This work was funded by the French National Agency of Research under the contract ANR-TLOG07-011-03 LIBRAERO. Numerical experiments have been done on the VARGAS supercomputer from GENCI-IDRIS (Grand Equipement National de Calcul Intensif -Institut du Développement et des Ressources en Informatique Scientique). Preliminary tests have been carried out using the GRID'5000 experimental testbed (https://www.grid5000.fr). We are very grateful to Bernard PHILIPPE and François PACULL for many suggestions and helpful discussions during this work. We would like to thank Guy A. Atenekeng for providing us the test matrices of the Convection-diusion equation. We are grateful to Roger B. Sidje for giving us the implementation of the RODDEC method.

In matrix form, we get

The matrices T (1) µ and D r are analogous to the matrices in equation 3.5. At this point, to avoid loss of orthogonality, the vectors K (1) µ+r+1 = [k µ , k µ+1 , . . . , k µ+m+r ] should be orthogonalized against the previous vectors V (0) µ . This can be done with a block Gram-Schmidt method which is equivalent to write 3

µ+r+1 .

(3.23)

Note that the same stability problems may arise just as in the classical Gram-Schmidt process. We discuss on this issue at the end of this subsection. So far, the vectors

µ+r are orthogonal to the basis vectors but it remains to orthogonalize them between each other. This can be done with a dense QR factorization to produce

µ+r+1 .

(3.24)

From equations 3.20 and 3.22, we get

where

Knowing that a QR factorization update have been perfomed on K (1) µ+r+1 , we get from equations 3.23 and 3.24 that,

(3.26) Substituting 3.26 in 3.25, we get

where

and

Cµ+r

.

From the fact that K (1) µ+r+1 e 1 is orthogonal to V (0) µ and that R (1) µ+r+1 e 1 = e 1 , we get

µ+r+1 Cµ+r (3.28) 3 Note that the rst vector kµ is already orthogonal to V (0) µ but we choose to orthogonalize it again. subdomains D and a small restart m with our method AGMRES(m, r). Clearly, AGMRES(m, r) is faster than GMRES(m). In order to compare the methods with similar memory requirements, we choose m = 24 for AGMRES and m = 48 for GMRES, since AGMRES needs to store the two systems W s and V s . For all but one values of D, AGMRES(24, r) is faster than GMRES(48), for both matrices VV11R and IM07R. It is also true for AGMRES(32, r) compared with GMRES(64) for the matrix RM07R.

Analysis of parallelism. The other advantage of AGMRES over GM-

RES is the communication volume. In Tables 4.5 and 4.6 and in Figure 4.7, we have reported the number of MPI messages exchanged. The counts are done on the Send/Receive routines as well as the collective communications (Reduce and Broadcast). We do not take into account the MPI message lengths. It appears rst that the number of messages is a function of the number of subdomains. This is generally reduced by allowing many subdomains to be assigned to a unique CPU. In the problems under study, this is not feasible in practise as it will induce more iterations as the subdomains increase. The communication volume is obviously proportional to the number of iterations as well. The second observation is that AGMRES communicates less than GMRES for the same number of subdomains and the same basis size. As more subdomains are used, the gap between the two methods increase. For instance, in Table 4.5, GMRES on 64 subdomains produces nearly a ratio of 1.5 more messages than AGMRES. In the augmented basis, the situation is dierent. At each cycle, AGMRES(m, r) communicates more than AGMRES(m) because of the computation of the eigenvectors. However, since a substantial number of iterations is saved by using the augmented basis, we observe actually a better communication in AGMRES(m, r).