Laura Grigori

Sophie Moufawad

Frederic Nataf

Enlarged Krylov Subspace Conjugate Gradient Methods for Reducing Communication rapport de recherche Inria

Keywords: algèbre linéaire, méthodes itératives ˚INRIA minimizing communication, linear algebra, iterative methods

Dans cet article, nous présentons deux nouveau methodes iterative pour la résolution des systèmes linéaires d'équations de très grande taille en minimisant les communications. Ces deux methodes sont basées sur l'enrichissement de l'espace de Krylov en décomposant le domaine de A.

Introduction

Krylov subspace methods are among the most practical and popular iterative methods today. They are polynomial iterative methods that aim to solve systems of linear equations (Ax " b) by finding a sequence of vectors x 1 , x 2 , x 3 , x 4 , ..., x k that minimizes some measure of error over the corresponding spaces

x 0 `Ki pA, r 0 q, i " 1, ..., k where x 0 is the initial iterate, r 0 is the initial residual, and K i pA, r 0 q " spantr 0 , Ar 0 , A 2 r 0 , ..., A i´1 r 0 u is the Krylov subspace of dimension i. Conjugate Gradient (CG) [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF], Generalized Minimal Residual (GMRES) [START_REF] Saad | Gmres: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF], bi-Conjugate Gradient [START_REF] Lanczos | Solution of systems of linear equations by minimized iterations[END_REF][START_REF] Fletcher | Conjugate gradient methods for indefinite systems[END_REF], and bi-Conjugate Gradient Sabilized [START_REF] Van Der | Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems[END_REF] are some of the most used Krylov subspace methods.

These methods are governed by Blas1 and Blas2 operations as dot products and sparse matrix vector multiplications. Parallelizing dot products is constrained by communication since the performed computation is negligible. If the dot products are performed by one processor, then there is a need for a communication before and after the computation. In both cases, communication is a bottleneck. This problem has been tackled by different approaches. First, block methods that solve system with multiple right-hand sides AX " B were introduced, as block CG [START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF]. Then, s-step methods that compute s basis vectors per iteration were proposed, examples are s-step CG [START_REF] Van Rosendale | Minimizing inner product data dependence in conjugate gradient iteration[END_REF][START_REF] Chronopoulos | s-step Iterative Methods For Symmetric Linear Systems[END_REF] and s-step GMRES [START_REF] Walker | Implementation of the gmres method using householder transformations[END_REF][START_REF] Erhel | A parallel gmres version for general sparse matrices[END_REF]. Both methods, block and s-step, use Blas2 and Blas3 operations. Recently, communication avoiding methods, based on s-step methods, that aim at avoiding communication at the expense of performing some redundant flops were introduced, as CA-CG, CA-GMRES [START_REF] Mohiyuddin | Minimizing communication in sparse matrix solvers[END_REF][START_REF] Hoemmen | Communication-Avoiding Krylov Subspace Methods[END_REF] and CA-ILU0 preconditioner [START_REF] Grigori | Communication Avoiding ILU0 Preconditioner[END_REF]. Another approach is to hide the cost of communication by overlapping it with other computation, like pipelined CG [START_REF] Demmel | Parallel Numerical Linear Algebra[END_REF][START_REF] Gropp | Update on Libraries for Blue Waters[END_REF] and pipelined GMRES [START_REF] Ghysels | Hiding Global Communication Latency in the GMRES Algorithm on Massively Parallel Machines[END_REF].

In this paper we introduce a new approach that consists of enlarging the Krylov subspace by a maximum of t vectors per iteration. First, the input matrix is partitioned into t sub-domains by using a graph partitioning algorithm. At the beginning of the iterative method, the residual is split into t vectors corresponding to the t sub-domains. Then, the obtained t vectors are multiplied by A at each iteration to generate t new basis vectors. The obtained enlarged Krylov subspace K t,k pA, r 0 q is a superset of the Krylov subspace K k pA, r 0 q, K k pA, r 0 q Ă K t,k`1 pA, r 0 q. Thus it is possible to search for the solution of the system Ax " b in K t,k pA, r 0 q instead of K k pA, r 0 q. Moreover, we show in this paper that the enlarged Krylov projection subspace methods lead to faster convergence in terms of iterations and parallelizable algorithms with less communication, with respect to Krylov methods.

In this paper we focus on Conjugate Gradient (CG) [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF], a Krylov projection method for symmetric (Hermitian) positive definite matrices, which was introduced by Hestenes and Stiefel in 1952 (section 2.1). After giving a brief overview of related existing CG methods (section 2) such as block-CG [START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF], coop-CG [START_REF] Bhaya | A cooperative conjugate gradient method for linear systems permitting multithread implementation of low complexity[END_REF], and MSD-CG [START_REF] Gu | Multiple search direction conjugate gradient method I: methods and their propositions[END_REF], we discuss two new versions of Conjugate Gradient (section 3). The first method, multiple search direction with orthogonalization CG (MSDO-CG), is an adapted version of MSD-CG [START_REF] Gu | Multiple search direction conjugate gradient method I: methods and their propositions[END_REF]. MSD-CG has the same structure as the classical conjugate gradient method where first t new search directions are defined on the t subdomains, then the t step lengths are obtained by solving a t ˆt system, and finally the solution and the residual are updated. But unlike CG, the search directions are not A-orthogonal. Thus, in MSDO-CG we A-orthonormalize the search directions, to obtain a projection method that guarentees convergence at least as fast as CG. The idea of using more than one search direction was also exploited in Rixen's thesis [START_REF] Rixen | Substructuring and Dual Methods in Structural Analysis[END_REF] for two subdomaines in the context of domain decomposition methods, and further developed in [START_REF] Gosselet | Simultaneous-FETI and related block strategies: robust domain decomposition methods for engineering problems[END_REF]. The second method that we propose here, long recurrence enlarged CG (LRE-CG), is similar to GMRES in that we build an orthonormal basis for the enlarged Krylov subspace rather than finding search directions. Then, we use the whole basis to update the solution and the residual. We show that this method is a projection method and hence should converge at least as fast as CG. We compare the convergence behavior of both methods using different Aorthonormalization and orthonormalization methods and then we compare the most stable versions with CG and other related methods (section 4 ).

We have tested our methods on matrices arising from the dicretization of 2D poisson equations (POIS-SON2D), 3D elasctisity equations (ELASTICITY3D), and 2D and 3D convection-diffusion equations such as NH2D, SKY2D, SKY3D, and ANI3D as discussed in section 4. Both methods converge faster than CG in terms of iterations, but LRE-CG converges faster than MSDO-CG since it uses the whole basis to update the solution rather than only t search directions. And the more subdomains are introduced or the larger t is, the faster is the convergence of both methods with respect to CG in terms of iterations. For example, for t " 64 the MSDO-CG and LRE-CG methods converge in 75% to 82% less iteration with respect to CG for the matrices NH2D, POISSON2D, and ELASTICITY3D, and 95% to 98% less iteration with respect to CG for the matrices SKY2D, SKY3D, and ANI3D. But increasing t also means increasing the memory requirements. Thus, in practice, t should be relatively small, depending on the available memory, on the size of the matrix, and on the number of iterations needed for convergence, as explained in section 4. We present the parallel algorithms along with their expected performance based on the estimated run times in section 5. And in section 6, we introduce the preconditioned version with its convergence behavior.

Overview of Existing Conjugate Gradient (CG) Methods

The Krylov projection methods find a sequence of approximate solutions x k (k ą 0) of the system Ax " b, and are defined by the following two conditions:

1. Subspace condition: x k P x 0 `Kk pA, r 0 q 2. Petrov-Galerkin condition: r k K L k ðñ pr k q t y " 0, @ y P L k where x 0 is the initial iterate, r 0 is the initial residual, K k pA, r 0 q " spantr 0 , Ar 0 , A 2 r 0 , ..., A k´1 r 0 u is the Krylov subspace of dimension k, and L k is a well-defined subspace of dimension k. The classical conjugate gradient is a Krylov projection method where L k " K k pA, r 0 q.

In this section we briefly introduce the conjugate gradient versions related to our MSDO-CG and LRE-CG versions, starting with the 1952 Hestenes and Stiefel version (section 2.1). Since then, many differenet versions of CG have been introduced (refer to [START_REF] Golub | Some History of the Conjugate Gradient and Lanczos Algorithms: 1948-1976[END_REF] for a historical overview of CG till 1976). In 1980 O'Leary introduced a Block CG version [START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF] that solves a system with multiple right-hand sides AX " B (section 2.2). The cooperative-CG [START_REF] Bhaya | A cooperative conjugate gradient method for linear systems permitting multithread implementation of low complexity[END_REF] which was recently introduced, solves the system Ax " b by starting with t distinct initial guesses. This is equivalent to solving the system AX " b ˚1t (algorithmically very similar to Block CG ) where 1 is a vector of ones of size t. The authors also present a parallel implementation that needs 2 to 3 synchronizations per iteration. This method has faster convergence than CG (section 2.3). The multiple search directions CG (MSD-CG) [START_REF] Gu | Multiple search direction conjugate gradient method I: methods and their propositions[END_REF] solves Ax " b by decomposing A's domain into t subdomains and defining a search direction on each of the t subdomains. Then x k " x k´1 `Pk α k , where P k is a matrix containing all the t search directions and α k is a vector of size t (section 2.4). Unlike CG, block CG and coop CG, MSD-CG does not have the A-orthogonality condition of the search directions, i.e. P t k AP i is not equal to zero for all i not equal to k. Hence it is not a projection method. This causes MSD-CG to have slower convergence than CG, and in some cases not to converge at all. That is why in multiple search directions with orthogonalization CG (MSDO-CG), after defining a search direction on each of the t subdomains, we A-orthonormalize the search directions and this leads to better convergence than CG (section 3.2).

Note that in this paper we use matlab notation for matrices and vectors. For example, given a vector p of size n ˆ1 and a set of indices δ, ppδq is the vector formed by the subset of the entries of p whose indices belong to δ. For a matrix A, Apδ, :q is a submatrix formed by the subset of the rows of A whose indices belong to α. Similarly, Ap:, αq, is a submatrix formed by the subset of the columns of A whoses indices belong to α. And Apα, βq " rApα, :qsp:, βq, is formed by the β columns of the submatrix Apα, :q Inria

Conjugate Gradient (CG) Method

Conjugate Gradient [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF] is an iterative Krylov projection method for symmetric (Hermitian) positive definite (SPD) matrices of the form # Ax " b, A " A t , x t Ax ą 0, @x ‰ 0.

(1)

Given an initial guess or iterate x 0 , at the k th iteration CG finds the new approximate solution x k " x k´1 `αk p k that minimizes φpxq " 1 2 pxq t Ax ´bt x over the corresponding space x 0 `Kk pA, r 0 q, where k ą 0, p k P K k pA, r 0 q is the k th search direction, and α k is the step along the search direction.

The minimum of φpxq is given by φpxq " 0, which is equivalent to φpxq " Ax ´b " 0. Thus, by minimizing φpxq we are solving the system [START_REF] Bhaya | A cooperative conjugate gradient method for linear systems permitting multithread implementation of low complexity[END_REF]. As the name of the method indicates, the gradients φpx i q for all i should be conjugate. And since CG is a Krylov projection method, the residual r k " b ´Ax k should respect the Petrov-Galerkin condition r k K L k , where r k is orthogonal to some well-defined subspace L k Ď R n (or Ď C n ) of dimension k. In CG, the subspace L k is the same as the Krylov subspace K k . Thus, pr k q t y " 0, for all y P K k . Hence, the residuals form an orthogonal set, pr k q t r i " 0, for all i ă k.

Moreover, the Petrov-Galerkin condition r k K K k pA, r 0 q is equivalent to the conjugacy of the gradients φpx k q t φpx i q " 0, for alli ‰ k. Once x k has been chosen, either x k is the required approximate solution of Ax " b or a new search direction p k`1 ‰ 0 must be determined to compute the new approximation x k`1 " x k `αk`1 p k`1 . This procedure is repeated until convergence or untill the maximum number of allowed iterations has been reached without convergence. The convergence criteron is set as

||r k || 2 ď ||b|| 2 , f or some P R, where r k " b ´Ax k P K k`1 pA, r 0 q is the k th residual.
Theorem 2.1. The Petrov-Galerkin condition pr k q t y " 0, for all y P K k implies the A-orthogonality of the search directions p t i Ap j " 0, for all i ‰ j and i, j ď k. Proof. By definition, p i P K i and K i Ă K i`1 . Thus p i P K i`c for c ě 0. By the Petrov-Galerkin condition r t k´1 p i " 0 and r t k p i " 0 for i ď k ´1. Thus, r t k p i " r t k´1 p i ´αp t k Ap i " 0 for i ď k ´1. This implies that p t k Ap i " 0 for i ď k ´1 since α ‰ 0. Therefore, the A-orthogonality of the search directions.

This theorem means that the A-orthogonality of the search directions has to be ensured or else the Petrov-Galerkin condition won't be respected. On the other hand, the search direction p k P K k is chosen according to the following recursion relation:

" p 1 " r 0 p k " r k´1 `βk p k´1 (2)
where p 1 is set equal to r 0 since the initial residual is equal to negative the gradient ´ φpx 0 q which is the steepest descent from x 0 . But p k is not set to r k´1 , the steepest descent from x k´1 for k ą 1, since the residuals are not A-orthogonal. It can be shown that the search directions defined in (2) are A-orthogonal i.e. p t k Ap i " 0 for all i ď k ´1. For i ă k ´1, we have

p t k Ap i " r t k´1 Ap i `βk p t k´1 Ap i " β k p t k´1 Ap i (3) since r t
k´1 Ap i " 0 by Petrov-Galerkin condition. In addition, r t k´1 p i " r t k´2 p i ´αk´1 p t k´1 Ap i " 0 with r t k´2 p i " 0 since i ď k ´2. Thus, p t k´1 Ap i " 0. Therefore, p t k Ap i " 0 for i ă k ´1. As for i " k ´1, r t k´1 Ap k´1 ‰ 0 and p t k´1 Ap k´1 ‰ 0 for p k´1 ‰ 0. Thus, β k " ´pr k´1 q t Ap k´1 pp k´1 q t Ap k´1 is chosen so that p t k Ap k´1 " 0 At each iteration, the step α k "

pp k q t r k´1 pp k q t Ap k " ||r k´1 || 2 2 ||p k || 2
A is chosen such that, φpx k q " mintφpx k´1 `αp k q, @α P Ru. Using the definition of α k , β k " ´pr k´1 q t Ap k´1 pp k´1 q t Ap k´1 "

||r k´1 || 2 2 ||r k´2 || 2 2 . RR n°8597

Block Conjugate Gradient (B-CG) Method

In 1980 O'Leary introduced a Block CG version [START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF] that solves an SPD system with multiple right-hand sides

# AX " B, A " A t , x t Ax ą 0, @x ‰ 0 (4)
where A is an n ˆn matrix, X P R nˆt is a block vector, and B is a block vector of size n ˆt containing the multiple right hand sides.

Starting with an initial guess X 0 P R nˆt , initial residual R 0 " B ´AX 0 , P 1 " R 0 γ 1 with γ 1 a t ˆt full rank freely chosen matrix, the B-CG searches for an approximate solution X k`1 P X 0 `Kk`1 pA, R 0 q where K k`1 pA, R 0 q " block ´spantR 0 , AR 0 , A 2 R 0 , ..., A k R 0 u is the block Krylov subspace. Every

n ˆt block Z P K k`1 pA, R 0 q is defined as Z " ř k i"1 A i R 0 ζ i where ζ i is a t ˆt matrix. By the Petrov- Galerkin condition we have that R k`1 K K k`1 pA, R 0 q. Then, R t k`1 Y " 0 for all Y P K k`1 pA, R 0 q, which implies that R t k`1 R i " 0 and R t k`1 AP i " 0 for all i ă k `1.
Then, for k ě 0 the iterates are defined similarly to CG:

X k`1 " X k `Pk`1 α k`1 P K k`1 pA, R 0 q R k`1 " R k ´AP k`1 α k`1 P K k`2 pA, R 0 q P k`2 " pR k`1 `Pk`1 β k`2 qγ k`2 P K k`2 pA, R 0 q where α k`1 " pP t k`1 AP k`1 q ´1γ t k pR t k R k q β k`2 " γ ´1 k`1 pR t k R k q ´1pR k`1 R k`1 q
Note that α k`1 is chosen such that φpX k`1 q " mintφpX k `Pk`1 αq, for allα P R t,t u. As for β k`1 , it is chosen to ensure the A-orthogonality of the P k`1 and P k (pP k`1 q t AP k " 0). Whereas γ k is a t ˆt full rank matrix that can be chosen freely to decrease roundoff errors in the implementation. Moreover, the search direction P k`1 P K k`1 pA, R 0 q of the block conjugate gradient method is A-conjugate, pP k`1 q t AY " 0, for all Y P K k pA, R 0 q. This leads to the A-orthogonality of the search direction ùñ pP k`1 q t AP i " 0, for all i ă k `1. We present the Block-CG algorithm in Appendix A ( Algorithm ?? ).

Cooperative Conjugate Gradient (coop-CG) Method

Recently, in 2012, Bhaya et al. presented a new version of conjugate gradient which is similar in structure to the Block conjugate gradient method. The coop-CG [START_REF] Bhaya | A cooperative conjugate gradient method for linear systems permitting multithread implementation of low complexity[END_REF] solves the system Ax = b by starting with t different initial guesses and solving the same system t times in parallel, where t threads/agents cooperate to find the solution. This is equivalent to solving the system AX " b ˚onesp1, tq where X 0 is a blockvector containing the t initial guesses, R 0 " AX 0 ´b ˚1t is the block residual, P 1 " R 0 is the initial block search direction. Then the derivations and the algorithm of the coop-CG (Algorithm ??, Appendix A) are the same as the Block-CG with γ k " I.

Multiple search direction Conjugate Gradient (MSD-CG) Method

