A flexible Generalized Conjugate Residual method with inner orthogonalization and deflated restarting

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

situation is frequent when solving very large systems of linear equations resulting from the discretization of partial differential equations in three dimensions. Thus flexible Krylov subspace methods have gained a considerable interest in the recent years and are the subject of both theoretical and numerical studies [27]. We refer the reader to [29, section 10] for additional comments on flexible methods.
When nonvariable preconditioning is considered, the full GMRES method [23] is often chosen for the solution of nonsymmetric or non-Hermitian linear systems because of its robustness and its minimum residual norm property [26]. Nevertheless to control both the memory requirements and the computational cost of the orthogonalization scheme, restarted GMRES is preferred; it corresponds to a scheme where the maximal dimension of the approximation subspace is fixed. It means in practice that the orthonormal basis built is thrown away at the end of the cycle. Since some information is discarded at the restart, the convergence may stagnate and is expected to be slower compared to full GMRES. Nevertheless to retain the convergence rate a number of techniques have been proposed; they fall in the class of augmented and deflated methods; see, e.g., [4,10,11,16,25]. Deflated methods compute spectral information at a restart and use this information to improve the convergence of the subspace method. One of the most recent procedures based on a deflation approach is GMRES with deflated restarting (GMRES-DR) [18]. This method reduces to restarted GM-RES when no deflation is applied, but may provide a much faster convergence than restarted GMRES for well-chosen deflation spaces as described in [18].
Quite recently a new minimum residual norm subspace method based on GMRES allowing deflated restarting and variable preconditioning has been proposed in [15]. It mainly attempted to combine the numerical features of GMRES-DR and the flexibility property of FGMRES. Numerical experiments in [15] have shown the efficiency of FGMRES with deflated restarting (FGMRES-DR) on both academic and industrial examples. In this paper we study a new minimum residual norm subspace method based on the generalized conjugate method with inner orthogonalization (GCRO) [9] allowing deflated restarting and variable preconditioning. It is named flexible generalized conjugate residual method with inner orthogonalization and deflated restarting (FGCRO-DR) and can be viewed as an extension of GCRO-DR [22] to the case of variable preconditioning. A major advantage of FGCRO-DR over FGMRES-DR is its ability to solve sequences of linear systems (where both the left-and right-hand sides can change) through recycling [22]. In [22] Parks et al. mentioned that GCRO-DR and GMRES-DR were algebraically equivalent, i.e., both methods produce the same iterates in exact arithmetic when solving the same given linear system starting from the same initial guess. When variable preconditioning is considered, it seems therefore natural to ask whether FGCRO-DR and FGMRES-DR could also be algebraically equivalent. We address this question in this paper, and the main theoretical developments that are proposed will help us to answer this question. The main contributions of the paper are then twofold. First we prove that FGCRO-DR and FGMRES-DR can be considered algebraically equivalent if a collinearity condition between two certain vectors is satisfied at each cycle. When considering nonvariable preconditioning, these theoretical developments will also allow us to show the algebraic equivalence between GCRO-DR and GMRES-DR that was stated without proof in [22]. Second we carefully analyze the computational cost of FGCRO-DR and show that the proposed method is nearly as expensive as FGMRES-DR in terms of operations per cycle. Furthermore it is explained how to include subspace recycling into FGCRO-DR, and numerical experiments are reported showing the efficiency of FGCRO-DR. This paper is organized as follows. In section 2 we introduce the general background of this study. We briefly recall the main properties of FGMRES-DR and then introduce the FGCRO-DR method both from a mathematical and an algorithmic point of view. Section 3 is mainly devoted to the analysis of both flexible methods. Therein we show that both methods can be algebraically equivalent in the flexible case if a certain collinearity condition is satisfied at each cycle. In section 4 we compare FGCRO-DR and FGMRES-DR in terms of computational operations per cycle and storage and discuss the solution of sequences of linear systems through subspace recycling. Finally we draw some conclusions and perspectives in section 5.

Flexible
Krylov methods with restarting.