The multiple search directions CG (MSD-CG), introduced by Gu et al. [START_REF] Gu | Multiple search direction conjugate gradient method I: methods and their propositions[END_REF], solves the system Ax=b, and starts by having a decomposed domain and by defining at each iteration k a search direction p k i on each of the t subdomains (δ i , i " 1, 2, ..., t) such that p k i pδ j q " 0 for all j ‰ i. Then, the approximate solution at the k th iteration is defined as x k " x k´1 `Pk α k where P k " rp k 1 p k 2 p k 3 ... p k t s is a matrix containing all the k th search directions and α k is a vector of size t.

Inria

Given an initial guess x 0 , the residual is defined as r k " b ´Ax k for k ě 0. The first set of domain search directions is defined by the initial residual r 0 , such that p 1 i pδ i q " r 0 pδ i q for i " 1, 2, ..., t and zero otherwise. Then, for k ą 1 the domain search directions are defined as follows, p k i " T i pr k´1 q `βk i p k´1 i for i " 1, 2, ..., t where β k i is a scalar and T i is an operator that projects a vector onto the subdomain δ i (rT i pxqspδ j q " 0 for j ‰ i and rT i pxqspδ i q " xpδ i q ). The search directions block has the following sparsity pattern for all k,

P k " ¨˚0 0 0 . . . . . . . . . . . . ˚0 0 0 0 ˚0 0 . . . . . . . . . . . . 0 ˚0 0 . . . 0 0 ˚0 . . . . . . . . . . . . 0 0 ˚0 0 0 0 .. . . . . . . . . . . 0 0 0 ˚‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' nˆt .
As for α k " pP t k AP k q ´1P t k r k´1 , it is chosen such that it minimizes φpx k q " mintφpx k´1 `Pk αq, @α P R t u. Unlike CG, block CG and coop CG, MSD-CG does not have the A-orthogonality condition of the search directions, i.e. P t k AP i is not equal to zero for all i not equal to k. Thus β k " pP t k´1 AP k´1 q ´1P t k´1 Ar k´1 is chosen so that the global search direction p k " ř t i"1 p k i is A-orthogonal to the previous domain search direction p k´1 i , i.e. pp k q t AP k´1 " 0, for i " 1, 2, .., t. As for the convergence, it is shown that the rate of convergence of MSD-CG is at least as fast as that of the steepest descent method. Yet, steepest descent is known for its slow "zig-zagging" convergence. This causes the MSD-CG to have slower convergence than CG, and in some cases it does not converge at all with respect to the given stopping criteria as shown in section 4.

The New Conjugate Gradients

We will introduce two new conjugate gradient methods, MSDO-CG and LRE-CG, which are based on replacing the Krylov subspace K k with a larger subspace leading to better convergence. Thus we will first introduce the new enlarged Krylov subspace and its properties in the context of conjugate gradient methods in section 3.1. Then in section 3.2 and section 3.3 we introduce the multiple search direction conjugate gradient with orthogonalization (MSDO-CG) and the long recurrence enlarged conjugate gradient (LRE-CG).

As previously mentioned, MSDO-CG is an adapted version of MSD-CG, where the t newly defined search directions are A-orthonormalized against previous search directions and against each others. This A-orthonormalization guarentees a convergence behavior at least as good as CG. As for the LRE-CG, at each iteration, t new basis vectors are computed for the enlarged Krylov subspace. Then, rather than having short recurrences, x k is defined by all the basis vectors as in GMRES, where the basis vectors are orthonormalized. Both methods, MSDO-CG and LRE-CG, require saving at most tk vectors versus one search direction in CG. Yet LRE-CG converges faster than MSDO-CG (section 4) at the expense of solving growing systems of size tk. Several remedies to this problem are discussed in section 3.3.1.

The Enlarged Krylov Subspace

The enlarged Krylov subspace and methods are based on a partition of the unknowns, or alternatively the rows of the n ˆn matrix A. Assume that the index domain δ " t1, 2, .., nu is divided into t distinct subdomains δ i , where δ " Y t i"1 δ i . We define T i pxq to be the operator that projects the vector x onto the subdomain δ i . Let y " T i pxq, then ypδ i q " xpδ i q and zero elsewhere. Then, we define T pxq to be an operator that transforms the n ˆ1 vector x into t vectors of size n ˆ1 that correspond to the projection of x onto the subdomains δ i for i " 1, 2, .., t. If the obtained t vectors are assembled in increasing order into a block vector X, then we have Xpδ i , iq " xpδ i q for all i and zero elsewhere. We will refer to R 0 as the block containing the t vectors obtained from T pr 0 q. Note that R 0 ‰ T pr 0 q since R 0 is a matrix, whereas T pr 0 q " tT 1 pr 0 q, T 2 pr 0 q, .... , T t pr 0 qu is a set of vectors. But R 0 " rT 1 pr 0 q T 2 pr 0 q .... T t pr 0 qs, where the brackets r..s denote a matrix format. Definition 3.1. Let K t,k " spantT pr 0 q, AT pr 0 q, A 2 T pr0q, ..., A k´1 T pr 0 qu " spantT 1 pr 0 q, T 2 pr 0 q, ..., T t pr 0 q, AT 1 pr 0 q, AT 2 pr 0 q, ..., AT t pr 0 q, ..., A k´1 T 1 pr 0 q, ..., A k´1 T t pr 0 qu be an enlarged Krylov subspace of dimension k ď z ď tk generated by the matrix A and the vector r 0 , and associated to a given partition defined by δ i for i " 1, 2, .., t.

The enlarged Krylov subspaces K t,k pA, r 0 q are increasing subspaces, yet bounded. We denote by k max the upper bound k for which the dimension of the enlarged Krylov subspace K t,k pA, r 0 q stops increasing. For simplicity, we will denote the enlarged Krylov subspace generated by A and r 0 , K t,k pA, r 0 q, by K t,k , and the Krylov subspace generated by A and r 0 , K k pA, r 0 q by K k .

Theorem 3.2. The Krylov subspace K k is a subset of the enlarged Krylov subspace K t,k (K k Ă K t,k ).

Proof. Let y P K k where K k " spantr 0 , Ar 0 , .., A k´1 r 0 u. Then

y " k´1 ÿ j"0 a j A j r 0 " k´1 ÿ j"0 a j A j R 0 ˚1t " k´1 ÿ j"0 t ÿ i"1 a j A j T i pr 0 q P K t,k
since r 0 " R 0 ˚1t " rT 1 pr 0 q T 2 pr 0 q .... T t pr 0 qs ˚1t . Krylov subspace methods search for an approximate solution x k P x 0 `Kk . A corollary of theorem 3.2 is that we can search for an approximate solution x k in x 0 `Kt,k instead, since K k Ă K t,k .

In theorem 3.3, we do not use the direct sum ' since it is not guarenteed that the intersection of the two subspaces, K t,k and spantA k T 1 pr 0 q, A k T 2 pr 0 q, ..., A k T t pr 0 qu, is empty. Theorem 3.3. By definition 3.1 of the enlarged Krylov subspace, K t,k`1 " K t,k `spantA k T 1 pr 0 q, A k T 2 pr 0 q, ..., A k T t pr 0 qu If A k T v pr 0 q P K t,k for all 1 ď v ď t, then A k`q T i pr 0 q P K t,k for some 1 ď i ď t and for some q ą 0.

Proof. We prove this by induction. Base Case: Given that A k T v pr 0 q P K t,k for all 1 ď v ď t , we show that A k`1 T i pr 0 q P K t,k , where

1 ď i ď t. A k T i pr 0 q " ř k´1 u"0 ř t v"1 α u,v A u T v pr 0 q since A k T i pr 0 q P K t,k . Then A k`1 T i pr 0 q " k´1 ÿ u"0 t ÿ v"1 α u,v A u`1 T v pr 0 q " k´2 ÿ u"0 t ÿ v"1 α u,v A u`1 T v pr 0 q `t ÿ v"1 α k´1,v A k T v pr 0 q " k´2 ÿ u"0 t ÿ v"1 α u,v A u`1 T v pr 0 q `t ÿ v"1 α k´1,v p k´1 ÿ u"0 t ÿ y"1 β u,y A u T y pr 0 qq " k´1 ÿ u"0 t ÿ v"1 γ u,v A u T v pr 0 q P K t,k

Inria

Assume true for q Assume that A k`q T i pr 0 q P K t,k where 1 ď i ď t, that is A k`q T i pr 0 q " ř k´1 u"0 ř t v"1 α u,v A u T v pr 0 q Prove true for q+1 Show that A k`q`1 T i pr 0 q P K t,k

A k`q`1 T i pr 0 q " k´1 ÿ u"0 t ÿ v"1 α u,v A u`1 T v pr 0 q " k´2 ÿ u"0 t ÿ v"1 α u,v A u`1 T v pr 0 q `t ÿ v"1 α k´1,v A k T v pr 0 q " k´2 ÿ u"0 t ÿ v"1 α u,v A u`1 T v pr 0 q `t ÿ v"1 α k´1,v p k´1 ÿ u"0 t ÿ y"1 β u,y A u T y pr 0 qq " k´1 ÿ u"0 t ÿ v"1 γ u,v A u T v pr 0 q P K t,k
Given that K t,k ‰ K t,k´1 , then a corollary of Theorem 3.3 is that K t,k " K t,k`q for all q ą 0, where k max " k is the upper bound for which the dimension of the enlarged Krylov subspace stops increasing. Assume that A k T v pr 0 q P K t,k for all 1 ď v ď t, then by Theorem 3.3 A k`q T i pr 0 q P K t,k for all q ą 0 and for some 1 ď i ď t. Then for all 1 ď i ď t and for all q ą 0, A k`q T i pr 0 q P K t,k . Thus no new vector is added to the basis of K t,k`q for all q ą 0 and K t,k " K t,k`q . Moreover, since K t,k ‰ K t,k´1 then k max " k, the upper bound for which the dimension of enlarged Krylov subspace K t,k pA, r 0 q stops increasing.

Theorem 3.4. If A k T i pr 0 q P K t,k `spantA k T 1 pr 0 q, ..., A k T i´1 pr 0 q, A k T i`1 pr 0 q, ..., A k T t pr 0 qu, then A k`q T i pr 0 q P K t,k`q `spantA k`q T 1 pr 0 q, ..., A k`q T i´1 pr 0 q, A k`q T i`1 pr 0 q, ..., A k`q T t pr 0 qu for all 1 ď i ď t and q ą 0.

Proof. If A k T i pr 0 q P K t,k `spantA k T 1 pr 0 q, ..., A k T i´1 pr 0 q, A k T i`1 pr 0 q, ..., A k T t pr 0 qu, then A k T i pr 0 q " ř k´1 u"0

ř t v"1 α u,v A u T v pr 0 q `řt v"1 v‰i α k,v A k T v pr 0 q. Thus, A k`q T i pr 0 q " k´1 ÿ u"0 t ÿ v"1 α u,v A u`q T v pr 0 q `t ÿ v"1 v‰i
α j,v A k`q T v pr 0 q P K t,k`q `spantA k`q T 1 pr 0 q, ..., A k`q T i´1 pr 0 q, A k`q T i`1 pr 0 q, ..., A k`q T t pr 0 qu A corollary of Theorem 3.4 is that if t ´ik vectors of the form A k T y pr 0 q with y " i k `1, ..., t belong to the subspace K t,k `spantA k T 1 pr 0 q, A k T 2 pr 0 q, ..., A k T i k pr 0 qu, then the t ´ik vectors of the form A k`q T y pr 0 q belong to the subspace K t,k`q `spantA k`q T 1 pr 0 q, A k`q T 2 pr 0 q, ..., A k`q T ij pr 0 qu. Theorem 3.5. Let k max be the smallest integer such that K t,kmax " K t,kmax`q for all q ą 0. Then, for all k ă k max the dimension of the enlarged Krylov subspaces K t,k and K t,k`1 is stricltly increasing by some number i k and i k`1 respectively, where

1 ď i k`1 ď i k ď t.
Proof. By definition of k max , we have that for all q ą 0 K t,1 Ĺ ... Ĺ K t,kmax´1 Ĺ K t,kmax " K t,kmax`q . Then for all k ă k max , the dimension of the enlarged Krylov subspaces K t,k is stricltly increasing by some number i k ‰ 0 with respect to the dimension of K t,k´1 .

If the t new vectors are linearly independent and none of them belongs to K t,k´1 , then the t vectors are added to the basis of K t,k and dimpK t,k q " dimpK t,k´1 q `t, where i k " t and dimpq is the dimension of a subspace . In case the t new vectors are linearly dependent and none of them belongs to K t,k´1 , then only one vector is added to the basis of K t,k and dimpK t,k q " dimpK t,k´1 q `1, that is i k " 1. There are many other cases where 1 ă t ´ik ă t of the t vectors belong to K t,k´1 or are linearly dependant on the other i k vectors and K t,k´1 . Then i k vectors are added to the basis of K t,k and dimpK t,k q " pK t,k´1 q `ik , where 1 ă i k ă t.

In general, dimpK t,k q " dimpK t,k´1 q `ik , where 1 ď i k ď t. Similarly, dimpK t,k`1 q " dimpK t,k q `ik`1 , where 1 ď i k`1 ď t. Moreover, in K t,k 's basis we added i k new vectors of the form A k´1 T i pr 0 q while the other t ´ik either belong to K t,k´1 or are linearly dependant on the i k vectors and K t,k´1 . In both cases, the t´i k vectors of the form A k´1 T i pr 0 q belong to the subspace K t,k´1 spantA k´1 T 1 pr 0 q, , ..., A k´1 T i k pr 0 qu. Then by Theorem 3.4 and its corollary, the t ´ik vectors of the form A k`q T i pr 0 q belong to the subspace K t,k`q `spantA k`q T 1 pr 0 q, A k`q T 2 pr 0 q, ..., A j`q T i k pr 0 qu for q ą 0.Therefore, we have at least t´i k linearly dependent vectors added to K t,k`1 , hence i k`1 can never be greater than i k .

Theorem 3.6. Let p max and k max be such that K pmax " K pmax`q and K t,kmax " K t,kmax`q for q ą 0. Then k max ď p max .

Proof. Let K pmax " K pmax`q and A pmax`q´1 r 0 P K pmax`q where q ą 0. Then A pmax`q´1 r 0 P K pmax Ă K t,pmax , and A pmax`q´1 r 0 " ř pmax j"1 ř t i"1 α j,i A j´1 T i pr 0 q. Thus A pmax`q´1 ř t i"1 T i pr 0 q " ř pmax j"1 ř t i"1 α j,i A j´1 T i pr 0 q. Suppose that A pmax`q´1 T i pr 0 q R K t,pmax for all 1 ď i ď t. Then A pmax`q´1 ř t i"1 T i pr 0 q " ř pmax`q´1 j"1 ř t i"1 α j,i A j´1 T i pr 0 q. We may assume that there exists at least one α j,i ‰ 0 for j ą p max , then this leads to a contradiction. This implies that A pmax`q´1 T i pr 0 q P K t,pmax for all 1 ď i ď t.

Thus by definition of the T pq operator and since K p is a subset of K t,p , if K pmax " K pmax`q , then K t,pmax " K t,pmax`q . However, if K t,kmax " K t,kmax`q this does not imply that K kmax " K kmax`q . Therefore, since K t,k is a much larger subspace than K k , it is possible to reach stagnation earlier. Therefore k max ď p max .

Theorem 3.7. The solution of the system Ax " b belongs to the subspace x 0 `Kt,kmax , where K t,kmax`q " K t,kmax , for q ą 0.

Proof. The solution x sol P x 0 `Kpmax , where K pmax " spantr 0 , Ar 0 , .., A pmax´1 r 0 u and K pmax " K pmax`q for q ą 0. Since K pmax Ă K t,pmax , the solution x sol P x 0 `Kt,pmax , where p max ě k max by Theorem 3.6.

Suppose that x sol P x 0 `Kt,pmax , but x sol R x 0 `Kt,kmax . This implies that K t,kmax ‰ K t,pmax . However, by definition of k max and since k max ď p max , we have that K t,kmax " K t,pmax . This is a contradiction.

Krylov projection methods

The Krylov projection methods find a sequence of approximate solutions x k (k ą 0) of the system Ax " b from the subspace x 0 `Kk Ď R n (or Ď C n ) by imposing the Petrov-Galerkin constraint on the k th residual r k " b ´Ax k , that is r k is orthogonal to some well-defined subspace of dimension k.

We define our new enlarged Krylov projection methods based on CG by the subspace K t,k and the following two conditions:

1. Subspace condition:

x k P x 0 `Kt,k 2. Orthogonality condition: r k K K t,k
ðñ pr k q t y " 0, for all y P K t,k where K t,k is a well-defined subspace of dimension k ď z ď tk.

Inria

The minimization property

The new enlarged CG methods find the new approximate solution by minimizing the function φpxq over the subspace x 0 `Kt,k .

Theorem 3.8. If r k K K t,k , then φpx k q " mintφpxq, @x P x 0 `Kt,k u.

Proof. By the Petrov-Galerkin condition we have that r k K K t,k ùñ pr k q t y " 0, @y P K t,k pb ´Ax k q t y " 0, @y P K t,k b t y ´px k q t Ay " 0, @y P K t,k Let y " x k ´x0 P K t,k ùñ px k q t Apx k ´x0 q ´bt px k ´x0 q " 0

ùñ px k q t Ax k ´bt x k " px k q t Ax 0 ´bt x 0 ùñ φpx k q " 1 2 px k q t Ax k ´bt x k " ´1 2 px k q t Ax k `px k q t Ax 0 ´bt x 0
By showing that φpxq ě φpx k q, for all x P x 0 `Kt,k then we have proven that φpx k q " mintφpxq, @x P x 0 `Kk u.

(5)

φpxq ´φpx k q " 1 2 x t Ax ´bt x ´r´1 2 px k q t Ax k `px k q t Ax 0 ´bt x 0 s " 1 2 x t Ax ´bt z `1 2 px k q t Ax k ´px k q t Ax 0 , where z " x ´x0 P K t,k " 1 2 x t Ax ´px k q t Az `1 2 px k q t Ax k ´px k q t Ax 0 , since b t z " px k q t Az " 1 2 x t Ax ´px k q t Ax `1 2 px k q t Ax k " 1 2
px ´xk q t Apx ´xk q ě 0, since A is positive def inite Theorem 3.9. φpx k q " mintφpxq, @x P x 0 `Kt,k u if and only if ||x ˚´x k || A " mint||x ˚´x|| A , @x P x 0 `Kt,k u, where x ˚is the exact solution of (1).

Proof. f pxq " ||x ˚´x|| A " px ˚qt Ax ˚´2px ˚qt Ax `xt Ax " b t x ˚´2b t x `xt Ax " b t x ˚`2φpxq.
The minimum of f pxq is given by f 1 pxq " φpxq " 0.

Convergence analysis

The conjugate gradient method of Hestenes and Stiefel is known to converge in K iterations where K ď

n, if the matrix A P R n,n is SPD. Moreover, the k th error of CG e k " ||x ˚´x k || ď 2 ´?κ´1 ? κ`1 ¯k ||e 0 || A where κ " ||A|| 2 ||A ´1|| 2 " λmax
λmin is the L2-condition number of the matrix A. Assuming that the k th residual of the new conjugate gradient methods r k K K t,k , then by Theorem 3.8 and Theorem 3.9 we have that

||e k || A " ||x ˚´x k || A " mint||x ˚´x|| A , @x P x 0 `Kt,k u (6)
ď mint||x ˚´x|| A , @x P x 0 `Kk u since K k Ă K t,k (7)
ď ||e k || A (8) Therefore, our methods converge at least as fast as the Classical Conjugate Gradient method, assuming that the Petrov Galeriken condition is respected (r k K K t,k ). Hence, the enlarged Krylov subspace CG methods will converge in K iterations, where K ď K ď n.

Multiple search direction with orthogonalization Conjugate Gradient (MSDO-CG) Method

The MSD-CG method introduced in Gu et al [START_REF] Gu | Multiple search direction conjugate gradient method I: methods and their propositions[END_REF], can be viewed as an enlarged Krylov method where P 0 " rT pr 0 qs, and at the k th iteration p k i "

T i pr k´1 q `βk i p k´1 i for i " 1, 2, ..., t, P k " rp k 1 , p k 2 , ...., p k t s, x k " x k´1 `Pk α k and r k " r k´1 ´AP k α k with α k " pP t k AP k q ´1P t k r k´1 and β k " rβ k 1 , β k 2 , ..., β k t s " pP t
k´1 AP k´1 q ´1P t k´1 Ar k´1 . However, the P k 's are not A-othogonal implying that r k M K t,k . Thus, MSD-CG is not a projection method.

The multiple search directions with orthogonalization CG (MSDO-CG) is an enlarged Krylov projection method that solves the system (1) (Ax=b), by approximating the solution at the k th iteration with the vector x k " x k´1 `Pk α k such that φpx k q " mintφpxq, @x P K t,k u, where P k α k P K t,k , P k is an n ˆt block vector containing the t subdomain search directions, and α k is a vector of size t.

The minimum of φpxq is given by φpxq " 0, which is equivalent to Ax ´b " 0. Thus, by minimizing φpxq, we are solving the system (1). Note that since φpx k q " mintφpxq, @x P x 0 `Kt,k u, then φpx k q " φpx k´1 `Pk α k q " mintφpx k´1 `Pk αq, @α P R t u.

Once x k has been chosen, either x k is the desired solution of Ax " b, or t new domain search direction vectors P k`1 and a new approximation

x k`1 " x k `Pk`1 α k`1 are computed. Similarly to MSD-CG, P k`1 " rp k`1 1 p k`1 2 ... p k`1 t s, where p 1 i " T i pr 0 q and p k`1 i " T i pr k q `βk`1 i p k i for i " 1, 2, ..., t.
But unllike MSD-CG, MSDO-CG is a projection method. Hence, we A-orthonormalize all the search directions, P k`1 , to ensure that r k`1 K K t,k`1 as discussed in section 3.2.2. By imposing the orthogonality condition, r k`1 K K t,k`1 u, it is guarenteed that MSDO-CG converges at least as fast as CG as prooven in section 3.1.3.

This procedure is repeated until convergence. Thus, we need to find the recursion relations of r k , P k , α k , and

β k " rβ k 1 , β k 2 , ..., β k t s t .

The residual r k

By definition, the residual r k " b ´Ax k , where x k P K t,k . Thus r k P K t,k`1 . As for the recursion relation of r k , we simply replace x k by its expression and obtain the following:

r k " b ´Ax k " b ´Apx k´1 `Pk α k q " r k´1 ´AP k α k Moreover, if the orthogonality condition, r k K K t,k
, is ensured, then pr k q t r i " 0, for all i ă k. Hence, the residuals form an orthogonal set.

Theorem 3.10. The orthogonality condition pr k q t y " 0 for all y P K t,k , implies the A-orthogonality of the block search directions P t i AP j " 0, for all i ‰ j, and i, j ď k. Proof. By definition, P i P K t,i and K t,i Ă K t,i`1 . Thus P i P K t,i`c for c ě 0. By the Petrov-Galerkin condition r t k´1 P i " 0 for i ď k ´1 and r t k P i " 0. Thus, r t k P i " r t k´1 P i ´αt k P t k AP i " 0 for i ď k ´1. This implies that P t k AP i " 0 for i ď k ´1 since α ‰ 0. Therefore, the A-orthogonality of the search directions.

The domain search direction P k

By definition, the domain search direction is P k " rp k 1 p k 2 ... p k t s where p 1 i " T i pr 0 q and p k i " T i pr k´1 qβ k i p k´1 i for i " 1, 2, ..., t. p k i P K t,k for i " 1, 2, ..., t and P k α k P K t,k .

Inria

The recursion relation of P k is defined as follows P k " T pr k´1 q `Pk´1 diagpβ k q, (10) where diagpβ k q is a t ˆt matrix with the vector β k on the diagonal.

The domain search directions defined in [START_REF] Gosselet | Simultaneous-FETI and related block strategies: robust domain decomposition methods for engineering problems[END_REF] are not A-orthogonal to each others. To ensure that the orthogonality condition is valid, at each iteration k the block vector P k is A-orthonormalized against all the previous P i , where i " 1, 2, .., k ´1. Then the column vectors of P k are A-orthonormalized against each others. Thus, the obtained search directions r P k satisfy p r P k q t A r P i " 0 for all i ‰ k. Moreover, p r P k q t A r P k " I, where I is the identity matrix, assuming that the column vectors of P k are linearly independant with respect to each others and the previous directions, or alternatively none of the column vectors of r P k is zero. Note that, once P k " T pr k´1 q `Pk´1 diagpβ k q is defined, it is directly A-orthonormalized. Thus, in the sections that follow, we denote by P k the A-orthonormalized search directions and we do not use the r P k notation to be consistent with the initial definitions in the previous sections.

There are several A-orthonormalization methods. First, for A-orthonormalizing P k against all the previous P i , where i " 1, 2, .., k ´1, one can use classical Gram Schmidt (CGS), modified Gram Schmidt (MGS), or classical Gram Schmidt with reorthogonalization (CGS2) where we apply the CGS algorithm twice for numerical stability reasons. As for A-orthonormalizing P k , there are many methods that are discussed in [START_REF] Lowery | Stability Analysis of QR factorization in an Oblique Inner Product[END_REF][START_REF] Rozloznik | Numerical stability of orthogonalization with a non-standard inner product[END_REF], but we will only refer to CGS, CGS2, MGS, A-CholQR and Pre-CholeQR. We seek a combination of both A-orthonormalizations that is stable and parallelizable with reduced communication. For that reason, in section 4 we test the MSDO-CG method with the different combinations of the A-orthonormalization methods and we conclude that the MSDO-CG is numerically most stable when we use MGS, CGS2+A-CholQR, or CGS2+Pre-CholQR. In section 5.1 we discuss the parallelization of the MSDO-CG algorithm with the stable A-orthonormalization methods.

Note that in Appendix B, we discuss the A-orthonormalization using modified Gram Schmidt and classical Gram Schmidt. We also present versions of the algorithms that reduce communication along with their parallelizations. For example, Algorithm 17 is a block Gram Schmidt A-orthonormalization based on classical Gram Schmidt that A-orthonormalizes P k against previous vectors. And Algorithm 20 A-orthonormalizes P k 's vectors against each others using a classical Gram Schmidt.

Finding the expression of α k`1 and β k`1

At each iteration the step α k`1 is chosen such that φpx k`1 q " mintφpx k `Pk`1 αq, @α P R t u Let F pαq " φpx k `Pk`1 αq where φpxq " 1 2 x t Ax ´xt b. Then, F pαq "

1 2 px k `Pk`1 αq t Apx k `Pk`1 αq ´px k `Pk`1 αq t b " φpx k q `1 2 rpx k q t AP k`1 α `αt pP k`1 q t Ax k `αt pP k`1 q t AP k`1 αs ´αt pP k`1 q t b " φpx k q `1 2 rpx k q t AP k`1 α ´αt pP k`1 q t Ax k s `1 2 α t pP k`1 q t AP k`1 α ´αt pP k`1 q t r k " φpx k q `1 2 α t pP k`1 q t AP k`1 α ´αt pP k`1 q t r k , since A is SP D
The minimum of F pαq is given by

F 1 pαq " 0. ñ F 1 pαq " pP k`1 q t AP k`1 α ´pP k`1 q t r k " 0 Therefore, α k`1 " pP t k`1 AP k`1 q ´1pP t k`1 r k q . As for β k`1 , it should be chosen to ensure that P k`1 is A-orthogonal to P k . P k`1 " T pr k q Pk diagpβ k`1 q and P t k AP k`1 " P t k AT pr k q `P t k AP k diagpβ k`1 q. Since P k is an A-orthonormal matrix, P t k AP k " I, diagpβ k`1
q should be equal to ´P t k AT pr k q. But nothing guarantees that P t k AT pr k q is a diagonal matrix. So we choose β k`1 " pP t k AP k q ´1P t k Ar k which guarantees that P k`1 ˚1t is Aorthogonal to P k , similarly to MSD-CG. Moreover, in case P t k AT pr k q is a diagonal matrix, then our choice of β k`1 implies that P k`1 is A-orthogonal to P k . If t " 1, then MSDO-CG is reduced to the classical conjugate gradient.

Note that, since the vectors of P k`1 are A-orthonormalized (P t k`1 AP k`1 " I), then α k`1 and β k`1 systems are reduced to α k`1 " P t k`1 r k and β k`1 " ´P t k Ar k .

The MSDO-CG Algorithm

After deriving the recurrence relations of x k , r k , P k , α k , and β k , we present the MSDO-CG algorithm in Algorithm 1. We do not specify the A-orthonormalization methods, since this choice will be based first on the numerical stability of the method (section 4), then on its parallelization with the least communication possible (section 5.1 ).

Thus we present the MSDO-CG algorithm (Algorithm 1) and the computed flops per iteration except for the A-orthonormalizations. To reduce communication and computation in the A-orthonormalizations, be it MGS (Algorithms 14 and 15), CGS (Algorithms 18 and 21 ), A-CholQR (Algorithm 25), or Pre-CholQR (Algorithm 27), we replace

W k`1 " AP k`1 by $ & % W 1 " AP 1 W k`1 " AT pr k q `AP k diagpβq @k ě 1 " AT pr k q `Wk diagpβq

Algorithm 1 MSDO-CG algorithm Flops

Input: A, the n ˆn symmetric positive definite matrix Input: b, the n ˆ1 right-hand side; x 0 , the initial guess or iterate Input: , the stopping tolerance; k max , the maximum allowed iterations Output: x k , the approximate solution of the system Ax " b 1: r " b ´Ax 0 , ρ " ||r|| 2 2 , k " 1 2nnz `2n ´1 2: Let P 1 " T pr 0 q and W 1 " AP 1 2nnz ´pt ´1qn α " pP t k W k q ´1pP t k rq " P t k r p2n ´1qt

6:

x " x `Pk α p2t ´1qn `n 7:

r " r ´Wk α p2t ´1qn `n 8:

ρ " ||r|| 2 2 2n ´1 9:
β " ´pP t k W k q ´1pW t k rq " ´W t k r p2n ´1qt

10:

P k`1 " T prq `Pk diagpβq 2nt 11: W k`1 " AT pr k q `Wk diagpβq 2nnz ´pt ´1qn `2nt 12:
A-orthonormalize P k`1 against all P i 's for i ď k not included here 13:

A-orthonormalize P k`1 not included here 14:

k " k `1 1 15: end while The total number of flops computed sequentially after k c iterations, except for the A-orthonormalizations, is Total Flops " 4nnz ´nt `5n ´1 `kc r11nt ´2t `2n ´1 `2nnz `n `1s " 4nnz ´nt `5n ´1 `kc r11nt `3n ´2t `2nnzs « 4nnz `5n `kc r11nt `2nnzs, which is of the order of nnzk c `ntk c flops, where nnz is the number of nonzero entries in the n ˆn matrix A and t is the number of search directions computed at each iteratin.

Inria

It must be noted that since the P i 's are A-orthonormal to each others, then the t ˆt matrix P t k W k " P t k AP k is the identity matrix. Hence, solving for α k and β k is simply perfoming matrix vector multiplication.

Long Recurrence Enlarged Conjugate Gradient (LRE-CG) Method

In this section, we introduce the long recurrence enlarged conjugate gradient (LRE-CG) method which is an enlarged Krylov projection method that solves the system Ax " b by approximating the solution at the k th iteration with the vector x k " x k´1 `Qk α k P x 0 `Kt,k such that φpx k q " mintφpxq, @x P x 0 `Kt,k u, where Q k α k P K t,k and Q k is an n ˆtk matrix containing the orthonormal basis vectors of K t,k and φpxq " 1 2 x t Ax ´xt b. The LRE-CG method does not have short recurrences as MSDO-CG, but it has similarities with GMRES in that the whole basis is used to define the new approximate solution rather than t search directions. As mentioned earlier, the minimum of φpxq is given by φpxq " 0 which is equivalent to Ax ´b " 0. Thus, by minimizing φpxq we are solving the system Ax " b. Since φpx k q " mintφpxq, @x P x 0 `Kt,k u, then φpx k q " φpx k´1 `Qk α k q " mintφpx k´1 `Qk αq, @α P R tk u. (11) Once x k has been chosen, either x k is the exact solution of Ax " b, or t new basis vectors and the new approximation x k`1 " x k `Qk`1 α k`1 are computed. This procedure is repeated until convergence.

Thus, we need to find the recursion relations of r k and α k . By definition, the residual r k " b ´Ax k where x k P x 0 `Kt,k . Thus r k P K t,k`1 . The recursion relation of r k can be simply obtained by replacing x k by its expression as follows:

r k " b ´Ax k " b ´Apx k´1 `Qk α k q " r k´1 ´AQ k α k . At each iteration the step α k`1 is chosen such that φpx k`1 q " mintφpx k `Qk`1 αq, @α P R t pk `1qu. Let F pαq " φpx k `Qk`1 αq where φpxq " 1 2 x t Ax ´xt b. Then, F pαq " 1 2 px k `Qk`1 αq t Apx k `Qk`1 αq ´px k `Qk`1 αq t b " φpx k q `1 2 rpx k q t AQ k`1 α `αt pQ k`1 q t Ax k `αt pQ k`1 q t AQ k`1 αs ´αt pQ k`1 q t b " φpx k q `1 2 rpx k q t AQ k`1 α ´αt pQ k`1 q t Ax k s `1 2 α t pQ k`1 q t AQ k`1 α ´αt pQ k`1 q t r k " φpx k q `1 2 α t pQ k`1 q t AQ k`1 α ´αt pQ k`1 q t r k , since A is SP D.
The minimum of F pαq is given by

F 1 pαq " 0 ñ F 1 pαq " pQ k`1 q t AQ k`1 α ´pQ k`1 q t r k " 0.
Therefore, α k`1 " pQ t k`1 AQ k`1 q ´1pQ t k`1 r k q . By minimizing φpxq, the Petrov-Galerkin condition, r k K K t,k , is ensured (Theorem 3.11). Therefore, pr k q t r i " 0, for all i ă k, and the residuals form an orthogonal set.

Theorem 3.11. The Petrov-Galerkin condition in LRE-CG, r k K K t,k , is equivalent to x k being the minimum of φpxq in x 0 `Kt,k . Proof. 1. x k is the minimum of φpxq in x 0 `Kt,k implies r k K K t,k
The minimum of F pαq " φpx k q " φpx k´1 `Qk αq is given by F 1 pαq " pQ k q t AQ k α ´pQ k q t r k´1 " 0. Since x k is the minimum, then α " α k and

F 1 pαq " ´Qt k r k " 0. Thus r k K K t,k . 2. r k K K t,k implies x k is the minimum of φpxq in x 0 `Kt,k (Proof by contradiction) Assume that r k K K t,k and x k is not the minimum of φpxq in x 0 `Kt,k . Then F 1 pα k q ‰ 0. Hence Q t k r k ‰ 0 and r k is not orthogonal to K t,k
. This contradicts our assumption. Thus x k is the minimum of φpxq.

The LRE-CG Algorithm

After deriving the expressions and the recursion relations of

x k " x k´1 `Qk α k , r k " r k´1 ´AQ k α k , α k " pQ t k AQ k q ´1pQ t k r k´1 q,
we present in Algorithm (2) the LRE-CG algorithm and the performed flops, except for the orthonormalization. We refer to the cost of solving the tk ˆtk linear system from step 5 in Algorithm 2 as Solve α ptkq.

Algorithm 2 The LRE-CG algorithm Flops

Input: A, the n ˆn symmetric positive definite matrix Input: b, the n ˆ1 right-hand side; x 0 , the initial guess or iterate Input: , the stopping tolerance; k max , the maximum allowed iterations Output: x k , the approximate solution of the system Ax " b 1: r 0 " b ´Ax 0 , ρ 0 " ||r 0 || 2 2 , k = 1 2nnz `2n ´1 2: Let W " T pr 0 q, normalise its vectors and then let Q " W 3n 3: while ( ? ρ k´1 ą ||b|| 2 and k ă k max ) do 2n 4:

G " pQ t AQq p2nnz ´nqtk `p2n ´1qt 2 k 2

5:

α " G ´1pQ t rq p2n ´1qtk `Solve α ptkq 6:

x " x `Qα 2tkn 7:

r " r ´AQα 2tkn

8:

ρ k " ||r|| 2 2 2n ´1 9:
Let W " AW p2nnz ´nqt 10:

Orthonormalise the vectors of W against the vectors of Q not included in here 11:

Orthonormalise the vectors of W and let Q " rQ W s not included in here c corresponds to the multiplication AQ and the second term nt 2 k 2 c corresponds to the orthonormalization with respect to previous vectors. As for the memory requirements, we have to store the matrix A and tk c `2 vectors of size nˆ1. And there should be enough memory for the tk c ˆtk c matrix Q t AQ.

However, the multiplication

Q t k AQ k can be reduced since Q k " rQ k´1 W k s. Let Z k " AW k , then D k " AQ k " rAQ k´1 Z k s. At iteration k ´1, D k´1 " AQ k´1 is computed. Thus at iteration k, only Z k " AW k is computed. As for G k " Q t k AQ k , it is equal to G k " Q t k AQ k " Q t k D k " ˆQt k´1 D k´1 Q t k´1 Z k W t k D k´1 W t k Z k ˙" ˆGk´1 Q t k´1 Z k Z t k Q k´1 W t k Z k ˙" ˆGk´1 F k F t k E k ˙,
where

G k´1 is computed at iteration k ´1, F k " Q t k´1 Z k , and E k " W t k Z k . Thus computing G k " Q t k AQ k can be reduced to computing F k and E k .

Algorithm 3 The LRE-CG Algorithm Flops

Input: A, the n ˆn symmetric positive definite matrix Input: b, the n ˆ1 right-hand side; x 0 , the initial guess or iterate Input: , the stopping tolerance; k max , the maximum allowed iterations Output: x k , the approximate solution of the system Ax " b 1: r 0 " b ´Ax 0 , ρ 0 " ||r 0 || 2 2 , k = 1 2nnz `2n ´1 2: Let W " T pr 0 q, normalise its vectors and then let Q " W 3n 3: while ( ? ρ k´1 ą ||b|| 2 and k ă k max ) do 2n

4:

Z " AW p2nnz ´nqt

5:

E " W t Z p2n ´1qt 2 6: if k "" 1 then 7:
D " Z and G " E F " Qp:, 1 : tpk ´1qq t Z p2n ´1qt 2 pk ´1q 10:

G " ˜G F F t E ¸and D " rD Zs 11: end if 12:
α " G ´1pQ t rq p2n ´1qtk `Solve α ptkq 13:

x " x `Qα 2tkn 14:

r " r ´Dα 2tkn

15:

ρ k " ||r|| 2 2 2n ´1 16: Let W " Z 17:
Orthonormalise the columns of W against those of Q not included in here 18:

Orthonormalise the vectors of W and let Q " rQ W s not included in here 19: k = k+1 1 20: end while Then, the cost of k c iterations of LRE-CG using Algorithm 3, except for the orthonormalization, is Total Flops " 2nnz `7n ´1 `kc r2nnzt ´nt `p6n ´1qt kc`1 2 `p2n ´1qt

2 kc`1 2 `2ns `řkc k"1 Solve α ptkq " 2nnz `7n ´1 `kc rp2nnz `2n ´1 2 qt `pn ´1 2 qt 2 `p6n ´1qt kc 2 `p2n ´1qt 2 kc 2 `2ns `řkc k"1 Solve α ptkq « 3ntk 2 c `nt 2 k 2 c `nt 2 k c `řkc k"1
Solve α ptkq Note that tk c should be much smaller than n, or otherwise the cost of Algorithm 3 would be Opn 3 Solve α pnqq, and n vectors of size n have to be stored in addition to the n ˆn system Q t kc AQ kc . One remedy to this problem, that we do not address in this paper, is to restart LRE-CG after some iterations. But this restart might have an effect on convergence as in restarted GMRES. Another alternative is to choose a linearly independent subset of the t computed vectors at each iteration i. This reduces the size of the system solved at each iteration. A third alternative is to compute at each iteration i, t i vectors and then choose a linearly independent set of cardinality p t i , where t 0 " t, t i ď t, p t i ď t i , and p t i " t i`1 . This reduces not only the size of the system solved at each iteration, but also the memory requirements and the number of computed vectors per iteration. In exact precision, the second and third alternatives are equivalent by Theorem 3.4, since if a vector T j pr 0 q is linearly dependent on tT 1 pr 0 q, .., T j´1 pr 0 q, T j`1 pr 0 q, ..., T t pr 0 qu then AT j pr 0 q is linearly dependent on tAT 1 pr 0 q, .., AT j´1 pr 0 q, AT j`1 pr 0 q, ..., AT t pr 0 qu. However, this has to be tested in finite precision. Note that there is an additional cost for choosing a linearly independent subset of the t or t i vectors.

The tk ˆtk α system can be solved using iterative methods like Jacobi method or Krylov subspace methods. Moreover, the s-step or communication avoiding Krylov subspace methods can be used to reduce communication. We use matlab's backslash to solve the α systems in the convergence tests that follow.

Convergence Results

After introducing the new CG methods, MSDO-CG and LRE-CG, we compare their convergence behavior with respect to different A-orthonormalization and orthonormalization schemes respectively. Then we compare the convergence behavior of both methods with repect to CG, Coop-CG, MSD-CG on several matrices for different number of partitions (2, 4, 8, 16, 32 and 64 partitions) or number of initial guesses (for Coop-CG only). The matrices are first reordered using Metis's kway partitioning [START_REF] Karypis | Metis 4.0: Unstructured graph partitioning and sparse matrix ordering system[END_REF] that defines the subdomains δ i . Then x is chosen randomly using matlab's rand function and b " A ˚x, except for the ELASTICITY3Dmatrix where A and b are available. In tables 2, 3, and 4, "Iter" is the number of iterations, k c , needed for convergence and "Err" is the relative error

||x´x kc ||2 ||x||2
at convergence. The first matrix POISSON2Dis a block tridiagonal matrix obtined from Poisson's equation (sparse) using matlab's gallery('poisson',100) function. The matrices refered to as NH2D, SKY2D, SKY3D, and ANI3D, arise from boundary value problem of the convection diffusion equations ηpxqu `divpapxquq ´divpκpxq∇uq " f in Ω

u " 0 on BΩ D Bu Bn " 0 on BΩ N
where Ω " r0, 1s n , (n " 2, or 3) and BΩ N " BΩzBΩ D . The function η, the vector field a, and the tensor κ are the given coefficients of the partial differential operator. In the 2D case, we have BΩ D " r0, 1sˆt0, 1u, and in the 3D case, we have BΩ D " r0, 1s ˆt0, 1u ˆr0, 1s. We focus on the following cases:

• NH2D: A non-homogeneous problem with large jumps in the coefficients.The coefficient η and a are both zero. The tensor κ is isotropic and discontinuous. It jumps from the constant value 10 3 in the ring 1 2 ?

2 ď |x ´c| ď 1 2 , c " p 1 2 , 1 2 q T , to 1 outside.

Inria

• SKY2D and SKY3D Skyscraper problems : The tensor κ is isotropic and discontinuous. The domain contains many zones of high permeability which are isolated from each other κpxq "

" 10 3 ˚pr10 ˚x2 s `1q if r10x i s is odd, i " 1, 2 1,
otherwise. where we note rxs as the integer value of x. SKY2Dand SKY3Dare discretized on a 2D and 3D cartisian grids respetively.

• ANI3D Anisotropic layers : the domain is made of 10 anisotropic layers with jumps of up to four orders of magnitude and an anisotropy ratio of up to 10 3 in each layer. The domain is divided into 10 layers parallel to z " 0, of size 0.1, in which the coefficients are constant. We have κ y " 10κ x and κ z " 1000κ x . The velocity field is zero.

POISSON2D, NH2D and SKY2D are discretized on a 100 ˆ100 2D cartesian grid. SKY3D and ANI3D are discretized on a 20 ˆ20 ˆ20 grid.

As for the ELASTICITY3D matrix, it arises from the linear elasticity problem with Dirichlet and Neumann boundary conditions, defined as follows:

divpσpuqq `f " 0 on Ω, u " u D on BΩ D , σpuq.n " g on BΩ N , where Ω is a 3D 30ˆ10ˆ10 parallelepiped, Ω D is the Dirichlet boundary, Ω N is the Neumann boundary, u is the unknown displacement field, f is some body force, σpuq is the Cauchy stress tensor given by Hooke's law. The ELASTICITY3D matrix was discretized with P1 finite elements and a triangular mesh using FreeFem++ [START_REF] Hecht | New development in freefem++[END_REF]. For a detailed description of the problem refer to [START_REF] Grigori | Robust algebraic Schur complement preconditioners based on low rank corrections[END_REF]. Table 1 briefly describes the test matrices. In table 2 we compare the convergence behvior of the MSDO-CG method (Algorithm 1) with different A-orthonormalization schemes for A-orthonormalizing P k against previous P i 's (MGS, CGS, CGS2) and then A-orthonormalizing P k against itself (MGS, CGS, CGS2, A-CholQR, Pre-CholQR) and for different number of partitions t " 2, 4, 8, 16, 32, 64 that correspond to the maximum number of vectors added at each iteration to the enlarged Krylov subspace, K t,k . We have tested different combinations of A-orthonormalization, but we only show MGS (MGS+MGS), CGS+A-CholQR, CGS+Pre-CholQR, CGS2+A-CholQR, and CGS2+Pre-CholQR. Note that MSDO-CG with CGS Aorthonormalization (CGS+CGS) did not converge neither with CGS2 A-orthonormalization (CGS2+CGS2) nor with CGS2+CGS or CGS+CGS2 A-orthonormalization. The reason is that the seach directions are not A-othogonal to satisfactory precision. And by Theorem 3.10, this implies that r k M K t,k . Thus, nothing guarentees convergence since we have shown in section 3.1.3 that MSDO-CG converges faster than CG if r k K K t,k . Moreover, we did not test combinations of MGS and QR factorizations since MGS is expensive in terms of communication compared to the other methods (section 5.1). But we tested MSDO-CG with MGS for comparison purposes since MGS is known for its numerical stability. As shown in table 2, MSDO-CG with MGS A-orthonormalization converges for all the tested matrices and as we increase t, the number of iterations needed for convergence decreases. As we mentioned earlier, MSDO-CG with CGS A-orthonormalization did not converge. Therefore, we replaced CGS with CGS+A-CholQR and with CGS+Pre-CholQR A-orthonormalization. We notice that MSDO-CG with CGS+A-CholQR A-orthonormalization and MSDO-CG with CGS+Pre-CholQR A-orthonormalization have the same convergence behavior. For the matrices POISSON2D and NH2D, both methods converge with the same number of iterations as MSDO-CG with MGS A-orthonormalization. However, for the matrix SKY2D, both methods did not converge. As for the matrices SKY3D and ANI3D, both methods converged only for t " 2 partitions, and t " 2, 4, 8 partitions respectively. The reason for this difference in behavior for different matrices is the condition number (cond 2 " ||A|| 2 ||A ´1|| 2 ). The condition number of the matrices POISSON2D and NH2D is 6 ˆ10 3 , whereas that of the matrices SKY3D, ANI3D Inria of iterations needed for convergence decreases when increasing the number of partitions t. However, as mentioned in section 5.1, MGS is expensive in terms of communication (Optklogptqq messages per iteration, where t processors A-orthonormalized t vectors against tk vectors). Thus, we tested the LRE-CG method with Classical Gram Schmidt (CGS) orthogonalization which requires sending Optlogptqq messages per iteration. The LRE-CG with CGS converges in the same number of iterations as LRE-CG with MGS for the matrices POISSON2D and NH2D. However, for the other matrices, it does not converge for the given stopping criteria except for t " 2 as shown in table 3. The reason is that the the matrix C " Q t AQ is becoming close to singular, with rankpCq ă tk, as the iterations proceed due to the loss of orthogonality in the CGS orthogonalization. The number of iterations in parentheses in table 3 is not the number of iterations for convergence but it denotes the iteration at which the C matrix becomes close to singular.

In CA-GMRES [START_REF] Mohiyuddin | Minimizing communication in sparse matrix solvers[END_REF], the authors use a parallelizable tall and skinny QR (TSQR) factorization [START_REF] Demmel | Communication-avoiding parallel and sequential qr factorizations[END_REF] for orthonormalizing the n ˆt tall and skinny matrix instead of CGS. They have shown that the combination of CGS for orthonormalizing W against Q and TSQR for orthonormalizing W is stable. We have tested LRE-CG with CGS and TSQR (CGS+TSQR) orthonormalization, and it has the same convergence behavior as LRE-CG with MGS (MGS+MGS) orthonormalization (table 3). Thus, we concude that MGS, and CGS+TSQR orthonormalizations are stable enough to be used in the LRE-CG method from Algorithm 3.

In table 4, we compare the convergence behavior of MSDO-CG with MGS A-orthonormalization, LRE-CG with MGS orthonormalization, Coop-CG and MSD-CG with respect to CG for several matrices with different number of partitions t " 2, 4, 8, 16, 32, 64 . The MSDO-CG, COOP-CG and LRE-CG have better convergence than CG, and LRE-CG has the best convergence. The MSD-CG converges, but requires more iterations than CG, three times more iterations for the matrices SKY2D, SKY3D, ANI3D, and ELASTICITY3D. As for Coop-CG, which starts with t different initial guesses and solves two systems of fixed size t ˆt, its convergence is slightly better than MSDO-CG for the matrices POISSON2D, NH2D, and ELASTICITY3D. But it requires much more iterations than both MSDO-CG and LRE-CG for the other matrices ( SKY2D, SKY3D, ANI3D). Moreover, the results may vary depending on the t initial guesses that are used for the different matrices. For the tested matrices, LRE-CG has slightly better convergence than MSDO-CG, since it uses the whole basis to define the new approximate solution rather than t search directions. For the matrices POISSON2D and NH2D, LRE-CG and MSDO-CG have almost the same convergence as CG for t " 2, and then as t is doubled the iterations needed for convergence is decreased by 20% to 30%. For t " 2, LRE-CG requires 35% and 40% less iterations than CG for the matrices ELASTICITY3D and SKY3D respectively. And as t is doubled the number of iterations needed for convergence is decreased by 25% to 30%, and 32% to 45% respectively. For t " 2, LRE-CG requires 60% and 80% less iterations than CG for the matrices SKY2D and ANI3D respectively. And as t is doubled, the number of iterations needed for convergence is decreased by 45% to 50% and 25% to 40% respectively.

As it is clear from the convergence tests, by doubling t the number of iterations needed for convergence is not always reduced by 50% for all the matrices. However, as shown in the previous sections the memory requirements for MSDO-CG and LRE-CG, except for the matrix A, is ptk c `2qn `ptk c q 2 and ptk c `t `2qn `t2 respectively, where n is the size of the matrix, k c is the number of computed iterations. Thus by doubling t, the memory requirements for MSDO-CG and LRE-CG for performing k iterations is at least doubled. But, when t is doubled k c is decreased. Thus, the memory requirements increase and at most double, when t is doubled as shown in Table 5. Hence, t should be relatively small depending on the size of the matrix, on the performed iterations and on the available memory.

Parallel Model and Expected Performance

In this section we describe the parallelization of the MSDO-CG method (section 5.1) and the LRE-CG method (section 5.2) with computed flops, number of messages and words sent and the estimated parallel runtime. The estimated time for computing z flops is γ c z, where γ c is the inverse floating-point rate, also called the floating-point throughput (seconds per floating-point operation). The estimated time for sending a mesages of size k is α c `βc k, where α c is the latency (with units of seconds) and β c is the inverse bandwidth (seconds per word). Hence, the estimated runtime of an algorithm with a total of z computed flops and s sent mesages each of size k is the sum of their corresponding estimated times γ c z `αc s `βc .

For simplicity, we assume that the algorithms are executed on a distributed memory machine formed by t processors, where t corresponds to the number of vectors computed at each iteration. We partition the graph of A into t subdomains using k-way partitioning or another graph partitioning. We denote by δ i , for i " 1, 2, .., t the subsets of indices obtained from the partitioning. That is δ i X δ h " φ for all i ‰ h, Y t h"1 δ h " t1, 2, 3, ..., nu, and |δ i | « n t . Then each processor i is assigned the n t ˆn rowwise part of the matrix A (Apδ i , :q " Ap:, δ i q since A is SPD), the n t ˆ1 rowwise part of the vector b (bpδ i q), and the vector x 0 pδ i q, where δ i " AdjacentpGpAq, δ i q is the adjacent of δ i in the graph of A. Processor i computes x k pδ i q.

However, for performance reasons and due to the multicore nature of most architectures, it is possible to use a number of processors greater than t, preferably a multiple of t. In this case, we start by partitioning the graph of A into t subdomains using k-way partitioning or another graph partitioning, where δ i for i " 1, 2, .., t are the subsets of indices obtained from the partitioning. This partitioning is used to define the T p.q operator and eventually the enlarged Krylov subspace. Assuming that we have ct processors, then every c processors are assigned an n t ˆn rowwise part of the matrix A, Apδ i , :q, n t ˆ1 rowwise part of the vector b (bpδ i q) and the vector x 0 pδ i q, and should output x k pδ i q. In other words, we partition each of our t subdomains into c non-overlapping subdomains to obtain a total of ct subdomains with set of indices δ i,j , where i " 1, 2, .., t, j " 1, 2, .., c, and δ i " Y c j"1 δ i,j . Then, in Algorithms 4 and 5, logptq is replaced by logpctq, and n t is replaced by n ct .

Inria

MSDO-CG

In this section we describe the parallelization of the MSDO-CG algorithm and we estimate its runtime in terms of flops, number of messages, and words sent. As mentioned in section 4, MGS, CGS2+A-CholQR, and CGS2+Pre-CholQR A-orthonormalizations are numerically the most stable and allow the convergence of MSDO-CG for the matrices in our test set. As discussed in Appendix B, the most parallelizable versions of MGS, Algorithms 14 and 15, require sending ptk `1qlogptq and 2pt ´1qlogptq messages respectively. Whereas CGS2, Algorithm 22, requires sending 4logptq messages. On the other hand, Algorithm 25 of A-CholQR requires sending logptq messages, and Pre-CholQR Algorithm 27 requires sending 3logptq messages. The CGS2+A-CholQR and CGS2+Pre-CholQR A-orthonormalizations can be called communication avoiding since they require sending 5logptq and 7logptq messages respectively, unlike the MGS A-orthonormalization. Since both methods are stable and CGS2+A-CholQR requires less communication, we present the Parallel MSDO-CG with CGS2+A-CholQR A-orthonormalization in Algorithm 4.

Algorithm 4 Parallel MSDO-CG with CGS2+A-CholQR A-orthonormalization Estimated Time Input: δ i , the set of indices assigned to processor i Input: Apδ i , :q, the n t ˆn row part of A Input: bpδ i q, the n t ˆ1 row part of b; x 0 pδ i q, the |δ i | « n t `ci row part of r 0 Input: , the stopping tolerance; k max , the maximum allowed iterations Output: x k pδ i q, the row part of the approximate solution of Ax " b 1: for each processor i " 1 : t in parallel do 2:

Processor i computes rpδ i q " bpδ i q ´Apδ i , δ i qx 0 pδ i q, and let k " 1 γ c p2 nnz t q 3:

Processor i computes rpδ i q t rpδ i q and receives the full ρ " ||r|| 2 2 γ c p2 n t q via an all reduce (overlapped with the next computation)

`pα c `βc qlogptq 4:

Processor i sends P 1 pδ i , :q " rT prqspδ i , :q " r0, .., 0, rpδ i q, 0, .., 0s to its γ c p2 nnz t q m M B neighboring processors and receives from them the corresponding `αc m M B blocks to obtain P 1 pδ i , :q. Then it computes W 1 pδ i q " Apδ i , δ i qP 1 pδ i , :q `βc n t m M B