General setting.
Notation. Throughout this paper we denote by . the Euclidean norm, by I k ∈ C k×k the identity matrix of dimension k, and by 0 i×j ∈ C i×j the zero rectangular matrix with i rows and j columns. Given N ∈ C n×m , Π N ⊥ = I n − N N † will represent the orthogonal projector onto range(N ) ⊥ , where the superscript † refers to the Moore-Penrose inverse. Finally, given Z m = [z 1 , . . . , z m ] ∈ C n×m , we will usually decompose Z m into two submatrices defined as Setting. We focus on minimum residual norm based subspace methods that allow flexible preconditioning for the iterative solution of given an initial vector x 0 ∈ C n . In this paper A is supposed to be nonsingular. Flexible methods refer to a class of methods where the preconditioner is allowed to vary at each iteration. We refer the reader to, e.g., [29] for a general introduction on Krylov subspace methods and to [29, section 10] and [26, section 9.4] for a review on flexible methods. The minimum residual norm GMRES method [23] has been extended by Saad [24] to allow variable preconditioning. The resulting algorithm known as FGMRES(m) relies on the Arnoldi relation where Z m ∈ C n×m , V m+1 ∈ C n×(m+1) has orthonormal columns andH m ∈ C (m+1)×m is upper Hessenberg. We denote by M j the preconditioning operator at iteration j and remark that M j may be a nonlinear preconditioning function. We will then denote by M j (v) the action of M j on a vector v. In (2.2), the columns of V m+1 form an orthonormal basis of the subspace spanned by the vectors At the end of the cycle an approximate solution x m ∈ C n is then found by minimizing the residual norm r 0 − AZ m y over the space x 0 + range(Z m ). Thus we obtain that where y * is the solution of the following least-squares problem of size (m + 1) × m: where e 1 is the first canonical vector of C m+1 . Flexible subspace methods with restarting are based on a procedure where the construction of the subspace is stopped after a certain number of steps (denoted by m in this paper with m < n). The method is then restarted mainly to control both the memory requirements and the cost of the orthogonalization scheme. In FGMRES(m) the restarting consists in taking as an initial guess the last iterate of the cycle (x m ).
The main focus of this paper is to present minimum residual norm subspace methods with deflated restarting that allow flexible preconditioning. Deflated restarting aims at determining an approximation subspace of dimension m as a direct sum of two subspaces of smaller dimension, where one of these subspaces will contain relevant spectral information that will be kept for the next cycle. We refer the reader to, e.g., [25] and [29, section 9] for a review of augmented and deflated methods. Flexible methods with deflated restarting will notably satisfy the following flexible Arnoldi relation: whereH m ∈ C (m+1)×m is not necessarily of upper Hessenberg form. In this paper we call this relation a flexible Arnoldi-like relation due to its similarity to relation (2.2). Stagnation and breakdown. We refer the reader to [27, section 6] for general comments and a detailed discussion on the possibility of both breakdown and stagnation in flexible inner-outer Krylov subspace methods. Although important, these issues are not addressed in this paper, and we assume that no breakdown occurs in the inner-outer subspace methods that will be proposed.

Flexible GMRES with deflated restarting.
A number of techniques have been proposed to compute spectral information at a restart and use this information to improve the convergence rate of the Krylov subspace methods; see, e.g., [16,17,18,25]. These techniques have been exclusively developed in the case of a fixed preconditioner. GMRES-DR is one of these methods. It focuses on removing (or deflating) the eigenvalues of smallest magnitude. A full subspace of dimension k, k < m (and not only the approximate solution with minimum residual norm) is now retained at the restart, and the success of this approach has been demonstrated in many academic examples [16]. Approximations of eigenvalues of smallest magnitude are obtained by computing harmonic Ritz pairs of A with respect to a certain subspace [18]. We present here a definition of a harmonic Ritz pair equivalent to the one introduced in [21,30]; it will be of key importance when defining appropriate deflation strategies.
Definition 2.1 (harmonic Ritz pair). Consider a subspace U of C n . Given B ∈ C n×n , θ ∈ C, and y ∈ U, (θ, y) is a harmonic Ritz pair of B with respect to U if and only if By − θ y ⊥ B U or, equivalently, for the canonical scalar product, We call y a harmonic Ritz vector associated with the harmonic Ritz value θ.
As in the case of fixed preconditioning, deflated restarting may also improve the convergence rate of flexible subspace methods. In [15] a deflated restarting procedure has been proposed for the FGMRES algorithm. The ith cycle of the resulting algorithm, called FGMRES-DR, is now briefly described, and we denote by , V m+1 ,H m , and Z m the residual and matrices obtained at the end of the (i − 1)th cycle.
Based on the Arnoldi-like relation (2.3), the deflation procedure proposed in [15, Proposition 1] relies on the use of k harmonic Ritz vectors Y k = V m P k of AZ m V H m with respect to range(V m ), where Y k ∈ C n×k and P k ∈ C m×k . In Lemma 2.2 shown in [15,Lemma 3.1], we recall a useful relation satisfied by the harmonic Ritz vectors.
Lemma 2.2. In FGMRES-DR, the harmonic Ritz vectors are given by Y k = V m P k with corresponding harmonic Ritz values λ k . P k ∈ C m×k satisfies the following relation: Next, the QR factorization of the (m + 1) × (k + 1) matrix appearing on the right-hand side of relation (2.4) is performed as where Q ∈ C (m+1)×(k+1) has orthonormal columns and R ∈ C (k+1)×(k+1) is upper triangular, respectively. We write the matrix Q obtained in relation (2.6) as where Q m×k ∈ C m×k andρ ∈ C m+1 is defined as