5:

Call A-CholQR algorithm 25 to A-orthonormalize P 1 γ c 4nt `pα c `βc t 2 qlogptq 6:

while ( ? ρ ą ||b|| 2 and k ă k max ) do 7:

Processor i computes P k pδ i , :q t rpδ i q and receives the full γ c p2 n t ´1qt α " P t k r via an all reduce `pα c `tβ c qlogptq 8:

Processor i computes xpδ i q " xpδ i q `Pk pδ i , :qα 4γ c n t t and rpδ i q " rpδ i q ´Wk pδ i , :qα 9:

Processor i computes rpδ i q t rpδ i q and receives ρ " ||r|| 2 2 γ c p2 n t ´1q via an all reduce (overlapped with the next communication)

`pα c `βc qlogptq 10:

Processor i computes ´Wk pδ i , :q t rpδ i q and receives the full γ c p2 n t ´1qt β " ´W t k r via an all reduce `pα c `tβ c qlogptq 11:

Processor i computes P k`1 pδ i , :q " rT prqspδ i , :q `Pk pδ i , :qdiagpβq 2γ c n t t

12:

Processor i sends rT prqspδ i , :q to its m M B neighboring processors γ c p2 nnz t q and receives from them the corresponding blocks to obtain `αc m M B rT prqspδ i , :q. Then it computes Zpδ i q " Apδ i , δ i qrT prqspδ i , :q `βc n t m M B

13:

Processor i computes W k`1 pδ i , :q " Zpδ i , :q `Wk pδ i , :qdiagpβq 2γ c n t t 14:

Call CGS2 Algorithm 22 to A-orthonormalize P k`1 against P i 12γ c ntk for all i ď k `2p2α c `βc t 2 kqlogptq 15:

Call A-CholQR Algorithm 25 to A-orthonormalize P k`1 γ c 4nt `pα c `βc t 2 qlogptq 16:

k " k `1

17:

end while 18: end for

In Algorithm 4 we have two types of communication. The first is an "all reduce" communication that requires synchronization between all the processors and is equivalent to logptq messages, each of the same size (refer to [START_REF] Thakur | Improving the performance of Collective Operations in MPICH[END_REF]). For example, in line 10 of Algorithm 4, the "all reduce" is equivalent to logptq messages each of size t words, since β is a vector of size t.

The second type of communication is a point-to-point communication between each processor i and its m i neighboring processors for computing a matrix -block of vectors muliplication, specifically ArT prqs. We denote by m M B " maxtm i | i " 1, 2, .., tu the largest number of neighboring processors, where m i ď m M B ď pt ´1q for all i. Note that processor i has to compute Apδ i , δ i qrT prqspδ i , :q where δ i " AdjacentpGpAq, δ i q. Then, the neighboring processors of a given processor i are defined as all the processors j from which processor i needs some rows of rT prqs to compute its part of ArT prqs. In other words, neighboring processors are all the processors j for which δ i X δ j ‰ φ. Moreover, rT prqspδ i , :q is all zeros except for the i th column which is equal to rpδ i q. Thus, processor i sends rpδ i q of size n t ˆ1 to its neighboring processors once rpδ i q is computed at step 8. Since rpδ i q is used in the computation at step 12, this communication is overlapped with the computations from step 9 to 11. Simultaneously, processor i receives rpδ j q from all its neighboring processors j for j " 1, 2, .., m i . Then it computes Apδ i , δ i qrT prqspδ i q by performing approximately 2nnz´n t flops. In summary, without the A-orthogonalization at steps 14 and 15, the estimated time of k c iterations of Algorithm 4, where we ignore lower order terms, is

γ c p11n `2n nnz t qk c `αc p2logptq `mMB qk c `βc p n t m M B `2tlogptqqk c .
At iteration k, the CGS2+A-CholQR A-orthonormalization requires sending 5logptq messages with pt 2tk `2qtlogptq words and performing approximately 12ntk `4nt `6n flops. After k c iterations the estimated time for the A-orthonormalization is γ c p12ntk c `16nt `6nqk c `αc p5logptqqk c `βc pt 2t pkc`1q 2 `2qtlogptqk c . Thus, the estimated time of k c iterations of algorithm 4 is Time MSDO-CG pk c q « γ c p2n nnz t `12ntk c `10nt `17nqk c `αc p7logptq `mMB qk c `βc p n t m M B `t2 k c logptqqk c .

LRE-CG

In this section we describe the parallelization of the LRE-CG algorithm and we estimate its runtime in terms of flops, number of messages, and words sent. As mentioned in section 4, MGS and the CGS+TSQR orthonormalizations are numerically the most stable and allow the convergence of LRE-CG for the matrices in our test set. The parallel version of MGS orthonormalization is similar to that of the Aorthonormalization discussed in Appendix B, Algorithms 14 and 15, and requires sending ptk `1qlogptq and 2pt ´1qlogptq messages respectively. Whereas the CGS orthonormalization can be parallelized in a block format like Algorithm 18, and requires sending 2logptq messages. On the other hand, the TSQR orthonormalization using binary trees as discussed in [START_REF] Demmel | Communication-avoiding parallel and sequential qr factorizations[END_REF] requires sending logptq messages. The combination of BCGS and TSQR was discussed in [START_REF] Mohiyuddin | Minimizing communication in sparse matrix solvers[END_REF] and it requires sending only 3logptq messages as compared to the ptk `2t ´1qlogptq messages of MGS. We present the Parallel LRE-CG with BCGS+TSQR orthonormalization in Algorithm 5.

In Algorithm 5 there are two types of communication, similarly to Algorithm 4. The first type of communication is a point-to-point communication between each processor i and its m i neighboring processors for computing the matrix -block of vectors multiplication Z " AW in line 7. We denote by m M B " maxtm i | i " 1, 2, .., tu the largest number of neighboring processors, where m M B ď pt ´1q for all i. Note that processor i has to compute Apδ i , δ i qW pδ i , :q, where δ i " AdjacentpGpAq, δ i q. Then, the neighboring processors of a given processor i are defined as all the processors j from which processor i needs some rows of W to compute its part of AW . In other words, neighboring processors are all the processors j for which δ i X δ j ‰ φ. Note that for the first iteration, W pδ i , :q " rT prqspδ i , :q is all zeros except for the i th column which is equal to rpδ i q. Thus, processor i sends rpδ i q of size n t ˆ1 to Inria Algorithm 5 Parallel LRE-CG with BCGS+TSQR Algorithm Estimated time

Input: δ i , the set of indices assigned to processor i; Apδ i , :q, the n t ˆn row part of A Input: bpδ i q, the n t ˆ1 row part of b; x 0 pδ i q, the |δ i | ˆ1 row part of r 0 Input: , the stopping tolerance; k max , the maximum allowed iterations Output: x k pδ i q, the row part of the approximate solution of Ax " b 1: for each processor i " 1 : t in parallel do 2:

Processor i computes rpδ i q " bpδ i q ´Apδ i , δ i qx 0 pδ i q γ c p2 nnz t q 3:

Processor i computes rpδ i q t rpδ i q and receives the full ρ " ||r|| 2 2 γ c p2 n t ´1q via an all reduce (overlapped with the next computation)

`pα c `βc qlogptq 4:

Let Qpδ i , :q " W pδ i , :q " rT prqspδ i , :q, and normalise its vectors γ c p2 n t tq 5:

Processor i sends W pδ i , :q " rT prqspδ i , :q to its m i neighboring `αc m M B processors and receives from them the corresponding blocks `βc

n t m M B to obtain W pδ i , :q. Let k " 1 6:
while ( ? ρ k´1 ą ||b|| 2 and k ă k max ) do 7:

Processor i computes Zpδ i , :q " Apδ i , δ i qW pδ i , :q 2γ c nnz t t

8:

Processor i computes W pδ i , :q t Zpδ i , :q and receives γ c 2 n t t 2 E " W t Z via an "all reduce"

`pα c `t2 β c qlogptq 9:

if k "" 1 then Processor i computes Qpδ i , 1 : tpk ´1qq t Zpδ i , :q and γ c 2 n t t 2 k receives F " Qp:, 1 : tpk ´1qq t Z via an "all reduce" `αc logptq 13:

G " ˜G F F t E ¸and D " rD Zs `βc t 2 pk ´1qlogptq Processor i computes Qpδ i , :q t rpδ i q and receives Q t r via an γ c 2 n t tk "all reduce"

`pα c `tkβ c qlogptq 16:

α " G ´1pQ t rq T ime α ptkq 17:

Processor i computes xpδ i q " xpδ i q `Qpδ i , :qα γ c p2tk n t q 18:

Processor i computes rpδ i q " rpδ i q ´Dpδ i , :qα γ c p2tk n t q

19:

Processor i computes rpδ i q t rpδ i q and receives ρ k " ||r|| 2 2 γ c p2 n t ´1q via an "all reduce" (overlapped with the next computation)

`pα c `βc qlogptq 20:

Let W " Z and orthogonalise the columns of W against γ c 4ntk `2α c logptq those of Q using BCGS `t2 kβ c logptq 21:

Orthonormalise the vectors of W using TSQR γ c p2 n t t 2 `2 3 t 3 logptqq and let Q " rQ W s `pα c `βc Processor i sends W pδ i , :q to its m i neighboring processors `αc m M B and receives from them the corresponding blocks to obtain W pδ i , :q `βc end while 25: end for its neighboring processors once rpδ i q is computed. However, for the next iteration the W is no longer sparse, therefore W pδ i , :q of size n t ˆt is sent. The second type of communication is an "all reduce" that requires synchronization between all the processors, and it is equivalent to logptq messages each of the same size (refer to [START_REF] Thakur | Improving the performance of Collective Operations in MPICH[END_REF]). For example, in lines 3 and 19 of Algorithm 5, the "all reduce" is equivalent to logptq messages each of size 1 word.

As mentioned, this communication can be overlapped with the next computation. The reception of E " W t Z, F " Qp:, 1 : tpk ´1qq t Z and Q t r via an "all reduce" in lines 8, 12 and 15 of Algorithm 5 is equivalent to logptq messages each of size t 2 words, logptq messages each of size t 2 pk ´1q words, and logptq messages each of size tk words respectively. However, the three computations are independent. Thus, each processor can compute its part of the three aforementioned computations and then receive the full matrices and vectors via logptq messages each of size ktpt `1q words, assuming that it is possible to send t 2 k words in one message. Another alternative is to compute Qpδ i , 1 : tpk ´1qq t Zpδ i , :q in several steps and overlap the communication with the next computation. The number of steps depends on the machine's architecture and on the values of t and k.

In Algorithm 5, we show the estimated time for each computation and communication, where T ime α ptkq is the estimated time for solving the tk ˆtk α system in line 16. At the k th iteration of Algorithm 5 the total flops, except for the α system, is 2nnz `p6nt ´2t 2 `6n ´tqk `2nt `2 3 t 3 logptq `2n `2 n t ´1. And 4logptq `mMB messages are sent with pp2t 2 `tqk `t2 2 `tqlogptq `nm M B words. Then, by ignoring the lower order terms the estimated time of k c iterations of Algorithm 5, where

t ą 1, k c ą 1, is T ime LRE-CG pk c q « γ c p2nnz `3ntk c `2 3 t 3 logptqqk c `αc p4logptq `mMB qk c `βc rt 2 k c logptq `3 2 t 2 logptq `mMB nsk c `řkc k"1 T ime α ptkq
6 Preconditioned enlarged Krylov subspace methods

After introducing the enlarged Krylov subspace methods and proving, theoretically and numerically, that these methods converge, we describe the preconditioned enlarged Krylov methods. A system Ax " b can be left, right, or split preconditioned. In the case of conjugate gadient methods, the matrix A is symmetric positive definite (SPD). Hence, the preconditioned matrix should also be SPD. For left and right preconditioning, it is not easy to find some matrix M such that M ´1A or AM ´1 is SPD. But assuming that M " LL t , then the split preconditioned matrix L ´1AL ´t is SPD with L ´t " pL t q ´1. Given an n ˆn SPD matrix A, n ˆ1 vector b and some preconditioner M " LL t , then the split preconditioned enlarged Krylov subspace corresponding to the system L ´1AL ´ty " L ´1b with y " L t x and M " LL t , is defined by K t,k pL ´1AL ´t, r 0 q " spantT pr 0 q, L ´1AL ´tT pr 0 q, pL ´1AL ´tq 2 T pr0q, ..., pL ´1AL ´tq k´1 T pr 0 qu " spantT 1 pr 0 q, T 2 pr 0 q, ..., T t pr 0 q, L ´1AL ´tT 1 pr 0 q, L ´1AL ´tT 2 pr 0 q, ..., L ´1AL ´tT t pr 0 q, ..., pL ´1AL ´tq k´1 T 1 pr 0 q, ..., pL ´1AL ´tq k´1 T t pr 0 qu, where r 0 " L ´1pb ´AL ´ty 0 q " L ´1pb ´Ax 0 q, y 0 " L t x 0 , and x 0 is the initial guess. Consequently, the split preconditioned enlarged conjugate gradient methods are defined by the orthogonality condition and the subspace condition associated with preconditioned enlarged Krylov subspace. For example, given a split preconditioned system L ´1AL ´ty " L ´1b with y " L t x and M " LL t , the the enlarged CG Krylov projection methods are defined by y k P y 0 `Kt,k pL ´1AL ´t, r 0 q (the subspace condition), and r k K K t,k pL ´1AL ´t, r 0 q (orthogonality condition), where r k " L ´1pb ´AL ´ty k q " L ´1pb ´Ax k q. Assuming p A " L ´1AL ´t , p b " L ´1b, and p x " y, then all the theorems and properties discussed in section 3.1 are valid for the system p Ap x " p b. Given an SPD matrix M , then the Cholesky factorization M " LL t can be used for split preconditioning the system Ax " b, where the matrix L ´1AL ´t is SPD. As the Cholesky factorization of an n ˆn matrix can be expensive, another alternative is to use block Jacobi preconditioner with Cholesky factorization of the diagonal blocks. A four blocks Jacobi preconditioner with Cholesky decomposition Inria has the following form.

M " ¨M1,1 0 0 0 0 M 2,2 0 0 0 0 M 3,3 0 0 0 0 M 4,4 ‹ ‹ ' " ¨L1,1 0 0 0 0 L 2,2 0 0 0 0 L 3,3 0 0 0 0 L 4,4 ‹ ‹ ' ¨Lt 1,1 0 0 0 0 L t 2,2 0 0 0 0 L t 3,3 0 0 0 0 L t 4,4 ‹ ‹ ' (12)
We present in Algorithm 6 the split preconditioned MSDO-CG with CGS2+ p A-CholQR p A-orthonormalization of the system p Ap x " p b, where p A " L ´1AL ´t , p b " L ´1b, p

x " y, and y " L t x. We omit the W recursion due to numerical errors since W " p AP " L ´1AL ´tP consists of performing backward and forward substitution in addition to the matrix vector multiplication. Thus, we use a version of the CGS2 p A-orthonormalization (Algorithm 23) that computes W " p AP " L ´1AL ´tP and outputs it. As for the p A-CholQR, by assuming that W " p AP " L ´1AL ´tP is computed, then we can use Algorithm 26 with input W " p AP . The additional cost of preconditioning is computing at each iteration k, four times W k`1 " L ´1AL ´tP k`1 which is equivalent to a backward and forward substitution with t right hand sides and a matrix vector multiplication. The difference between the preconditioned and unpreconditioned MSDO-CG is in the A-orthonormalization, the computation of W , and the backward substitution L t x " y. Note that the split preconditioned MSDO-CG with MGS p A-orthonormalization did not converge. This might be due to numerical errors in solving the tk backward and forward substitutions in the MGS p A-orthonormalization. However, the MSDO-CG with CGS2+ p A-CholQR p A-orthonormalization converges very well as shown in section 6.1.

Algorithm 6 Split preconditioned MSDO-CG algorithm with CGS2+ p A-CholQR p A-orthonormalization

Input: A, the n ˆn symmetric positive definite matrix; L, n ˆn split preconditioner Input: b, the n ˆ1 right-hand side; x 0 , the initial guess or iterate Input: , the stopping tolerance; k max , the maximum allowed iterations Output: x, the approximate solution of Ax " b by solving for L ´1AL ´ty " L ´1b and y " L t x 1: y " L t x 0 , r " L ´1pb ´Ax 0 q , ρ " ||r|| 2 2 , nb " ||L ´1b|| 2 , k " 1 2: Let P 1 " T prq and L ´1AL ´t-orthonormalize it using Algorithm 26 which outputs W 1 " L ´1AL ´tP 1 3: while ( ? ρ ą pnbq and k ă k max ) do 4:

α " pP t k W k q ´1pP t k rq " P t k r 5:

y " y `Pk α and r " r ´Wk α 6:

ρ " ||r|| 2 2 7:

β " ´pP t k W k q ´1pW t k rq " ´W t k r 8:

P k`1 " T prq `Pk diagpβq 9:
L ´1AL ´t-orthonormalize P k`1 against all P i 's for i ď k using Algorithm 23 10:

L ´1AL ´t-orthonormalize P k`1 using Algorithm 26 which outputs W k`1 " L ´1AL ´tP k`1 11:
k " k `1 12: end while 13: Solve L t x " y In Algorithm 7, we present the split preconditioned LRE-CG algorithm with BCGS+TSQR orthonormalization of the system p Ap x " p b, where p A " L ´1AL ´t , p b " L ´1b, p

x " y, and y " L t x. At first glance, it might appear to the reader that the additional cost at iteration k in Algorithm 7 is solving a forward and backward substitution with tk right hand sides ( L ´1AL ´tQ) and a forward and backward substitution with t right hand sides ( L ´1AL ´tW ). However, by taking a quick look at Algorithm 8, it is clear that the additional cost of preconditioning at iteration k is solving only a forward and backward substitution

Algorithm 7 Split preconditioned LRE-CG Algorithm with BCGS+TSQR orthonormalization

Input: A, the n ˆn symmetric positive definite matrix; L, n ˆn split preconditioner Input: b, the n ˆ1 right-hand side; x 0 , the initial guess or iterate Input: , the stopping tolerance; k max , the maximum allowed iterations Output: x, the approximate solution of Ax " b by solving for L ´1AL ´ty " L ´1b and y " L t x 1: y " L t x 0 , r " L ´1pb ´Ax 0 q, ρ 0 " ||r|| 2 2 , nb " ||L ´1b|| 2 ,k = 1 2: Let W " T pr 0 q, normalise its vectors and then let Q " W 3: while ( ? ρ k´1 ą pnbq and k ă k max ) do 4:

α " pQ t L ´1AL ´tQq ´1pQ t rq 5:

y " y `Qα

6:
r " r ´L´1 AL ´tQα and ρ k " ||r|| 2 2 7:

Let W " L ´1AL ´tW 8:

Orthonormalise the vectors of W against the vectors of Q

9:
Orthonormalise the vectors of W and let Q " rQ W s and k " k `1 10: end while 11: Solve L t x " y Algorithm 8 Split preconditioned LRE-CG Pseudo Code with BCGS+TSQR orthonormalization Input: A, the n ˆn symmetric positive definite matrix Input: b, the n ˆ1 right-hand side; x 0 , the initial guess or iterate Input: , the stopping tolerance; k max , the maximum allowed iterations Output: x k , the approximate solution of the system Ax " b 1: y " L t x 0 , r " L ´1pb ´Ax 0 q, ρ 0 " ||r|| 2 2 , nb " ||L ´1b|| 2 ,k = 1 2: Let W " T pr 0 q, normalise its vectors and then let Q " W 3: while ( ? ρ k´1 ą pnbq and k ă k max ) do 4:

Z " L ´1AL ´tW , E " W t Z 5:

if k "" 1 then 6:

D " Z and G " E α " G ´1pQ t rq 11:

y " y `Qα

12:
r " r ´Dα and ρ k " ||r|| 2 2 13:

Let W " Z 14:

Orthonormalise the columns of W against those of Q using BCGS 15:

Orthonormalise the vectors of W using TSQR and let Q " rQ W s and k " k `1 16: end while 17: Solve L t x " y with t right hand sides ( L ´1AL ´tW q. And this is the only difference with the unpreconditioned version in addition to the backward substitution L t x " y at the end.

Inria

Cholesky decomposition to obtain a block L. Then the matrix A is partitioned once again using kway into t " 2, 4, 8, 16, 32 or 64 parts that define the enlarged Krylov subspace. The preconditioner is permuted accordingly. The preconditioned LRE-CG converges faster than the preconditioned MSDO-CG and PCG. As the number of partitions or the maximum basis vectors added at each iteration is doubled, the preconditioned LRE-CG and MSDO-CG converge faster. However, for some matrices, like POISSON2D, NH2D, and ANI3D, as the number of partitions is doubled, the number of iterations till convergence decreases by only 10% ´20%. Thus, it is efficent to use a maximum of t " 4 partitions which correspond to adding at most t vectors to the basis of the enlarged Krylov subspace. For the matrices SKY2D and SKY3D, the number of iterations decreases by 33% ´45% and 25% ´45% respectively. Hence it is possible to use a maximum of t " 8 partitions or at most t " 16.

Conclusion

In this paper we have introduced two new iterative methods, MSDO-CG and LRE-CG, which are based on the enlarged Krylov subspace. After introducing the related existing methods (Block-CG, Coop-Cg and MCD-CG), we have defined the properties of the enlarged Krylov subspace, derived the new methods in the context of projection CG versions, provided parallel versions that reduce communication, and shown that the methods converge at least as fast as Classical CG in exact precision arithmetic. The convergence results show that they also converge faster than CG in finite precision arithmetic. We have also presented the preconditioned versions and tested their convergence with block Jacobi preconditioner. MSDO-CG is a variation of the MSD-CG version, where we A-orthonormalize the t search directions against previous directions and against each others. Due to the A-orthonormalization, we lose the short recurrence property of CG and we are obliged to save all the tk c search directions, where k c is the number of iterations till convergence. In LRE-CG we start by building an orthonormal basis for the enlarged Krylov subspace, then we use the whole basis to update the solution. The main difference between both methods in terms of performance, is that at each iteration of MSDO-CG, we use t search directions to update the new approximate solution. Whereas in LRE-CG, at each iteration i, we use the entire basis formed by ti vectors, to update the approximate soltion and we solve a ti ˆti system. However, this use of the whole basis leads to a relatively faster convergence than MSDO-CG. One way to limit this increasing cost is by restarting LRE-CG after some iterations. Another alternative is to choose at each iteration i, a linearly independant subset of the t computed vectors. This adds an extra cost, but reduces the size of the system that has to be solved at each iteration. A third alternative is to compute t i vectors at each iteration i, where t 0 " t, and t i ď t. Then choose p t i linearly independant vectors where p t i ď t i , and t i`1 " p t i .

Although each iteration of the MSDO-CG and LRE-CG methods is at least t times more expensive than the CG iteration in terms of flops, as shown in section 5, both methods use less communication, and Blas2 and Blas3 operations that can be parallelized in a more efficent way than the dot products in CG. Moreover, the MSDO-CG and LRE-CG methods converge faster than CG in terms of iterations as shown in section 4.

Our future work will focus on testing the LRE-CG versions discussed above, that are less expensive in terms of flops and memory requirements than LRE-CG, like restarted LRE-CG or LRE-CG with selected basis vectors. Then, the most stable version will be implemented in a parallel environment. We will also test LRE-CG on other real applications' matrices, and with different preconditioners. Moreover, we would also like to compare the runtime of the LRE-CG version with the MSDO-CG method on a parallel environment. We will also derive and test other enlarged Krylov methods, like enlarged GMRES which has been derived but not tested yet. only logptq messages. We also present the Pre-CholQR version that was introduced in [START_REF] Lowery | Stability Analysis of QR factorization in an Oblique Inner Product[END_REF] and requires 3logptq messages.

B.1 Modified Gram Schmidt A-orthonormalization

We start by introducing A-orthonormalization of the vectors of P k`1 against the vectors of all the previous P i 's for i ă k `1, then against each others in the first section. In the second section,we discuss versions that save flops and reduce communication. Finally, in the last section the parallelization of both kernels is described.

B.1.1 The A-orthonormalization using MGS

Assuming that the vectors of P i are A-normalized, i.e. P i p:, jq t AP i p:, jq " 1 for all j " 1, 2, .., t and i " 1, 2, .., k, then the A-orthonormalization of the vectors of P k`1 against the vectors of all the previous P i 's for i ă k `1 is defined as follows:

Algorithm 10 A-orthonormalization against previous vectors with MGS Input: A, the n ˆn symmetric positive definite matrix Input: P 1 , P 2 ,.. , P k`1 , the k `1 sets of search directions Output: P k`1 , the search directions A-orthonomalized against P 1 , P 2 ,.. , P k 1: for o " 1 : t do %loop over the vectors of P k`1

2:

for i " 1 : k do %loop over the different Pi's 3:

for j " 1 : t do %loop over the vectors of Pi

4:

P k`1 p:, oq " P k`1 p:, oq ´pP i p:, jq t AP k`1 p:, oqqP i p:, jq At each inner iteration, one matrix-vector multiplication has to be computed (AP k`1 p:, oq), 1 dot product, and 1 saxpy, which costs 2nnz ´n `p2n ´1q `2n " 2nnz `3n ´1 flops. Then, at each outermost iteration, one matrix-vector multiplication is computed (AP k`1 p:, oq), 1 dot product, 1 square root and 1 division which costs 2nnz ´n `p2n ´1q `2 " 2nnz `n `1. The total cost of Algorithm 10 is p2nnz `3n ´1qt 2 k `p2nnz `n `1qt, which is of the order of nnzt 2 k `nt 2 k.

As for the A-orthonormalization of the vectors of P k`1 against each others, it is defined as follows:

Algorithm 11 A-orthonormalization against each others using MGS Input: A, the n ˆn symmetric positive definite matrix Input: P k`1 , the search directions to be A-orthonormalized Output: P k`1 , the A-orthonomalized search directions 1: for i " 1 : t do %loop over the vectors of P k`1 2:

for j " 1 : pi ´1q do%A-orthogonalize against the vectors P k`1 p:, 1 : i ´1q

3:

P k`1 p:, iq " P k`1 p:, iq ´pP k`1 p:, jq t AP k`1 p:, iqqP k`1 p:, jq P k`1 p:, iq " P k`1 p:,iq ||P k`1 p:,iq|| A " P k`1 p:,iq P k`1 p:,iq t AP k`1 p:,iq %A-normalize 6: end for Similarly, the cost of the inner loop is 2nnz `3n ´1 flops and that of the outer loop is 2nnz `n `1, but the total cost is p2nnz `3n ´1qpt ´1q t 2 `p2nnz `n `1qt flops, which is of the order of nnzt 2 `nt 2 .

B.1.2 Saving flops in the A-orthonormalization using MGS

Since the A-orthonormalizations are expensive in term of flops, we present another alternative for computing the A-orthonormalizations that reduces the computed flops at the expense of storing more vectors.

In Algorithm [START_REF] Grigori | Communication Avoiding ILU0 Preconditioner[END_REF] and Algorithm 10, some matrix vector multiplications are repeatedly computed. For example, in Algorithm 11 AP k`1 p:, 1q is computed t ´1 times, AP k`1 p:, 2q is computed t ´2 times and generally, AP k`1 p:, iq is computed t ´i times, which means that the matrix A is accessed pt ´1q t 2 times for every call of the algorithm. Thus, it is possible after A-orthogonalizing a vector P k`1 p:, iq to compute and store w i " AP k`1 p:, iq. This eliminates the redundant flops and reduces the number of accesses of A to t times, but there is a need to store t extra vectors (w i ).

Moreover, it is possible to further reduce the computations and the number of times A is accessed at the expense of storing tk vectors as shown in Algorithm 12, where the multiplication W k`1 " AP k`1 is first performed by only reading the matrix A once. Then the vectors W k`1 p:, iq are updated and stored.

The A-orthonormalization against previous vectors with flops reduction can be performed as follows: for j " 1 : pi ´1q do%A-orthogonalize against the vectors P k`1 p:, 1 : i ´1q

3:

P k`1 p:, iq " P k`1 p:, iq ´pW k`1 p:, jq t P k`1 p:, iqqP k`1 p:, jq 4n ´1 for j " 1 : t do %loop over the vectors of Pi

3:

P k`1 " P k`1 ´rP i p:, jqs t diagpW i p:, jq t P k`1 q p4n ´1qt 4:

W k`1 " W k`1 ´rW i p:, jqs t diagpW i p:, jq t P k`1 q 2nt 5:

end for 6: end for 7: for o " 1 : t do 8: pap k`1 poq " W k`1 p:, oq t P k`1 p:, oq 2n ´1 9: end for 10: pap k`1 " p ? pap k`1 q t 11: P k`1 " P k`1 diagppap k`1 q ´1 and W k`1 " W k`1 diagppap k`1 q ´1 p2n `2qt

In Algorithm 13, rather than A-orthonormalizing each vector P k`1 p:, iq against all previous vectors P k`1 p:, jq, we can A-orthogonalize P k`1 p:, i`1 : tq against the A-normalised vector P k`1 p:, iq as shown in Algorithm 15. Let rP k`1 p:, jqs t´i be an n ˆpt ´iq block containing t ´i duplicates of the vector P k`1 p:, jq. Then P k`1 p:, i`1 : tq " P k`1 p:, i`1 : tq´rP k`1 p:, jqs t´i diagpW k`1 p:, iq t P k`1 p:, i`1 : tqq Algorithm 15 A-orthonormalization against each others with MGS Input: P k`1 , the search directions to be A-orthonormalized Input: W k`1 , AP k`1 Output: P k`1 , the A-orthonomalized search directions Output: W k`1 , AP k`1 where P k`1 is the A-orthonomalized search directions 1: for i " 1 : pt ´1q do %A-orthogonalize against the vectors P k`1 p:, 1 : i ´1q

2:

P k`1 p:, i `1 : tq " P k`1 p:, i `1 : tq ´rP k`1 p:, jqs t´i diagpW k`1 p:, iq t P k`1 p:, i `1 : tqq 3:

W k`1 p:, i `1 : tq " W k`1 p:, i `1 : tq ´rW k`1 p:, jqs t´i diagpW k`1 p:, iq t P k`1 p:, i `1 : tqq 4:

pap k`1 " W k`1 p:, i `1q t P k`1 p:, i `1q

5:

P k`1 p:, i `1q " P k`1 p:,i`1q ? pap k`1 and W k`1 p:, i `1q " W k`1 p:,i`1q ? pap k`1 6: end for Then the parallelization of Algorithms 14 and 15 goes as follows. We assume that we have t proces-RR n°8597 sors with distributed memory, and each processor pi is assigned a rowwise part of all W j (W j pδ pi , :q) for j " 1, 2, .., k `1 and the same rowwise part of all P j (P j pδ pi , :q) for j " 1, 2, .., k `1 where δ pi Xδ h " φ for all pi ‰ h and Y t h"1 δ h " t1, 2, 3, ..., nu. At each inner iteration of Algorithm 14, each processor pi has to compute P k`1 pδ pi , :q " P k`1 pδ pi , : q ´rP i pδ pi , jqs t diagpW i p:, jq t P k`1 q. First, each processor pi computes a part of the matrix vector multiplication W i pδ pi , jq t P k`1 pδ pi , :q. Then, a communication of the form "all reduce" is performed to send the 1 ˆt W i p:, jq t P k`1 's value to all the processors. Finally, processor pi computes P k`1 pδ pi , :q and W k`1 pδ pi , :q.

Finally, each processor pi computes its corresponding part of the dot product W k`1 pδ i , oq t P k`1 pδ i , oq for all o " 1, 2, .., t and an "all reduce" is used to send pap k`1 's value to all the processors. Then, each processor A-normalizes P k`1 pδ pi , oq and W k`1 pδ pi , oq. All the communication in Algorithm 14 is of the form "all reduce" of a t ˆ1 vector which is equivalent to sending logptq messages and tlogptq words. So, in total ptk `1qlogptq messages and ptk `1qtlogptq words are sent in Algorithm 12. Hence, by ignoring lower order terms we obtain T ime M GS1 Aort « γ c 6nt 2 k `αc tklogptq `βc t 2 klogptq As for the parallelization of algorithm 15, it is similar to that of algorithm 14 where at each inner iteration processor pi computes a part of the matrix vector multiplication W k`1 pδ pi , iq t P k`1 pδ pi , i`1 : tq and then receives the whole 1 ˆt vector, W k`1 p:, iq t P k`1 p:, i `1 : tq, using an "all reduce". Then, it computes P k`1 pδ pi , :q, W k`1 pδ pi , :q and a part of the dot product pap k`1 and receives the whole dot product by an "all reduce". Finally, each processor A-normalizes its part of P k`1 p:, iq and W k`1 p:, iq. Thus, at each iteration 2 "all reduce" communications are performed, where t words are sent in the first and one word in second. So, in total 2pt ´1qlogptq messages are sent in Algorithm 15 where pt ´1qpt `1qlogptq words are sent. Hence, by ignoring lower order terms we obtain

T ime M GS2 Aort « γ c 3nt 2 `αc 2tlogptq `βc t 2 logptq

B.2 Classical Gram Schmidt A-orthonormalization

Since the MGS A-orthonormalization is costly in terms of communication, we introduce the classical Gram Schmidt (CGS) A-orthonormalization and show that it is equivalent to a QR decomposition with A inner product rather than the usual L2 inner product. Then we present the parallelization of the introduced algorithms. In section B.2.1, the A-orthonormalization against previous vectors using CGS is discussed, whereas in section B.2.2 we discuss the A-orthonormalization of the vectors using CGS. Then in section B.2.3 we introduce the CGS A-orthonormalization with reorthogonalization.

B.2.1 A-orthonormalization against previous vectors using CGS

The A-orthonormalization of P k`1 against the vectors of all the previous P i 's for i ă k `1 is defined as in Algorithm 16. More precisely, r P k`1 p:, oq " P k`1 p:, oq ´řk i"1 ř t j"1 pP i p:, jq t AP k`1 p:, oqqP i p:, jq " P k`1 p:, oq ´řk

i"1 P i pP t i AP k`1 p:, oqq If we let W k`1 " AP k`1 , then r P k`1 p:, oq " P k`1 p:, oq ´řk i"1 P i pP t i W k`1 p:, oqq. Moreover, r P k`1 " P k`1 ´řk i"1 P i pP t i W k`1 q. Let Q k " rP 1 , P 2 , ..., P k s, then r P k`1 " P k`1 ´Qk pQ t k W k`1 q.

This represents a Block classical gram schmidt (BCGS) version of the A-orthonormalization.

The total flops performed in Algorithm 17 is 2p2nnz ´nqt `p2n ´1qt 2 k `2t 2 kn `3nt " 4nnzt ǹt `r4nt 2 ´t2 sk « 4nnzt `4nt 2 k. As for the parallelization of Algorithm 17, it is straightforward due to the block format. Assuming that we have t processors with distributed memory, and each processor pi is assigned a rowwise part of A (Apδ pi , :q), a rowwise part of Q k (Q k pδ pi , :q) and a rowwise part of P k`1 , where δ pi X δ h " φ for all pi ‰ h and Y t h"1 δ h " t1, 2, 3, ..., nu. First, processor pi computes Ap:, δ pi qP k`1 pδ pi , :q and receives the full n ˆt matrix W k`1 via an "all reduce". Then it computes Q k pδ pi , :q t W k`1 pδ pi , :q and obtains the full tk ˆt matrix Q t k W k`1 using an "all reduce". Then, processor pi computes r P k`1 pδ pi , :q " P k`1 pδ pi , :q ´Qk pδ pi , :qpQ t k W k`1 q. Another "all reduce" is needed so that processor pi has the full r P k`1 needed to compute Ă W k`1 pδ pi , :q " Apδ pi , :q r P k`1 . Processor pi computes t partial dot products of the form r P k`1 pδ pi , oq t W k`1 pδ pi , oq and obtains the full dot products via an all reduce. Finally each processor A-normalizes its part of P k`1 , i.e r P k`1 pδ pi , iq " r P k`1 pδpi,iq r P k`1 p:,iq t W k`1 p:,iq for all i " 1, 2, .., t. So in total there is a need to perform 4 all reduce for parallelizing algorithm 17.

1: W k`1 " AP k`1 p2nnz ´nqt 2: r P k`1 " P k`1 ´Qk pQ t k W k`1 q p2n ´1qt 2 k `p2tk ´1qnt `nt 3: Ă W k`1 " A r P k`1 p2nnz ´nqt 4: for i " 1 : t do
It is possible to reduce the communication to only two by assuming that W k`1 " AP k`1 has already been computed and it is an input to Algorithm 18 along with W k " AQ k " rW 1 , W 2 , ..., W k s. The only communication is an "all reduce" of the tk ˆt matrix Q t k W k`1 and another "all reduce" of the vector of size t containing the norms of the columns of r P k`1 . We assume that it is possible to send a message of size t 2 k words at once. Thus, 2logptq messages are sent with ptk `1qtlogptq words where

; Ă W k`1 " A r P k`1 1: r P k`1 " P k`1 ´Qk pQ t k W k`1 q p2n ´1qt 2 k `p2tk ´1qnt `nt 2: Ă W k`1 " W k`1 ´Wk pQ t k W k`1 q 2nt 2 k 3: for i " 1 :

7: end for

The CGS A-orthonormalization can be reformulated as a QR factorization

P k`1 " r P k`1 R
where r P k`1 is an A-orthonormal matrix, and R is a t ˆt upper triangular matrix defined by the entries r j,i for all j " 1, 2, .., i and i " 1, 2, .., t. R "

¨r1,1 r 1,2 r 1,3 ¨¨¨r 1,t r 2,2 r 2,3 ¨¨¨r 2,t r 3,3 ¨¨¨r 3,t . . . . . . r t,t ‹ ‹ ‹ ‹ ‹ ' " ¨||p 1 || A ă p1 , p 2 ą A ă p1 , p 3 ą A ¨¨¨ă p1 , p t ą A r 2,2 ă p2 , p 3 ą A ¨¨¨ă p2 , p t ą A r 3,3 ¨¨¨ă p3 , p t ą A . . . . . . r t,t ‹ ‹ ‹ ‹ ‹ '
Although the CGS A-orthonormalization is equivalent to a QR factorization with the A inner poduct, we were not able parallelize it using reduction trees with the same communication pattern as in TSQR [START_REF] Demmel | Communication-avoiding parallel and sequential qr factorizations[END_REF]. But we can optimize the communication in algorithm 19 by noticing that once a vector p i is orthonormalized, we can compute the corresponding entries of the matrix R, i.e Rpi, i `1 : tq. By taking this into consideration, algorithm 19 can be restructured as in algorithm 20.