Proposition 1. In FGMRES-DR, the flexible Arnoldi relation
holds at the ith cycle with matrices Z k , V k ∈ C n×k andH k ∈ C (k+1)×k defined as where V m+1 , Z m , andH m refer to matrices obtained at the end of the (i − 1)th cycle.
FGMRES-DR then carries out m − k Arnoldi steps with flexible preconditioning and starting vector v k+1 while maintaining orthogonality to V k , leading to We note thatH m−k ∈ C (m−k+1)×(m−k) is upper Hessenberg. At the end of the ith cycle this gives the flexible Arnoldi-like relation We note thatH m is no longer upper Hessenberg due to the leading dense (k +1)×k submatrixH k . At the end of the ith cycle, an approximate solution x (i) 0 ∈ C n is then found by minimizing the residual . We refer the reader to [15] for the complete derivation of the method and numerical experiments showing the efficiency of FGMRES-DR in both academic and industrial examples.

Flexible GCRO with deflated restarting. GCRO-DR [22]
-a combination of GMRES-DR and GCRO-is a Krylov subspace method that allows deflated restarting and subspace recycling simultaneously. This latter feature is particularly interesting when solving sequences of linear systems with possibly different left-or right-hand sides. As pointed out in [22], GCRO-DR is attractive because any subspace may be recycled. In this paper we restrict the presentation to the case of a single linear system as proposed in (2.1).
GCRO and GCRO-DR belong to the family of inner-outer methods [2,Chap. 12] where the outer iteration is based on GCR, a minimum residual norm method proposed by Eisenstat, Elman, and Schultz [13]. To this end GCR maintains a correction subspace spanned by range(Z m ) and an approximation subspace spanned by The optimal solution of the minimization problem min b − Ax over the subspace . In [9] de Sturler proposed an improvement to GMRESR [34], an inner-outer method based on GCR in the outer part and GMRES in the inner part, respectively. He suggested that the inner iteration takes place in a subspace orthogonal to the outer Krylov subspace. In this inner iteration the projected residual equation is solved only approximately. If a minimum residual norm subspace method is used in the inner iteration to solve this projected residual linear system, the residuals over both the inner and outer subspaces are minimized. This leads to the GCRO Krylov subspace method [9]. Numerical experiments [9] indicate that the resulting method may perform better than other inner-outer methods (without orthogonalization) in some cases. The GCRO method with deflated restarting (named GCRO-DR) based on harmonic Ritz value information was proposed in [22]. An approximate invariant subspace is used for deflation following closely the GMRES-DR method. We refer the reader to [22] for a description of this method, algorithms, and implementation details. We present now a new variant of GCRO-DR that allows flexible preconditioning by explaining the different steps occurring during the ith cycle. Again we denote by r ,H m , and Z m the residual and matrices obtained at the end of the (i − 1)th cycle.
We suppose that a flexible Arnoldi-like relation of type (2.3) holds. As in section 2.2 an important point is to specify which harmonic Ritz information is selected. Given a certain matrix W m ∈ C n×m to be specified later on, such as range(W m ) = range(V m ), the deflation procedure relies on the use of k harmonic Ritz vectors all belong to a subspace of dimension m + 1 (spanned by the columns of V m+1 ) and are orthogonal to the same subspace of dimension m (spanned by the columns of AZ m subspace of range(V m+1 )), so they must be collinear. Consequently there exist k coefficients noted β j ∈ C with 1 ≤ j ≤ k such that , the collinearity expression (2.17) can be written in matrix form as Due to the flexible Arnoldi-like relation (2.3), relation (2.16) can be also expressed as If required, β 1×k can be deduced from (2.18) by Next, the QR factorization of the (m + 1) × k matrixH m P k appearing in relation (2.18) is performed asH m P k = QR with Q ∈ C (m+1)×k and R ∈ C k×k . Proposition 2. In FGCRO-DR, the relation where V m+1 and Z m refer to matrices obtained at the end of the (i − 1)th cycle. In addition V H k r (i−1) 0 = 0 holds during the ith cycle. Proof. By using information related to the QR factorization ofH m P k and the flexible Arnoldi relation (2.3) exclusively, we obtain This finally shows that V H k r (i−1) 0 = 0 since R is supposed to be nonsingular. To complement the subspaces, the inner iteration is based on the approximate solution of where the last equality is due to Proposition 2. For that purpose FGCRO-DR then carries out m − k steps of the Arnoldi method with flexible preconditioning, leading to . At the end of the cycle this gives the flexible Arnoldi-like relation where Z m ∈ C n×m , V m+1 ∈ C n×(m+1) , andH m ∈ C (m+1)×m . At the end of the ith cycle, an approximate solution x (i) 0 ∈ C n is then found by minimizing the residual