Inria

Algorithm 20 QR factorization with A inner product using CGS Flops Input: A, the n ˆn symmetric positive definite matrix Input: P k`1 , the search directions to be A-orthonormalized Output: r P k`1 , the A-orthonomalized search directions r P t k`1 A r P k`1 " I Output: R, the upper triangular matrix such that P k`1 " r The parallelization of Algorithm 20 starts by distributing the data similarly to Algorithm 17. Processor pi computes Ap:, δ pi qP k`1 pδ pi , :q and receives W k`1 via an "all reduce", and computes P k`1 pδ pi , 1q t W pδ pi , 1q, and receives the full dot product P k`1 p:, 1q t W p:, 1q, needed to compute Rp1, 1q, via an "all reduce". Then, it computes r P k`1 pδ pi , 1q "

P k`1 R 1: W k`1 " AP k`1 p2nnz
P k`1 pδpi,1q
Rp1,1q . At each iteration , processor pi computes r P k`1 pδ pi , i´1q t W k`1 pδ pi , i : tq and receives the full Rpi1 , i : tq by an all reduce. Then, it computes r P k`1 pδ pi , iq " P k`1 pδ pi , iq´r P k`1 pδ pi , 1 : i´1qRp1 : i´1, iq and recieves the P k`1 pβ pi , iq from m M B adjacent processors where β pi " AdjacentpGpAq, δ pi q. Then it computes Ă W k`1 pδ pi , iq " Apδ pi , β pi q r P k`1 pβ pi , iq and r P k`1 pδ pi , iq t Ă W k`1 pδ pi , iq and receives the full r P k`1 p:, iq t A r P k`1 p:, iq, needed to compute Rpi, iq, via an all reduce. Finally, it computes r P k`1 pδ pi , iq "

r P k`1 pδpi,iq Rpi,iq
. So, there is a need for a total of 2t ´1 "all reduce" and t communications with the m M B neighboring processors.

It is possible to reduce the communication to only 2t ´1 "all reduce" by assuming that W k`1 " AP k`1 has already been computed and it is an input to Algorithm 21. Then, at each iteration i, an "all reduce" of the vector Rpi ´1, i : tq of size t ´i `1 is performed and another "all reduce" of the entry Rpi, iq is performed. Thus, a total of p2t ´1qlogptq messages are sent with p1 `řt i"2 t `2 ´iqlogptq " The CGS with reorthogonalization (CGS) consists of calling the CGS algorithms twice, be it for Aorthonormalizing P k`1 against previous vectors of Q k (Algorithm 22), or A-orthonormalizing P k`1 .

Algorithm 22 A-orthonormalization of P k`1 against previous vectors of Q k using CGS2 Input: Q k , the tk search directions Input: P k`1 , the t search directions to be A-orthonormalized Input: W k`1 , " AP k`1 ; W k , " AQ k Output: r P k`1 , the search directions A-orthonomalized against P 1 , P 2 ,.. , P k ; Ă W k`1 " A r P k`1 1: Call Algorithm 18 with P k`1 and W k`1 as input and with P In the case of L ´1AL ´t-orthonormalization of P k`1 against previous vectors of Q k where L ´t " pL t q ´1, the CGS2 algorithm is defined in Algorithm 23. Note that we have to solve 6 systems with multiple right hand sides. If L is a lower triangular matrix, then we perform three backward substitutions and three forward substitutions.

B.3 Cholesky QR A-orthonormalization

A-orthonormalizing the n ˆt full rank matrix P k`1 is equivalent to a QR factorization P k`1 " r P k`1 R where r P t k`1 A r P k`1 " I. Thus, P t k`1 AP k`1 " p r P k`1 Rq t A r P k`1 R " R t R and R can be obtained by performing a Cholesky factorization of the SPD matrix P t k`1 AP k`1 . Then, r P k`1 " P k`1 R ´1 is obtained by solving the lower triangular system R t r P t k`1 " P t k`1 with multiple right-hand sides. This procedure is called A-CholQR and summarized in Algorithm 24 [START_REF] Rozloznik | Numerical stability of orthogonalization with a non-standard inner product[END_REF][START_REF] Lowery | Stability Analysis of QR factorization in an Oblique Inner Product[END_REF]. Similarly to the other A-orthonormalization methods, we may assume that W k`1 is already computed, then the obtained A-CholQR is described in Algorithm 25.

The parallelization of Algorithm 25 assumes that we have t processors and each is assigned a rowwise part of P k`1 and W k`1 corresponding to the δ i subset of indices defined previously, P k`1 pδ i , :q and W k`1 pδ i , :q. And each processor i should compute r P k`1 pδ i , :q and Ă W k`1 pδ i , :q. Then each processor i computes W k`1 pδ i , :q t P k`1 pδ i , :q and receives the t ˆt matrix C via an "all reduce" or equivalently Inria Algorithm 23 L ´1AL ´t-orthonormalization against previous vectors of Q k with CGS2 Input: A, the n ˆn symmetric positive definite matrix; L, n ˆn preconditioner Input: Q k , the tk search directions Input: P k`1 , the t search directions to be L ´1ApL t q ´1-orthonormalized Output: p P k`1 , the search directions L ´1ApL t q ´1-orthonomalized against Q k Output: x W k`1 " L ´1AL ´t p P k`1 1: W k`1 " L ´1AL ´tP k`1 2: r P k`1 " P k`1 ´Qk pQ t k W k`1 q 3: Ă W k`1 " L ´1AL ´t r P k`1 4: for i " 1 : t do %loop over the vectors of P k`1 and L ´1 AL ´t-normalize logptq messages and t 2 logptq words. Finally, each processor i can compute the Cholesky factorization of the matrix C to obtain R which is needed to solve R t r P k`1 pδ i , :q t " P k`1 pδ i , :q t and R t Ă W k`1 pδ i , : q t " W k`1 pδ i , :q t . Thus, it is possible to parallelize the A-CholQR A-orthonormalization, Algorithm 25, by sending logptq messages with t 2 logptq words and performing t 2 `p4n´1qt 2 t « 4nt flops in parallel. Hence, by ignoring lower order terms we obtain T ime A´CholQR « γ c 4nt `αc logptq `βc t 2 logptq

In the case of p A-orthonormalization of P k`1 , where p A " L ´1AL ´t and L ´t " pL t q ´1, the p A-CholQR algorithm is defined in Algorithm 26. Note that we have to solve 2 systems with multiple right hand sides. If L is a lower triangular matrix, then we perform a backward and forward substitution. [START_REF] Lowery | Stability Analysis of QR factorization in an Oblique Inner Product[END_REF], which they call Pre-CholQR (Algorithm 27). It consists in performing a Euclidean QR factorization with L2 before calling the A-CholQR A-orthonormalization, Algorithm 24. The QR factorization of P k`1 can be done using the TSQR [START_REF] Demmel | Communication-avoiding parallel and sequential qr factorizations[END_REF], which requires sending logptq messages, each of size t 2 2 words and computing 2nt 2 3 t 3 logptq. Then, parallelizing Algorithm 24 requires performing two "all reduce" or 2logptq messages with pnt `t2 qlogptq words. In total, parallelizing Algorithm 27 requires sending 3logptq messages with pnt `1.5t

Algorithm

D

  " rD Zs, F " Qp:, 1 : tpk ´1qq t Z, and G "

  1 p:, oq " P k`1 p:,oq ||P k`1 p:,oq|| A " P k`1 p:,oq ? P k`1 p:,oq t AP k`1 p:,oq %A-normalize 8: end for

  %loop over the vectors of P k`1 and A-normalize 5: r P k`1 p:, iq " r P k`1 p:,iq || r P k`1 p:,iq|| A " r P k`1 p:,iq ? r P k`1 p:,iq t W k`1 p:,iq 3n 6: end for

1 k`1 and W 1 k`1 as output 2 : 1 k`1 and W 1 k` 1

 112111 Call Algorithm 18 with P as input and with r P k`1 and Ă W k`1 as output

  The cost of the LRE-CG, using Algorithm 2, except for the orthonormalization in steps 10 and 11, is Solve α ptkq, where the first term nnztk2

	Total Flops " 2nnz `7n ´1 `kc rp2nnz `5n ´1qt kc`1 2 `2n `2nnzt ´nts `řkc k"1 Solve α ptkq `p2n ´1q t 2 6 pk c `1qp2k c `1q " 2nnz `7n ´1 `kc rp2nnz `5n ´1qt kc`1 2 `2n `2nnzt ´nts `p2n ´1q t 2 6 p2k 2 c `3k c `1q " 2nnz `7n ´1 `p2n ´1q t 2 6 `kc rp2nnz `5n ´1qt kc`1 2 `2n `2nnzt ´nt `řkc k"1 Solve α ptkq `p2n ´1q t 2 6 p2k c `3qs `řkc k"1 Solve α ptkq
		« nnztk 2 c `nt 2 k 2 c	`řkc k"1
	12:	k = k+1	1
	13: end while
			Inria

Table 1 :

 1 The test matrices

	Matrix	Size	Nonzeros	Symetric 2D/3D	Problem
	POISSON2D	10000	49600	Yes	2D	Poisson equations
	NH2D	10000	49600	Yes	2D	Boundary value
	SKY2D	10000	49600	Yes	2D	Boundary value
	SKY3D	8000	53600	Yes	3D	Skyscraper
	ANI3D	8000	53600	Yes	3D	Anisotropic Layers
	ELASTICITY3D 11253	373647	Yes	3D	Linear Elasticity P1 FE

Table 2 :

 2 Comparison of the convergence of MSDO-CG with different A-orthonormalization schemes, with respect to number of partions (Pa) with x 0 " 0.

				MSDO-CG with different A-Orthonormalization Methods
		MGS	CGS+A-CholQR CGS+Pre-CholQR CGS2+A-CholQR CGS2+Pre-CholQR
		Pa Iter Err Iter	Err	Iter	Err	Iter	Err	Iter	Err
		2 200 4E-5 204	3E-5	204	3E-5	204	3E-5	204	3E-5
		4 167 2E-5 167	2E-5	167	2E-5	167	2E-5	167	2E-5
	POISSON2D 8 139 1E-5 139	1E-5	139	1E-5	139	1E-5	139	1E-5
	tol " 10 ´6 16 121 5E-6 121	5E-6	121	5E-6	121	5E-6	121	5E-6
		32 94 2E-6 94	2E-6	94	2E-6	94	2E-6	94	2E-6
		64 69 2E-6 69	2E-6	69	2E-6	69	2E-6	69	2E-6
		2 256 1E-7 256	1E-7	256	1E-7	256	1E-7	256	1E-7
		4 208 1E-7 208	1E-7	208	1E-7	208	1E-7	208	1E-7
	NH2D	8 169 8E-8 169	8E-8	169	8E-8	169	8E-8	169	8E-8
	tol " 10 ´8 16 138 6E-8 138	6E-8	138	6E-8	138	6E-8	138	6E-8
		32 107 2E-8 107	2E-8	107	2E-8	107	2E-8	107	2E-8
		64 77 1E-8 77	1E-8	77	1E-8	77	1E-8	77	1E-8
		2 1559 8E-4 -	-	-	-	1562	8E-4	1559	9E-4
		4 917 4E-4 -	-	-	-	917	4E-4	917	4E-4
	SKY2D	8 532 3E-4 -	-	-	-	531	2E-4	534	2E-4
	tol " 10 ´8 16 307 1E-4 -	-	-	-	307	1E-4	307	1E-4
		32 178 6E-5 -	-	-	-	178	6E-5	178	6E-5
		64 126 3E-6 -	-	-	-	124	2E-6	124	2E-6
		2 610 4E-5 611	4E-5	611	4E-5	611	4E-5	638	1E-5
		4 420 2E-5 -	-	-	-	424	1E-5	418	2E-5
	SKY3D	8 228 1E-5 -	-	-	-	230	1E-5	228	2E-5
	tol " 10 ´8 16 134 1E-5 -	-	-	-	134	1E-5	134	1E-5
		32 87 1E-6 -	-	-	-	83	1E-5	83	1E-5
		64 53 6E-6 -	-	-	-	51	1E-5	51	1E-5
		2 893 6e-5 893	6e-5	893	6e-5	893	6e-5	893	6e-5
		4 749 8e-5 749	8e-5	749	8e-5	749	8e-5	749	8e-5
	ANI3D	8 498 8e-5 506	9e-5	511	8e-5	498	7e-5	503	7e-5
	tol " 10 ´8 16 328 1e-4 -	-	-	-	326	1e-4	326	1e-4
		32 192 2e-4 -	-	-	-	192	1e-4	192	1e-4
		64 122 5e-5 -	-	-	-	122*	4e-5	122*	4e-5

Table 5 :

 5 The memory requirements in words, for MSDO-CG and LRE-CG with a varying number of partitions (Pa) or t and for different matrices.

		POISSON2D	NH2D		SKY2D	SKY3D	ANI3D	ELASTICITY3D
		Pa Iter Mem	Iter Mem	Iter	Mem	Iter Mem	Iter Mem	Iter	Mem
		2 200 4040004 256 5160004 1559 31220004 610 9792004 893 14320004 764 17239600
		4 167 6740016 208 8380016 917 36740016 420 13488016 749 24016016 622 28064998
	MSDO-CG	8 139 11220064 169 13620064 532 42660064 228 14672064 498 31952064 472 42603922 16 121 19540256 138 22260256 307 49300256 134 17296256 328 42128256 343 61959274
		32 94 30421024 107 34581024 178 57301024 87 22545024 192 49425024 234 84646090
		64 69 44824096 77 49944096 126 81304096 53 27668096 122 62996096		111111111
		2 193 4028996 245 5160100 1415 36328900 557 10168996 875 17078500 634 15899134
		4 153 6514544 188 8105504 757 39468784 373 14178064 673 28798864 480 25314666
	LRE-CG	8 123 10828256 149 13360864 398 41997856 211 16369344 449 41654464 334 37230106 16 95 17530400 112 21151264 220 47610400 119 18873216 253 48786304 235 56471386
		32 70 27437600 82 33145376 126 56597024 69 22555264 148 60333696 170 90832426
		64 52 44375584 60 53165600 75 71060000 43 29605504 92 81788544 117 140355114

  Algorithm 12 A-orthonormalization against previous vectors with MGS Flops Input: P 1 , P 2 ,.. , P k`1 , the k `1 sets of search directions Input: W 1 , W 2 ,.. , W k`1 , the k `1 sets of AP i Output: P k`1 , the search directions A-orthonomalized against P 1 , P 2 ,.. , P k 1: for o " 1 : t do %loop over the vectors of P k`1 P k`1 p:, oq " P k`1 p:, oq ´pW i p:, jq t P k`1 p:, oqqP i p:, jq 4n ´1 W k`1 p:, oq " W k`1 p:, oq ´pW i p:, jq t P k`1 p:, oqqW i p:, jq2nThen, the cost of the A-orthonormalization against previous vectors in Algorithm 12 is p6n ´1qt 2 k p4n `1qt of the order of 6nt 2 k flops.Algorithm 13 Flops reduction in A-orthonormalization against each others with MGS FlopsInput: P k`1 , the search directions to be A-orthonormalizedInput: W k`1 , AP k`1Output: P k`1 , the A-orthonomalized search directions Output: W k`1 , AP k`1 where P k`1 is the A-orthonomalized search directions 1: for i " 1 : t do %loop over the vectors of P k`1

	2:	for i " 1 : k do	%loop over the different Pi's
	3:	for j " 1 : t do %loop over the vectors of Pi
	4:
	5:
	6:	end for
	7:	end for
	8:	pap k`1 " W k`1 p:, oq t P k`1 p:, oq		2n	´1
	9:	P k`1 p:, oq "	P k`1 p:,oq ? pap k`1 and W k`1 p:, oq "	W k`1 p:,oq ? pap k`1	2n	`2
	10: end for
	2:

  InriaAs for the A-orthonormalization of P k`1 's vectors against each others using MGS with flops reductions, it can be performed as in Algorithm 13. The cost of this version of A-orthonormalization (Algorithm 13) is p6n ´1qpt ´1q t 2 `p4n `1qt, which is of the order of 3nt 2 .B.1.3 Parallelization of the A-orthonormalization using MGSIn Algorithm 12, at each inner iteration we are A-orthonormalizing the updated vectors P k`1 p:, oq against the vector P i p:, jq, where the vector P k`1 p:, oq is changed at each inner iteration. Thus it is not possible to have a block MGS by eliminating all the for loops. However, it is possible to eliminate one for loop in Algorithm 12 as shown in Algorithm 14, by A-orthonormalizing the whole block P k`1 against the vector P i p:, jq, where P k`1 p:, oq " P k`1 p:, oq ´pP k`1 p:, oq t W i p:, jqqP i p:, jq for all o " 1, 2, ..., t. Let rP i p:, jqs t be an n ˆt block containing t duplicates of the vector P i p:, jq. Then, P k`1 " P k`1 ´rP i p: , jqs t diagpP t k`1 W i p:, jqq. Algorithm 14 A-orthonormalization against previous vectors with MGS Flops Input: P 1 , P 2 ,.. , P k`1 , the k `1 sets of search directions Input: W 1 , W 2 ,.. , W k`1 , the k `1 sets of AP i Output: P k`1 , the search directions A-orthonomalized against P 1 , P 2 ,.. , P k 1: for i " 1 : k do

			%loop over the different Pi's
	2:
	4:	W k`1 p:, iq " W k`1 p:, iq ´pW k`1 p:, jq t P k`1 p:, iqqW k`1 p:, jq		2n
	5:	end for
	6:	pap k`1 " W k`1 p:, iq t P k`1 p:, iq		2n	´1
	7:	P k`1 p:, iq "	P k`1 p:,iq ? pap k`1 and W k`1 p:, iq "	W k`1 p:,iq ? pap k`1	2n	`2
	8: end for

  InriaAlgorithm 16 A-orthonormalization against previous vectors with CGS Input: A, the n ˆn symmetric positive definite matrix Input: P 1 , P 2 ,.. , P k`1 , the k `1 sets of search directions Output: r P k`1 , the search directions A-orthonomalized against P 1 , P 2 ,.. , P k P k`1 p:, oq ´pP i p:, jq t AP k`1 p:, oqqP i p:, jq

	1: Let r P k`1 " P k`1
	2: for o " 1 : t do		%loop over the vectors of P k`1
	3:	for i " 1 : k do	%loop over the different Pi's
	4:	for j " 1 : t do	%loop over the vectors of Pi
	5:	r P k`1 p:, oq " r
	6:	end for
	7:	end for
	8:	r P k`1 p:, oq "	r P k`1 p:,oq || r P k`1 p:,oq|| A	"	?	r P k`1 p:,oq P k`1 p:,oq t A r r P k`1 p:,oq	%A-normalize
	9: end for
	Algorithm 17 A-orthonormalization against previous vectors with BCGS	Flops

Input: A, the n ˆn symmetric positive definite matrix Input: Q k " rP 1 , P 2 , ..., P k s, the tk search directions Input: P k`1 , the t search directions to be A-orthonormalized Output: r P k`1 , the search directions A-orthonomalized against P 1 , P 2 ,.. , P k

  BCGS Aort « γ c 6ntk `αc 2logptq `βc t 2 klogptqAlgorithm 18 A-orthonormalization against previous vectors with BCGS FlopsInput: Q k " rP 1 , P 2 , ..., P k s, the tk search directions Input: P k`1 , the t search directions to be A-orthonormalizedInput: W k`1 " AP k`1 ; W k " AQ k Output: r P k`1, the search directions A-orthonomalized against P 1 , P 2 ,.. , P k

	p6n´1qt 2 k`3nt t	" p6n ´1qtk `3n flops are performed in parallel. Hence, by ignoring lower order terms
	we obtain
		T ime

  t do Given a set of vectors P k`1 that are A-normalized, i.e the diagonal of P t k`1 AP k`1 is equal to ones, we A-orthonormalize it (P t k`1 AP k`1 " I) using a classical Gram Schmidt procedure as in algorithm 19.

			%loop over the vectors of P k`1 and A-normalize
	4:	r P k`1 p:, iq "	r P k`1 p:,iq || r P k`1 p:,iq|| A	"	?	r P k`1 p:,iq P k`1 p:,iq t r Ă W k`1 p:,iq	3n
	5: end for
	B.2.2 A-orthonormalization of a set of vectors using CGS
			%loop over the vectors of P k`1
	3:	for j " 1 : pi ´1q do	%A-orthogonalize against the vectors P k`1 p:, 1 : i ´1q
	4:
	5:	end for
	6:	r P k`1 p:, iq "	r P k`1 p:,iq || r P k`1 p:,iq|| A	"	?	r P k`1 p:,iq P k`1 p:,iq t A r r P k`1 p:,iq	%A-normalize