Algorithms.
Details of the FGCRO-DR method are given in Algorithm 1, where MATLAB-like notations are adopted (for instance, in step 7b, Q(1 : m, 1 : k) denotes the submatrix made of the first m rows and first k columns of matrix Q noted Q m×k in (2.7)). For the sake of completeness the FGMRES-DR algorithm has also been described with notation chosen as closely as possible to FGCRO-DR to make code comparison easier. Concerning Algorithm 1 we make the following comments: • As will be discussed later, the computation of W † m in step 5a is not required thanks to the definition of the harmonic Ritz pair (see Definition 2.1). • As pointed out by Morgan [18] and Parks et al. [22] we might have to adjust k during the algorithm to include both the real and imaginary parts of complex eigenvectors. 1: Choose m, k, tol, and x 0 2:

Analysis of FGMRES-DR and FGCRO-DR.
We compare now the flexible variants of GMRES-DR and GCRO-DR introduced in sections 2.2 and 2.3, respectively. In the following we use to denote quantities related to the FGMRES-DR algorithm, e.g., Y k denotes the set of harmonic Ritz vectors computed in the FGMRES-DR algorithm. When analyzing both algorithms we will suppose that identical preconditioning operators are used in steps 10a and 10b, respectively, i.e.,

Equivalent preconditioning matrix.
obtained during a cycle of a flexible method with (standard or deflated) restarting (with 1 ≤ p ≤ m < n) are both of full rank, i.e., rank V p = rank Z p = p. We will then denote by M V p ∈ C n×n a nonsingular equivalent preconditioning matrix defined as Such a matrix represents the action of the nonlinear operators M j on the set of vectors v j (with j = 1, . . . , p). It can be chosen, e.g., as the equivalent preconditioning matrices used in the initialization phase of both algorithms (step 3 in Algorithm 1). With this notation we remark that the following relations hold:

Relations between
We first analyze the relation between Z m and V m .
Proof. The initialization phase leads to the relation Z m = M We suppose that at the end of the (i − 1)th cycle the following relation holds: At step 9b of the ith cycle, Z k is defined as The proof is then complete since The next lemma details a relation between Z m and W m that is satisfied in the FGCRO-DR method.
Proof. The initialization phase leads to the relation Z m = M W m W m . We suppose that at the end of the (i − 1)th cycle the following relation holds: At step 9a of the ith cycle, Z k is defined as Proof. We use the relation AZ k = V k satisfied in the FGCRO-DR method shown in Proposition 2. The proof is then complete since We conclude this section by presenting a technical lemma related to the FGMRES-DR method.
Lemma 3.5. During the ith cycle of the FGMRES-DR method, v k+1 satisfies the relation denotes the residual obtained at the end of the (i − 1)th cycle.
Proof. Using Proposition 1 and relation (2 Since V m Q m×k has orthonormal columns, this last expression now becomes Because Q m×k is the orthogonal factor of the QR decomposition of P k , we obtain range( V m P k ) = range( V m Q m×k ).
Since from Lemma 2.3 Y k = V m P k , the proof is then complete.

3.3.
Analysis of the FGMRES-DR and FGCRO-DR methods. Lemma 3.3 has already described an important property satisfied by W m in the FGCRO-DR method proposed in Algorithm 1. We will analyze further the relation between the FGMRES-DR and FGCRO-DR methods. The next theorem states that the two flexible methods generate the same iterates in exact arithmetic under some conditions involving notably two vectors.
. Under this assumption the harmonic Ritz vectors Y k and Y k can be chosen equal during the (i +1)th cycle. If in addition there exists a real-valued positive coefficient η i+1 such that in the FGCRO-DR algorithm, then both algorithms generate the same iterates in exact arithmetic and
Proof. The whole proof is performed in three parts assuming that we analyze the (i + 1)th cycle of each algorithm. Suppose that at the beginning of the (i + 1)th cycle (step 4) there exist a unitary matrix Q ∈ C (k+1)×(k+1) and a nonsingular matrix X ∈ C (k+1)×(k+1) such that the following relations hold: (3.14) [ We will then prove the existence of a unitary matrix Q ∈ C (k+1)×(k+1) and of a nonsingular matrix X ∈ C (k+1)×(k+1) such that at the end of the (i + 1)th cycle Regarding FGCRO-DR we assume that at the beginning of the (i + 1)th cycle (step 4) range(W m ) = range(V m ). (3.21) We will also prove that relation (3.21) holds at the end of the (i+1)th cycle. Note that relations (3.9), (3.10), and (3.21) are obviously satisfied before the first cycle, because steps 1 to 3 are identical in both algorithms, yielding V m+1 = V m+1 , Z m = Z m , and W m = V m . Finally a consequence of (3.13), (3.15), (3.14), and (3.16) is that the residuals of the linear system Ax = b in both algorithms are equal at the beginning of the (i + 1)th cycle, i.e., r (i) 0 = r (i) 0 . We will denote by r 0 this residual for ease of notation.
Part I: Steps 5a and 5b. In this part, we prove that we can choose FGCRO-DR. Let y j = W m p j be the jth column of Y k . Since y j is a harmonic Ritz vector of AZ m W † m with respect to range(W m ), the following relation holds (see Definition 2.1): Due to (3.14) and (3.16) there exists a nonsingular matrix X ∈ C m×m that relates Z m and Z m :

From Lemma 3.3 and relation (3.23) we deduce
where we have used explicitly the assumption on the equivalent preconditioning matrix obtained at the end of the ith cycle, i.e., M (i) Next, the application of Lemma 3.2 leads to Since X is nonsingular the last equality proves that V m Xp j is a harmonic Ritz vector of A Z m V H m with respect to range( V m ) associated to the Ritz value θ j . From relations (3.22) and (3.24) we deduce that the harmonic Ritz vectors can be chosen to be equal and correspond to the same harmonic Ritz values. In this case they notably satisfy the following equality: ∀j ∈ {1, . . . , k}, V m Xp j = W m p j , i.e., p j = Xp j . (3.25) We will then denote by Y = Y k = Y k the k harmonic Ritz vectors computed in either FGCRO-DR or FGMRES-DR. We assume that the harmonic Ritz values θ j (1 ≤ j ≤ k) are nonzero.
Part IIa: Steps 6a-10a, 6b-10b. We show that at the end of steps 10a and 10b the following relations hold: range(V k+1 ) = range( V k+1 ) = range([Y, r (i) 0 / r (i) 0 ]). This result will help us to prove the existence of the matrix Q introduced in relation (3.17).
FGCRO-DR. Since AZ m P k = V k R (Proposition 2), we deduce from Lemma 2.3 This relation leads to the following result: ], using Y = W m P k , can be written as Since both V k+1 and V k+1 have orthonormal columns, we deduce from (3.27) and (3.30) that there exists a unitary matrix Q such that which proves the relation proposed in (3.17). Part IIb: Steps 6a-10a, 6b-10b. We show that at the end of steps 10a and 10b the following relation holds: range(Z k+1 ) = range( Z k+1 ). This result will help us to prove the existence of the matrix X introduced in relation (3.18).

Lemma 3.3 and Lemma 3.2, respectively, give us two useful relations for
Using the assumption on the equivalent preconditioning matrix obtained at the end of the ith cycle, i.e., M (i) The fact that identical (possibly nonlinear) preconditioning operators are used in steps 10a and 10b of Algorithm 1 (see relation (3.1)) allows us to write Relations (3.37) and (3.38) finally show the relation (3.34). Consequently from relations (3.32), (3.33), and (3.34) we deduce that there exists a nonsingular matrix X ∈ C (k+1)×(k+1) such that This proves the relation proposed in (3.18). Since T and T are both triangular, we note that X = T T −1 is also triangular.
Part IIIa: Steps 10a and 10b. We first show that v k+2 = v k+2 by expressing these two quantities as a function of r (i) 0 and Y . FGCRO-DR. The Arnoldi relation (step 10a) yields v k+2 =v k+2 /||v k+2 ||, ). Since from Proposition 2 V H k r (i) 0 = 0 in the (i + 1)th cycle, (I n − v k+1 v H k+1 ) and (I n − V k V H k ) commute, and from Part IIa of the proof, the following expression can be derived: ).
Part IIIb: Steps 10a and 10b. In this part we continue the analysis of the Arnoldi procedure with flexible preconditioning and show that v k+2+j = v k+2+j for j = 1, . . . , m − k − 1.
For the case j = 1, we introducev k+3 and v k+3 such that v k+3 =v k+3 /||v k+3 || and v k+3 = v k+3 /|| v k+3 ||. The application of the Arnoldi procedure in both algorithms leads tov  (3.20). This finally shows the main relations (3.9) and (3.10) of Theorem 3.6 that are satisfied at the end of the (i + 1)th cycle.  can be chosen equal. We deduce that FGCRO-DR and FGMRES-DR are algebraically equivalent.