Algorithm 19 A-orthonormalization against each others using CGS Input: A, the n ˆn symmetric positive definite matrix Input: P k`1 , the search directions to be A-orthonormalized Output: P k`1 , the A-orthonomalized search directions 1: Let r P k`1 " P k`1 2: for i " 1 : t do r P k`1 p:, iq " r P k`1 p:, iq ´p r P k`1 p:, jq t AP k`1 p:, iqq r P k`1 p:, jq

  The total flops of Algorithm 20 is of the order of nnzt `nt 2 .

					´nqt
	2: Rp1, 1q "	a P k`1 p:, 1q t W p:, 1q	2n
	3: r P k`1 p:, 1q "	P k`1 p:,1q Rp1,1q	n
	4: for i " 2 : t do
	5:	Rpi ´1, i : tq " r P k`1 p:, i ´1q t W k`1 p:, i : tq	p2n ´1qpt ´i `1q
	6:	r P k`1 p:, iq " P k`1 p:, iq ´r P k`1 p:, 1 : i ´1qRp1 : i ´1, iq	p2pi ´1q ´1qn	`n
				b
	7:	Rpi, iq "	r P k`1 p:, iq t A r P k`1 p:, iq	p2nnz ´nq `2n
	8:	r P k`1 p:, iq "	r P k`1 p:,iq Rpi,iq	n
	9: end for
		Total " 2nnzt ´nt `3n	`řt i"2 rp2n ´1qpt ´i `1q `2pi ´1qn `p2nnz `2nqs
		" 2nnzt ´nt `3n	`řt i"2 rp2n ´1qpt `1q ´p2n ´1qi `2ni ´2n `2nnz `2ns
		" 2nnzt ´nt `3n	`řt i"2 rp2n ´1qpt `1q `i `2nnzs
		" 2nnzt ´nt `3n `rp2n ´1qpt `1q `2nnzspt ´1q `tpt`1q 2	´1
		" 4nnzt ´2nnz ´nt `3n `p2n ´1qpt 2 ´1q	`t2 `t 2	´1
		" 4nnzt ´2nnz ´nt `n `2nt 2 ´t2 `t2 `t 2 " 4nnzt ´2nnz ´nt `n `2nt 2 `´t 2 `t 2

  T ime QRCGS Aort « γ c 3nt `αc 2tlogptq `βc t 2 logptqRR n°8597Algorithm 21 QR factorization with A inner product using CGS FlopsInput: P k`1 , the search directions to be A-orthonormalized; W k`1 , " AP k`1 Output: r P k`1 , the A-orthonomalized search directions; Ă W k`1 Output: R, the upper triangular matrix such that P k`1 " r P k`1 R 1: Rp1, 1q " Rpi ´1, i : tq " r P k`1 p:, i ´1q t W k`1 p:, i : tq p2n ´1qpt ´i `1q 5: r P k`1 p:, iq " P k`1 p:, iq ´r P k`1 p:, 1 : i ´1qRp1 : i ´1, iq Ă W k`1 p:, iq " W k`1 p:, iq ´Ă W k`1 p:, 1 : i ´1qRp1 : i ´1, iq 2pi ´1qn

		a P k`1 p:, 1q t W p:, 1q	2n
	2: r P k`1 p:, 1q "	P k`1 p:,1q Rp1,1q and Ă W k`1 p:, 1q "	W k`1 p:,1q Rp1,1q	2n
	3: for i " 2 : t do
	4:
						p2pi ´1q ´1qn	`n
	6:
				b
	7:	Rpi, iq "	r P k`1 p:, iq t	Ă W k`1 p:, iq	2n
	8:	r P k`1 p:, iq "	r P k`1 p:,iq Rpi,iq and Ă W k`1 p:, iq "	Ă W k`1 p:,iq Rpi,iq	2n
	9: end for
	B.2.3 CGS with Reorthogonalization (CGS2)
	tpt`1q 2 logptq words where	p1´tq 3nt 2 `nt`t 2 t	" 3nt `n `p1´tq 2	flops are performed in parallel. Hence, by
	ignoring lower order terms we obtain

  ´1AL ´t p P k`1 9: for i " 1 : t do %loop over the vectors of P k`1 and L ´1 AL ´t-normalize A, the n ˆn symmetric positive definite matrix Input: P k`1 , the search directions to be A-orthonormalized Output: rP k`1 , the A-orthonomalized search directions 1: Compute W k`1 " AP k`1 p2nnz ´nqt 2: Compute C " W t k`1 P k`1 p2n ´1qt 2 3: Compute the Cholesky factorization of C " R t R to obtain R P k`1 , the search directions to be A-orthonormalized, W k`1 " AP k`1 Output: r P k`1 , the A-orthonomalized search directions; Ă W k`1 " A r P k`1 1: Compute C " W t k`1 P k`1 p2n ´1qt 2 2:Compute the Cholesky factorization of C " R t R to obtain R

	5:	r P k`1 p:, iq "	r P k`1 p:,iq P k`1 p:,iq|| L ´1 AL || r ´t "	?	r P k`1 p:,iq P k`1 p:,iq t r Ă W k`1 p:,iq	and Ă W k`1 p:, iq "	?	Ă W k`1 p:,iq P k`1 p:,iq t r Ă W k`1 p:,iq
	6: end for
	7: p P k`1 " r P k`1 ´Qk pQ t k	Ă W k`1 q
	10:	p P k`1 p:, iq "	p P k`1 p:,iq P k`1 p:,iq|| L ´1 AL || p ´t "	?	p P k`1 p:,iq P k`1 p:,iq t p x W k`1 p:,iq	and x W k`1 p:, iq "	?	x W k`1 p:,iq P k`1 p:,iq t p x W k`1 p:,iq
	11: end for
	Algorithm 24 A-CholQR					Flops
		Input: t 2
	4: Solve R t	r P t k`1 " P t k`1					nt 2
	5: for i " 1 : t do 6: r P k`1 p:, iq "	r P k`1 p:,iq || r P k`1 p:,iq|| A	"	?	r P k`1 p:,iq P k`1 p:,iq t A r r P k`1 p:,iq	3n
	7: end for
	Algorithm 25 A-CholQR					Flops
		Input: t 2
	3: Solve R t	r P t k`1 " P t k`1					t 2 n
	4: Solve R t	Ă W t k`1 " W t k`1			t 2 n
	5: for i " 1 : t do 6: r P k`1 p:, iq "	r P k`1 p:,iq || r P k`1 p:,iq|| A	"	?	r P k`1 p:,iq P k`1 p:,iq t r Ă W k`1 p:,iq	and	4n
			Ă W k`1 p:, iq "	?	Ă W k`1 p:,iq P k`1 p:,iq t r Ă W k`1 p:,iq
	7: end for
	RR n°8597

8: x W k`1 " L

  26 p A-CholQRInput: A, the n ˆn symmetric positive definite matrix; L, n ˆn preconditioner Input: P k`1 , the search directions to be A-orthonormalized;W k`1 , " L ´1AL ´tP k`1 Output: r P k`1 , the A-orthonomalized search directions; Ă W k`1 " L ´1AL ´t r P k`1 1: Compute C " W t k`1 P k`1 2:Compute the Cholesky QR factorization of C to obtain R

	3: Solve R t	r P t k`1 " P t k`1
	4: Compute Ă W k`1 " L ´1AL ´t r P k`1
	5: for i " 1 : t do
	6:	r P k`1 p:, iq "	r P k`1 p:,iq || r P k`1 p:,iq|| A	"	?	r P k`1 p:,iq

r P k`1 p:,iq t Ă W k`1 p:,iq and Ă W k`1 p:, iq " Ă W k`1 p:,iq ? r P k`1 p:,iq t Ă W k`1 p:,iq 7: end for Recently, Lowery and Langou presented a new version of A-CholQR in

  2 qlogptq words. Hence, by ignoring lower order terms we obtain A, the n ˆn symmetric positive definite matrix Input: P k`1 , the search directions to be A-orthonormalized Output: r P k`1 , the A-orthonomalized search directions 1: Compute the QR factorization of P k`1 " P

	T ime P reCholQR « γ c p6nt	``2 3	t 3 logptqq `αc 2logptq `βc pnt	`3 2	t 2 logptqq
	Algorithm 27 Pre-CholQR
	Input: 1 k`1 R
	2: Call Algorithm 24 with A and P k`1 as input and with r 1 P k`1 as output

Inria

RESEARCH CENTRE PARIS -ROCQUENCOURT

Domaine de Voluceau, -Rocquencourt B.P. 105 -78153 Le Chesnay Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr

RR n°8597

and SKY2D is 1 ˆ10 6 , 2 ˆ10 6 , and 3 ˆ10 7 respetively. Although it was shown in [START_REF] Lowery | Stability Analysis of QR factorization in an Oblique Inner Product[END_REF] that Pre-CholQR A-orthonormalization is more stable than A-CholQR, however MSDO-CG with CGS+A-CholQR Aorthonormalization and MSDO-CG with CGS+Pre-CholQR A-orthonormalization are both numerically unstable.

Thus, we replace CGS with CGS2 where the A-orthonormalization is performed twice for numerical stability. Then, the MSDO-CG with CGS2+A-CholQR A-orthonormalization and MSDO-CG with CGS2+Pre-CholQR A-orthonormalization converge as fast as MSDO-CG with MGS A-orthonormalization for all t and all the tested matrices. Hence, we concude that MGS, CGS2+A-CholQR, and CGS2+Pre-CholQR A-orthonormalizations are stable enough to be used in the MSDO-CG method (Algorithm 1).

Convergence

We compare the convergence of split preconditioned MSDO-CG with CGS2+CholQR p A-orthonormalization and split preconditioned LRE-CG with CGS+TSQR orthonormalization to CG and split preconditioned CG (PCG). We use Block Jacobi with Cholesky factorization of the block diagonals as a preconditioner. In table 6, we use a different Block Jacobi preconditioner for the different partitions. First, the graph of A is partitioned into t parts that define the enlarged Krylov subspace using Metis's kway edge separator where t " 2, 4, 8, 16, 32, 64 . Then the Block Jacobi preconditioner M is defined as the t diagonal blocks of the permuted matrix A. Each of the t blocks is factorized using Cholesky decomposition to obtain a block L. The preconditioned LRE-CG converges faster than the preconditioned MSDO-CG and PCG for the different configurations. As the number of partitions or the maximum basis vectors added at each iteration is doubled, the Block Jacobi preconditioned CG needs more iterations to converge. However, for the matrices POISSON2D, NH2D, and SKY2D, the number of iteration of the preconditioned LRE-CG and MSDO-CG decreases. As for the matrices SKY3D and ANI3D, the number of iterations of LRE-CG increases then decreases back to the same number of iterations for t " 2, unlike preconditioned MSDO-CG. In table 7, we use a fixed Block Jacobi preconditioner for all the partitions to compare the convergence of the methods with respect to the doubling of the number of partitions t. First, the graph of A is partitioned into 64 parts using Metis's kway edge separator. Then the Block Jacobi preconditioner M is defined as the 64 block diagonals of the permuted matrix A. Each of the 64 blocks is factorized using Inria Appendix A MCD-CG Algorithm Algorithm 9 The MSD-CG Algorithm Input: A, the n ˆn symmetric positive definite matrix Input: b, the n ˆ1 right-hand side; x 0 , the initial guess or iterate Input: , the stopping tolerance; k max , the maximum allowed iterations Output: x k , the approximate solution of the system Ax " b 1: r 0 " b ´Ax , ρ " ||r 0 || 2 2 , k " 1 2: for i " 1, .., t do Let P p:, iq " T i pr 0 q 3: end for 4: while ( ? ρ ą ||b|| 2 and k ă k max ) do 5:

α " pP t AP q ´1pP t rq 6:

x " x `P α 7:

r " r ´AP α 8:

β " pP t AP q ´1pP t Arq 9:

for i " 1, .., t do Let P p:, iq " T i prq `βpiqP p:, iq k " k `1 13: end while Appendix B The A-orthogonalization of search directions

In this section we describe the A-orthogonalization of the t newly computed vectors P k`1 against all the previous vectors P i 's for i ă k `1, and then against each others. Here A is a symmetric positive definite matrix. The A-orthogonalization is simply an orthogonalization with the A inner product (ă . , . ą A " ă . , A . ą) rather than the L2 inner product (ă . , . ą). A-orthonormalizing a tall and skinny n ˆt matrix P k`1 , or alternatively computing the oblique QR factorization of P k`1 , has been discussed in [START_REF] Rozloznik | Numerical stability of orthogonalization with a non-standard inner product[END_REF] and [START_REF] Lowery | Stability Analysis of QR factorization in an Oblique Inner Product[END_REF] in terms of stability and ease of parallelization. The goal is to get a r P k`1 , such that r P t k`1 A r P k`1 " I. There are two main classes for computing this oblique QR factorization of P k`1 " r P k`1 R. The first class is to factorize the matrix A " B t B using Cholesky decomposition or eigenvalue decomposition, which is expensive. Then P t k`1 AP k`1 " pBP k`1 q t pBP k`1 q, where the oblique QR factorization of P k`1 is transformed into a Eucleadean QR factorization of the matrix BP k`1 " Q B R B with r P k`1 " B ´1Q B and R " R B . The second class consists of avoiding any factorization of A, like CGS, CGS2, MGS, and the Cholesky factorization of the t ˆt matrix P t k`1 AP k`1 . For A-orthonormalizing P k`1 against all the previous vectors P i with i ă k `1, it is possible to use CGS, CGS2, MGS and A-choleskyBGS which was discussed in Hoemmen's thesis ( [START_REF] Hoemmen | Communication-Avoiding Krylov Subspace Methods[END_REF], page 115).

Thus, we start by discussing the A-orthonormalization using modified Gram Schmidt in section B.1. However, this version is not easily parallelized and requires a lot of communication (ptk `1qlogptq 2pt ´1qlogptq messages) as compared to the classical Gram Schmidt version.Then, in section B.2 we adapt the A-orthonormalization of the vectors of P k`1 against P i 's for i ă k `1 using the classical Gram Schmidt (CGS) to obtain a Block Gram Schmidt (BGS) version (Algorithm 14) with A inner product that requires only 2logptq messages. As for the A-orthonormalization of the P k`1 vectors against each others, we introduce a parallelizable version with reduced communication (p2t ´1qlogptq messages ). Note that CGS2, section B.2.3, consists of calling the algorithm CGS twice. Thus its cost is twice the cost of CGS. In section B.3, we briefly discuss the A-orthonormalization of P k`1 using the Cholesky factorization (CholQR) of the t ˆt matrix P t k`1 AP k`1 which is referred to as A-CholQR and requires Inria