About GCRO-DR and GMRES-DR.
We propose a second consequence of Theorem 3.6 analyzed now with a fixed preconditioning matrix M .
Corollary 3.8. When a fixed right preconditioner is used, the GCRO-DR and GMRES-DR methods sketched in Algorithm 1 are algebraically equivalent.
Proof. We denote by M the fixed right preconditioning operator. A straightforward reformulation of Lemma 3.3 in this context leads to the relation Z m = M W m in GCRO-DR. Exploiting now Lemma 2.3 allows us to derive the following relation, which holds during the (i + 1)th cycle: (3.42) and Part IIIa of the proof of Theorem 3.6 we deduce the following development:v By induction it is possible to deduce the rest of the proof regardingv k+j , j > 2. Using range( V k+1 ) = range(V k+1 ) obtained in Part IIa we deduce that Consequently the minimization problem min r (i) 0 − AZ m y leads to the same solution for both algorithms at each cycle: GCRO-DR and GMRES-DR sketched in Algorithm 1 are thus algebraically equivalent.

A numerical illustration.
In this section we intend to illustrate the results shown in sections 3.3.1 and 3.3.2 on a simple numerical example. We consider a real symmetric positive definite matrix A = Q D Q T of size 200 with Q orthonormal and D diagonal with entries ranging from 10 −4 to 1. The spectrum of A contains eigenvalues of small magnitude, 1 and consequently the use of deflation techniques should improve the convergence rate of Krylov subspace methods if the harmonic Ritz values of smallest modulus are taken into account. In this experiment we consider a polynomial preconditioner represented by two iterations of unpreconditioned GMRES for the solution of Ax = b with b given by b = Ae Ae 2 (e ∈ R 200 denoting the vector with all components equal to one) starting from a zero initial guess.  over a subspace of same dimension, i.e., FGMRES(10), FGMRES-DR(10,6), FGCRO-DR(10,6), respectively, and full flexible GMRES with such a variable preconditioner. Flexible methods with deflated restarting are found to be efficient since they are close to the full flexible GMRES method in terms of performances. We also remark that the convergence histories of FGCRO-DR(10,6) and FGMRES-DR(10,6) are different. According to Corollary 3.7 we compute the scalar product of v k+2 and v k+2 (which are both vectors of unit norm) to determine the cosine of the angle between these two vectors. The values are reported in Table 3.1 for the first five cycles. With such a variable preconditioner it is found that the methods are not equivalent in the first cycle already since the collinearity condition between v k+2 and v k+2 is not fulfilled. The situation is similar during the following cycles, which explains why different convergence histories for FGMRES-DR(10,6) and FGCRO-DR(10,6) observed in Figure  3.1 are obtained in such a case. As expected from section 3.3.2, if a fixed right preconditioner is used, the convergence histories of GMRES-DR(10,6) and GCRO-DR (10,6) are found to be exactly the same (results not shown here). In such a case v k+2 and v k+2 fulfill the collinearity condition; this is confirmed in Table 3.1 when a diagonal preconditioning is used. (m, k). In this section we first compare FGCRO-DR(m, k) with FGMRES-DR(m, k) presented in Algorithm 1 from both a computational and a storage point of view. Then we detail how subspace recycling can be used in FGCRO-DR(m, k) when solving a sequence of linear systems.

Computational cost.
We first analyze the computational cost related to the generalized eigenvalue problem to deduce harmonic Ritz information and then detail the total cost of the proposed method.

Harmonic Ritz information.
The generalized eigenvalue problem (3.22) can also be written as can be decomposed at the end of the cycle as where the structure of the (k + 1) × (k + 1) block V H k+1 W k+1 is as follows: V H k W k is a k × k matrix that satisfies the following relation at the end of the ith cycle: where the superscript is related to the cycle index. Thus storing the (m + 1) × m matrix (V H m+1 W m ) (i−1) can be used to obtain (V H k W k ) (i) at a cost that is independent of n. Next we analyze how to compute efficiently v H k+1 W k during the ith cycle. From relation (2.18) shown in Lemma 2.3 and Proposition 2, respectively, we deduce the relation

Due to Proposition 2 and the definition of
where β 1×k is obtained from relation (2.19), which does only involve projected quantities. This allows us to deduce v H k+1 W k at a cost independent of n. From this development we draw two important consequences from a computational point of view. First, (V H m+1 W m ) (i) can be obtained recursively at a cost that is independent of the problem size n. Second, storing W m (which would represent m additional vectors of size n) is not mandatory; only V H m+1 W m -matrix of size (m + 1) × m-is required.  , k) and FGCRO-DR(m, k). C represents the total cost of FGCRO-DR(m, k) and corresponds to C = (m − k)(op M + op A ) + n(2(m + 1)k + 1 + 2mk + (m − k)(2m + 2k + 6)).

Cost of a cycle.
We summarize in Table 4.1 the main computational costs associated with each generic cycle of FGMRES-DR(m, k) and FGCRO-DR(m, k). In FGCRO-DR(m, k), an Arnoldi method based on the modified Gram-Schmidt procedure has been assumed. 2 We have included only the costs proportional to the size of the original problem n which is supposed to be much greater than m and k. We denote by op A and op M the floating point operation counts for the matrix-vector product and the preconditioner application, respectively.
The generalized eigenvalue problem in FGCRO-DR(m, k) has been ignored in Table 4.1 since it can be performed at a cost independent of n as outlined in section 4.1.1. Furthermore the computation of c (required at step 12 of Algorithm 1) has not been included in Table 4.1 since in both methods it can be obtained at a cost independent of n (see Proposition 3 in [15] for FGMRES-DR). From Table 4.1 we deduce that FGCRO-DR(m, k) requires slightly fewer operations per cycle than FGMRES-DR(m, k).

Storage requirements.
We consider only the storage proportional to the size of the original problem n. Similarly, as in FGMRES-DR(m, k) (see [15, section 3.2.2]), if the matrix multiplications V m+1 Q and Z m P k R −1 at steps 8a and 9a of Algorithm 1 are performed in place (i.e., overwriting V k and Z k , respectively), FGCRO-DR(m, k) requires only the storage of Z m and V m+1 , which corresponds to (2m + 1) vectors of length n. The same storage cost is needed in FGMRES-DR(m, k) as detailed in [15].

Solution of sequence of linear systems.
As advocated in [22], GCRO-DR(m, k) is suited for the solution of a sequence of slowly changing linear systems defined as A l x l = b l where both the matrix A l ∈ C n×n and the right-hand side b l ∈ C n change from one system to the next, and the linear systems may typically not be available simultaneously. Next, we analyze how subspace recycling can be used in FGCRO-DR(m, k). We suppose that FGCRO-DR(m, k) has been applied for the solution of a given linear system (indexed by s−1) in this sequence and that appropriate subspaces to be recycled, Z s−1 k and W s−1 k , have been selected during a given cycle. ) are then supposed to hold. Proposition 3 details how to consider subspace recycling in the initial phase of FGCRO-DR(m, k), when solving the new linear system A s x s = b s with x 0 as an initial guess.
R −1 is imposed to make sure that the relation shown in Lemma 3.3 will hold at the end of the initial phase of FGCRO-DR(m, k) with subspace recycling.
In the case of a sequence where only the right-hand sides are changing, we note that the reduced QR factorization (step 3 in Algorithm 2) is not required. The complete construction of the initial generation of subspaces    Subspace recycling can thus be easily used in FGCRO-DR(m, k) to solve sequences of linear systems.

A numerical illustration.
As a numerical illustration we consider sequences of linear systems arising from the finite difference discretization of multidimensional elliptic partial differential equations (isotropic Laplace operator) posed on the [0, 1] d hypercube with homogeneous Dirichlet boundary conditions. These sequences correspond to situations where only the right-hand sides are changing for a given dimension d. An efficient solution method is of primary interest in certain applications related to, e.g., financial engineering, molecular biology, or quantum dynamics [5,6]. In the numerical experiments reported here (performed in MATLAB) we have used second order finite difference discretization schemes leading to sparse matrices with at most 2d + 1 nonzero elements per row. We analyze the performances of various flexible methods used with four iterations of unpreconditioned GMRES as a preconditioner. This polynomial preconditioner is a variable nonlinear function which thus requires a flexible Krylov subspace method as an outer method [28]. Table  4.2 collects the number of matrix-vector products of some flexible methods minimizing over a subspace of the same dimension, i.e., FGMRES (20), FGMRES-DR(20,10), FGCRO-DR (20,10), and FGCRO-DR(20,10) with subspace recycling, respectively. Using deflation helps to improve the convergence rate of FGMRES in this application since a reduction of approximately 20% to 25% in terms of matrix-vector products is obtained for FGMRES-DR(20,10) independently of the dimension d. FGCRO-DR(20,10) leads to numbers of matrix-vector products which are similar to FGMRES-DR(20,10) although the convergence histories are found to be different. Finally, using both deflation and recycling in FGCRO-DR leads to a significant decrease in terms of matrix-vector products. A reduction in the range of 40% to 45% is indeed obtained versus another flexible Krylov subspace method with deflated restarting (FGMRES-DR(m, k)). This can be considered as a primary advantage over FGMRES-DR(m, k) since FGMRES-DR(m, k) does not allow subspace recycling. It nicely extends to the flexible setting the advantage of GCRO-DR versus GMRES-DR previously illustrated in [22]. We note that the resulting method is factorization free and mostly relies on matrix-vector products, a nice feature if distributed memory platforms are targeted to address numerical problems of larger size in higher dimension.

Conclusion and perspectives.
In this paper we have studied a new minimum residual norm subspace method with deflated restarting that allows flexible preconditioning based on the GCRO subspace method. The resulting method, named FGCRO-DR, has been presented together with FGMRES-DR, a recently proposed algorithm of the same family but based on the GMRES subspace method. A theoretical comparison analysis of both algorithms has been performed in section 3, where Theorem 3.6-the main result of this paper-proves the algebraic equivalence of FGMRES-DR and FGCRO-DR if a certain collinearity condition holds at each cycle. Corollary 3.8 has also proved that GMRES-DR and GCRO-DR are algebraically equivalent when a fixed right preconditioner is used. Furthermore we have carefully analyzed the computational cost of a given cycle of FGCRO-DR and have shown that FGCRO-DR is nearly as expensive as FGMRES-DR in terms of operations. FGCRO-DR offers the additional advantage of being suitable for the solution of sequences of slowly changing linear systems (where both the matrix and right-hand side can change) through subspace recycling.
In [8] variants of FGCRO-DR have been proposed which only differ in the formulation of the projected generalized eigenvalue problem. In future work we plan to investigate the numerical properties of these variants on realistic problems of large size for both single and multiple left-or right-hand side situations. Of interest are applications related to, e.g., steady or unsteady simulations of nonlinear equations [7] or stochastic finite element methods [12,33] in three dimensions where variable preconditioning using approximate solvers has to be usually considered. We also note that when all right-hand sides are available simultaneously and when the matrix is fixed, block subspace methods may be also suitable. Thus a perspective could be to propose a block variant of FGCRO-DR.