domain: spi phys

text:
Mehdi Jadoui
email: mehdi.jadoui@onera.fr.

Christophe Blondeau
email: christophe.blondeau@onera.fr

François-Xavier Roux
email: roux@ann.jussieu.fr

Recycling Krylov Subspaces for Efficient Partitioned Solution of Aerostructural Adjoint Systems

Robust and efficient solvers for coupled-adjoint linear systems are crucial to successful aerostructural optimization. Monolithic and partitioned strategies can be applied. The monolithic approach is expected to offer better robustness and efficiency for strong fluid-structure interactions. However, it requires a high implementation cost and convergence may depend on appropriate scaling and initialization strategies. On the other hand, the modularity of the partitioned method enables a straightforward implementation while its convergence may require relaxation. In addition, a partitioned solver leads to a higher number of iterations to get the same level of convergence as the monolithic one. The objective of this paper is to accelerate the partitioned solver by considering techniques borrowed from Krylov subspace recycling strategies adapted to sequences of linear systems with varying right-hand sides. Indeed, in a partitioned framework, the structural source term attached to the fluid block of equations affects the right-hand side with the nice property of quickly converging to a constant value. We will also consider error approximation and approximate eigenvectors deflation in conjunction with advanced inner-outer Krylov solvers for the fluid block equations. We demonstrate the benefit of these techniques by computing the coupled derivatives for an aeroelastic configuration of the ONERA-M6 fixed wing in transonic flow. For this exercise the fluid grid was coupled to a structural model specifically designed to exhibit a high flexibility. All computations are performed using RANS flow modeling and a fully linearized one-equation Spalart-Allmaras turbulence model.

I. Introduction

N

, aeronautical industries are deeply concerned about the environmental issues generated by the growth in air traffic. As such, they are constantly confronted with new technological challenges, either by improving the performance of an aircraft, or by making new insights. One of them consists in finding, in the early design stage, the best potential design configuration that minimizes a certain objective function (e.g the lift-to-drag ratio) with respect to a set of design variables leading potentially in the end to a great reduction of CO 2 emission. Fine control of an aerodynamic shape or of a structural layout leads to a high-dimensional parameter space. For high-fidelity simulations, gradient-based optimizers in conjunction with the adjoint approach are the methods of choice. However, the coupled adjoint linear system is inherently ill-conditioned as it embeds matrix blocks of different scales and structures. In addition, the fluid block coming from the exact linearization of the RANS equations associated to a turbulence model is often very stiff. Besides, a strong level of fluid-structure interaction is known to be detrimental to the robustness and efficiency of existing solution techniques. Resolving such linear systems is achieved by two main approaches. The partitioned (or segregated) approach and the monolithic approach. The first one simply consists in resolving, in an alternating way, the aerodynamic and the structural sub-problems by applying the Linear Block Gauss-Seidel algorithm (LBGS). It expresses the inder-disciplinary coupling as a source term to the right-hand side of each set of disciplinary adjoint equations. The modularity of this approach makes it rather interesting since it takes advantage of the specific routines designed for each sub-problem and does not demand a high implementation cost. Nevertheless, this approach becomes rapidly inefficient and could even diverge for strong fluid-structure coupling even though the addition of some degree of relaxation helps to mitigate this issue. On the other hand, the monolithic approach consists in solving the fluid and structural equations simultaneously making it more robust in the sense that it is less sensitive to the strength of the fluid-structure interaction. The coupled adjoint system is generally solved by using Krylov subspace methods. The challenging aspect of such an approach is then to develop advanced preconditioning strategies combined with numerical ingredients so that Krylov methods reach the best performances in terms of robustness and efficiency.

In a high-fidelity aerostructural optimization context, Zhang and Zingg [START_REF] Zhang | Efficient monolithic solution algorithm for high-fidelity analysis and optimization[END_REF] implemented a robust monolithic solution method for both aerostructural analysis and coupled adjoint problem. A three-dimensional formulation was adopted involving the mesh, the flow and the structural states. The performance of the monolithic method as well as the partitioned one was investigated through a comparative study by varying the level of fluid-structure coupling. For the coupled adjoint problem solution, a GCROT [START_REF] De Sturler | Truncation strategies for optimal Krylov subspace methods[END_REF] Krylov solver has been used in conjunction with a block Gauss-Seidel preconditioner. The monolithic adjoint solution has been 60 % more efficient than the partitioned one for strong coupling. For weak fluid-structure coupling, the monolithic solution still outperformed the partitioned one with a better efficiency of 40 %. In terms of computational time, the monolithic method showed 50 % to over 60 % faster than the partitioned method. A similar comparative study of both monolithic and partitioned approaches was performed by Kenway et al. [START_REF] Kenway | Scalable parallel approach for high-fidelity steady-state aeroelastic analysis and adjoint derivative computations[END_REF] except that a Block Jacobi preconditioner was applied to the coupled adjoint system. The aerodynamic and structural block preconditioners were solved by a preconditioned Krylov method (restarted Generalized Minimal RESidual -GMRES [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]) and a direct factorization method respectively. A Flexible Krylov method (e.g FGMRES [START_REF] Saad | A flexible inner-outer preconditioned GMRES algorithm[END_REF]) has been used for the coupled adjoint system solution. The Common Research Model (CRM) wing-body-tail configuration was sized by considering two critical load cases: 1g cruise condition with moderate elastic deformation and a 2.5g pull-up with significantly more deflection. For the same memory footprint, the best monolithic solution seems to outperform the best partitioned one by reducing the time to 19% for the 1g load and by 29% for the 2.5g case. These numerical experiments demonstrate the great benefit of using monolithic approach in the strong coupling case but at the price of a robust preconditioner for the Krylov solver. We note however that both studies only considered inviscid flow modeling. These conclusions about the monolithic solver efficiency might be mitigated by the added stiffness of adjoint system matrices produced by a RANS fluid model associated to a linearized turbulence model.

Although the satisfactory performance of monolithic solvers, advanced strategies that could accelerate the partitioned algorithm using black-box solvers have received less attention. We recall that a partitioned algorithm consists in approximately solving the aerodynamic adjoint block at each fluid-structure iteration resulting in a sequence of adjoint linear systems with varying right-hand sides. As already mentioned, the structural source term that affects the right-hand side of the fluid block has the nice property of rapidly converging to a constant value. The corollary of this property is that after several fluid-structure couplings, the subsequent fluid systems should greatly benefit from recycling spectral information from the previous fluid-structure cycles. At the start of this work, the current partitioned adjoint solver did not take advantage of recycling and at each update of the structural source term the Krylov solver did a cold start from the previous solution.

The principle of deflation is to remove the influence of a system's subspace on the iterative process. This is usually beneficial when directions of certain subspaces hamper convergence. Deflation of an eigenspace can be performed in two ways: the linear system (matrix and right-hand side) is left-multiplied by a projector P, i.e. deflation by projection [START_REF] Vuik | An efficient preconditioned CG method for the solution of a class of layered problems with extreme contrasts in the coefficients[END_REF][START_REF] Tang | On deflation and singular symmetric positive semi-definite matrices[END_REF][START_REF] Gaul | Recycling Krylov subspace methods for sequences of linear systems : Analysis and applications[END_REF], or some eigenvectors are added to the Krylov subspace, i.e. deflation by augmentation [START_REF] Morgan | GMRES with deflated restarting[END_REF]. A survey of deflation and augmentations techniques can be found in [START_REF] Coulaud | Deflation and augmentation techniques in Krylov linear solvers[END_REF]. A specific type of deflation preconditioning aims at solving a rank-deficient projected system, using a Krylov solver, in a certain subspace outside of the problematic subspace. The solution is then complemented with the solution in the latter subspace. In [START_REF] Gaul | Recycling Krylov subspace methods for sequences of linear systems : Analysis and applications[END_REF] the author makes the post-correction superfluous by using a projection as right-preconditioner instead. For deflation by augmentation, adding eigenvectors to the Krylov subspace can effectively deflate corresponding eigenvalues from the spectrum because when these directions are included in the solution approximation, the convergence of the Krylov solver continues according to the modified spectrum. The deflation by augmentation led to the well-known FGMRES-DR solver [START_REF] Giraud | Flexible GMRES with deflated restarting[END_REF] and its extension to inner-outer Krylov solver [START_REF] Simoncini | Flexible inner-outer Krylov subspace methods[END_REF][START_REF] Jadoui | Comparative Study of Inner-Outer Krylov Solvers for Linear Systems in Structured and High-Order Unstructured CFD Problems[END_REF].

Unfortunately, the deflated restarted GMRES framework based on subspace augmentation is not adapted for solving sequences of linear systems [START_REF] Parks | Recycling Krylov subspaces for sequences of linear systems[END_REF]. Fortunately, some authors have proposed new strategies in order to reuse information accumulated in previous fluid-structure cycles and use it to accelerate the solution of the next linear system. Krylov subspaces recycling methods seem to be the suitable choice. Historically, De Sturler suggested the Generalized Conjugate Residuals with inner Orthogonalization (GCRO) method [START_REF] De Sturler | Nested Krylov methods based on GCR[END_REF], an improved version of the recursive GMRES (GMRESR) solver [START_REF] Van Der Vorst | GMRESR: a family of nested GMRES methods[END_REF] by maintaining an orthogonality condition between the outer and the inner spaces generated by GMRESR. This way, it provides the optimal correction to the solution in a global search space. Later, Parks et al. formulated the GCRO-DR algorithm [START_REF] Parks | Recycling Krylov subspaces for sequences of linear systems[END_REF] that combines GCRO and deflation techniques by augmentation introduced by Morgan [START_REF] Morgan | GMRES with deflated restarting[END_REF]. They demonstrated better performances of GCRO-DR compared to the GMRES-DR in a long sequence of linear systems from a fracture mechanics problem. Carvhalo et al. extended GCRO-DR to the flexible case (FGCRO-DR) [START_REF] Carvalho | A flexible generalized conjugate residual method with inner orthogonalization and deflated restarting[END_REF] and they conducted in-depth analysis of both flexible methods. In particular, they showed that both methods can be algebraically equivalent if a certain colinearity condition is satisfied at each cycle. In 2013, Niu et al. introduced Loose GCRO-DR (LGCRO-DR) [START_REF] Niu | Accelerated GCRO-DR Method for Solving Sequences of Systems of Linear Equations[END_REF] for improving the convergence of GCRO-DR by recycling both spectral information and approximate error information. The error is defined as the distance between the current iterate and the exact solution of the system. It is not known by definition but a fair approximation to it can be computed. This idea was initially proposed by Baker et al. [START_REF] Baker | A technique for accelerating the convergence of restarted GMRES[END_REF] and mimics the idea behind GMRESR of including approximations to the error in the current approximation space. This error information is interesting since it represents in some sense the previous Krylov space generated in the previous cycle and subsequently discarded. In addition to that, LGCRO-DR is straightforward and economic to implement.

In this paper, we investigate advanced Krylov subspace methods using subspace recycling strategies for accelerating the partitioned solver. More specifically, we compare the performances of flexible GCRO-DR with and without subspace recycling and flexible GMRES-DR for varying degrees of fluid-structure coupling. This work will benefit from the recent achievements to improve efficiency of the fluid adjoint solution by applying nested Krylov subspace methods [START_REF] Jadoui | Comparative Study of Inner-Outer Krylov Solvers for Linear Systems in Structured and High-Order Unstructured CFD Problems[END_REF].

The numerical experiments are performed on an aeroelastic configuration of the ONERA M6 fixed wing in transonic viscous flow.

This paper is organized as follows. In section II we briefly recall the theoretical background of aeroelastic and coupled-adjoint equations. The partitioned algorithm is also outlined. Section III is devoted to the description of the GCRO algorithm with adaptations related to variable preconditioning and subspace recycling. The aeroelastic numerical test case is presented in section IV. In section V we detail our preliminary numerical experiments. First, we show that FGMRES-DR and FGCRO-DR are numerically equivalent for solving the fluid adjoint system as demonstrated in [START_REF] Carvalho | A flexible generalized conjugate residual method with inner orthogonalization and deflated restarting[END_REF]. Then, the benefit of adding more physics in the fluid preconditioner is demonstrated for the standard coupled-adjoint partitioned solver. Last, we present preliminary results illustrating the acceleration of the partitioned solver thanks to subspace recycling.

II. Aerostructural adjoint system

A. Aeroelastic equilibrium

Let us denote the state variables of the coupled system W and U, representing the fluid conservative variables and the structural displacements respectively. At the aeroelastic equilibrium, the state variables and the meshes satisfy the discretized equations of fluid and structural mechanics simultaneously:

R a (X a , W, U) = 0 R s (X s , W, U) = 0 (1)
where R a is the discrete aerodynamic residual and R s the discrete structural residual. These two blocks of equations are coupled through aerodynamic forces Q a loading the skin of the structure and the structural displacements U deforming the fluid mesh. The structural mesh is noted X s . In the following we introduce two aerodynamic grids X a and X a0 . X a denotes the deformed aerodynamic grid at the aeroelastic equilibrium at the outcome of the aeroelastic analysis. X a0 is called the reference mesh which supports the aerodynamic shape parametrization. Typically for an aircraft design study the reference mesh is chosen as the jig shape or the flight shape in reference nominal flight conditions. The load, displacement and mesh deformation operators then merely depend on X a0 for an aeroelastic or coupled-adjoint analysis.

The structural loads Q s are obtained with a suitable load transfer technique applied to Q a such that

Q s (Q a (W, X a ), X a0 , X s ) = T Q sur f (X a0 , X s )Q a (W, X a ) (2)
where T Q sur f represents a linear load transfer operator. The subscript sur f stipulates that the associated linear operators or data relate to the fluid-structure interface. The structural displacements alter the fluid grid locations through the relation:

X a = X a0 + δX a (δX a,sur f , X a0 ) = X a0 + T vol (X a0 )δX a,sur f (3)
with T vol (X a0 ) the volume operator performing the deformation of the fluid domain. The vector δX a,sur f corresponds to the displacements of the fluid nodes at the fluid-structure interface.

δX a,sur f = δX a,sur f (X a0 , X s , U) = T U sur f (X a0 , X s )U (4)
where T U sur f (X a0 , X s ) represents a linear displacement transfer operator.

B. Partitioned strategy for the coupled adjoint system

Let us consider a scalar aeroelastic objective function J(W, X a , X s ) and a design parameter p. One way to obtain the coupled adjoint equations is to formulate an augmented objective function by adding the total variation of the residuals R s and R a to the total derivative dJ/dp. More specifically, we define d J/dp as

d J dp = dJ dp + λ T a dR a dp + λ T s dR s dp (5)
where

dJ dp = ∂ J ∂W dW dp + ∂ J ∂X a dX a dp + ∂ J ∂X s dX s dp + ∂ J ∂U ∂U ∂p , (6)
In Eq. ( 5) the total variations of residuals are exactly zero since they represent constraints related to the satisfaction of the equilibrium equations at the outcome of the aeroelastic analysis. For simplicity, we restrict here to the specific case of a shape design parameter not affecting the structural geometry nor the structural stiffness. In addition, the explicit dependency of the objective function with respect to the structural states is dropped, i.e., we consider only derivatives of aerodynamic coefficients. The full derivation for the general case can be found in [START_REF] Achard | Techniques de calcul de gradient aéro-structure haute-fidélité pour l'optimisation des voilures flexibles[END_REF][START_REF] Achard | High-Fidelity Aerostructural Computation Techniques with Application to a Realistic Wing Sizing[END_REF]. It is worth to mention that these assumptions do not lead to any loss of generality of the work presented in this paper since we focus on solution techniques for the adjoint system. As X s does not depend on the design parameter p, we have dX s /dp = 0. Under the same assumption we also have dK/dp = 0. After some algebra manipulation we end up with the following expression for d J/dp in which the total derivatives dW/dp and dU/dp have been factored out:

d J dp = ∂ J ∂W + λ T a ∂R a ∂W -λ T s C dW dp + ∂ J ∂X a A + λ T a ∂R a ∂X a A + λ T s (K -D) dU dp + ∂ J ∂X a + λ T a ∂R a ∂X a B dX a0 dp -λ T s E dX a0 dp (7)
Constant matrices A to E are defined analytically with the following formulas (see [START_REF] Achard | Techniques de calcul de gradient aéro-structure haute-fidélité pour l'optimisation des voilures flexibles[END_REF]):

A = T vol T U sur f (8)
B = ∂X a ∂X a0 = I + ∂A ∂X a0 U (9)
C = T Q sur f ∂Q a ∂W (10)
D = T Q sur f ∂Q a ∂X a T U sur f (11)
E = T Q sur f ∂Q a ∂X a B + ∂Q s ∂X a0 (12)
The coupled adjoint linear system is obtained by canceling factors related to dW/dp and dU/dp in Eq. ( 7) to give

         ∂R a ∂W T -C T A T ∂R a ∂X a T K T -D T                 λ a λ s        =         -∂J ∂W T -A T ∂J ∂X a T         (13)
The process for solving the adjoint system follows an iterative block scheme. Algorithm 1 details the Linear Block Gauss Seidel (LBGS) scheme applied for the solution of system [START_REF] Jadoui | Comparative Study of Inner-Outer Krylov Solvers for Linear Systems in Structured and High-Order Unstructured CFD Problems[END_REF]. In this derivation, we use the structural flexibility S which is a small reduced matrix relating the set of structural forces to the set of structural displacments pertaining to the fluid-structure coupling. The relaxation factor θ s has been introduced on the adjoint vector λ s . Assuming that the coupled system is solved to machine accuracy, the total derivative reconstruction is given by

dJ dp = ∂ J ∂X a + λ T a ∂R a ∂X a B dX a0 dp -λ T s E dX a0 dp (14)
In the expression above the computation of the product of the geometrical sensitivities with the matrix B is not trivial. If one already has at hand a linearized version of the operator A, i.e. of T U sur f and T vol , it can be applied to dX a0 /dp as many times as the number of design variables. This is the most straightforward manner but the benefit of the adjoint formulation is then mitigated by the cost of the gradient assembly. The other way is to transpose the first term in the right-hand side of Eq. ( 14) and compute products like [∂T U sur f /∂X a0 ] T v and [∂T vol /∂X a0 ] T v, where v has fluid grid size. We call this mode the geometrical adjoint of the mesh deformation and displacement transfer operators. These two modes of gradient assembly have been implemented in the coupled-adjoint module of the elsA software.

Algorithm 1 Partitioned LBGS strategy for coupled-adjoint solution

1: Initialize: U, W X a , X s , X a0 , λ 0 a , λ 0 s , θ s , tol 2: if λ 0 s 0 then 3: RHS stru ← T Q sur f ∂Q a ∂W T S T λ 0
A Xs,sur f ← ∂Q a ∂X a, sur f T T Q sur f T S T λ k-1 s Structural geometric adjoint 12: A X a ← -[λ k-1 a ] T ∂R a ∂X a + ∂J ∂X a T Aerodynamic geometric adjoint 13:
A X a,sur f ← (T vol ) T A X a Mesh deformation adjoint 14:

λ k s ← (T U sur f ) T (A Xs,sur f + A X a,sur f ) Structural adjoint vector 15: λ k s ← (1 -θ s )λ k-1 s + θ s λ k s Relaxation 16: RHS stru ← T Q sur f ∂Q a ∂W T S T λ k

III. Recycling strategies for sequences of linear systems

As already explained, the current implementation of our partitioned solver does not take advantage of spectral information produced during the Krylov solver instances applied to the sequence of preceding right-hand sides. An advanced inner-outer GMRES-DR solver is used for the approximate solution of the fluid block between consecutive fluid-structure couplings but a simple restart is performed when the structural source term is updated. Unfortunately, deflated restarting based on subspace augmentation by appending approximate Ritz vectors to the Krylov subspace is not suitable for sequences of linear systems. Indeed, if the Ritz vectors are obtained from a previous linear system with another matrix or even another right-hand side, the concatenation of the recycled subspace to the current Krylov subspace does not form a Krylov subspace for the current problem. It is necessary to introduce a new Krylov solver that uses recycling of any given subspace without restriction. A famous one is the generalized conjugate residual with inner orthogonalization (GCRO). It belongs to the family of inner-outer methods [START_REF] Axelsson | Iterative solution methods[END_REF] where the outer method is based on GCR, a minimum residual norm method proposed by Eisenstat, Elman and Schultz [START_REF] Eisenstat | Variational iterative methods for nonsymmetric systems of linear equations[END_REF]. The inner solver is GMRES applied to a projected system matrix.

A. Flexible Generalized Conjugate Residual with inner Orthogonalization and Deflated Restarting: FGCRO-DR

As previously mentioned, we focus on the solution of a sequence of linear systems with a varying right-hand side only. This section briefly introduces the Flexible GCRO algorithm with deflated restarting. In this framework, deflation can reuse spectral information from a previous cycle or from a previous linear system. We start by recalling the original formulation of the generalized conjugate residual method (GCR) [START_REF] Eisenstat | Variational iterative methods for nonsymmetric systems of linear equations[END_REF]. The idea is to introduce the concept of optimality for the solution residual. We want to solve

Ax (i) = b (i) , i = 1, 2, ... (15)
where A ∈ R n×n and b (i) ∈ R n changes from one system to the next. The GCR method relies on a given full-rank matrix Z k ∈ R n×k and an orthonormal matrix C k ∈ R n×k as the image of Z k by A satisfying the relations

AZ k = C k , (16)
C T k C k = I k . ( 17
)
For the sake of understanding, suppose that we want to solve the first system of the sequence of linear systems, i.e., i = 1. Given an initial guess x 1 0 , the principle is to compute an approximation to the solution x 1 ∈ x 1 0 + range(Z k ) that minimizes the corresponding residual norm in the approximation space range(Z k ). More precisely, GCR solves the following minimization problem

x (1) = argmin

x ∈x

(1)

0 +range(Z k ) b (1) -Ax 2 , (18)
The optimal solution of (18) over the subspace x (1) 0 + range(Z k ) is defined by

x (1) = x (1) 0 + Z k C T k r (1) 0 ( 19
)
Consequently, the corresponding residual vector satisfies

r (1) k = b (1) -Ax (1) = (I -AZ k C T k )r (1) 0 = (I -C k C T k )r (1) 0 , r (1) k ⊥ range(C k ). ( 20
)
The orthogonality of the residual r (1) k to the subspace AK k (A, r (1) 0 ) spanned by the columns of C k is known as the optimality property of the residual. In practice, GCR is not considered as a means to solve the linear system. It is replaced by a GMRES solver which rather computes an implicit representation of the matrices Z k and C k [START_REF] De Sturler | Truncation strategies for optimal Krylov subspace methods[END_REF].

Given an orthonormal basis C k of an outer subspace and the corresponding residual r k = (I -C k C T k )r 0 after k steps of GCR, GCRO(m) obtains the next iterates r k+1 and x k+1 by performing m steps of GMRES applied to the projected operator A C k = (I -C k C T k )A, thereby maintaining optimality of the inner residual to the outer space. Let V m+1 ∈ R n×(m+1) be an orthonormal basis for K m+1 (A C k , r k ) with v 1 = r k / r k , if we consider a flexible inner GMRES we have the following Arnoldi relation:

A C k Z m = (I -C k C T k )AZ m = V m+1 Hm with Hm = V T m+1 A C k Z m (21)
This equation can be expanded as

AZ m = C k B m + V m+1 Hm with B m = C T k AZ m (22)
At the end of the inner GMRES cycle, the residual is expressed as

r k+1 = r k -A C k Z m y = r k -V m+1 Hm y m (23)
with, using Hm = Qm R m , y m such that

y m = argmin y ∈R m r k -A C k Z m y 2 = R -1 m QT m r k 2 e 1 ( 24
)
In the above Arnoldi process [START_REF] Axelsson | Iterative solution methods[END_REF] of the inner GMRES, the vectors Az i are first orthogonalized against C k , thus constructing V m+1 such that C T k V m+1 = 0. So far, we have not yet explained how the matrix C k is built in the context of subspace recycling. First, we slightly change the notation in ( 21) and [START_REF] Axelsson | Iterative solution methods[END_REF] to point out that we perform mk steps (instead of m) of inner FGMRES. This produces the Arnoldi relation

A C k Z m = (I -C k C T k )AZ m-k = V m-k+1 Hm-k (25)
Combining ( 16) and [START_REF] Cambier | Verification and validation of the Onera elsA flow solver on RANS benchmarks[END_REF], where C k is renamed V k for better readability, we have

A[Z k Z m-k ] = [V k V m-k+1 ] I k B m-k 0 Hm-k with B m-k = V T k AZ m-k (26)
From now on, we define

Z m = [Z k Z m-k ] and V m+1 = [V k V m-k+1
] leading to the following compact expression of ( 26)

AZ m = V m+1 Hm . ( 27
)
We have now completed the process for the first linear system which formulates the solution from the complementary subspaces V k and V m-k+1 . At this stage, we introduce deflated restarting within the FGCRO framework to build the outer space V k for the next linear system. The deflation subspace is spanned by the k harmonic Ritz vectors corresponding to harmonic Ritz values of smallest magnitude. This leads to the FGCRO-DR algorithm that is a combination of FGCRO and FGMRES-DR (see Algorithm 2). To compute the harmonic Ritz vectors, we solve the m × m generalized eigenvalue problem defined in step 28 of Algorithm 2. We then form the matrix P k ∈ R m×k which stores the retained eigenvectors. This matrix is expanded by Z m to form the basis of the deflation subspace Z m P k . Then, the image of the deflation space by A gives

AZ m P k = V m+1 Hm P k (28)
then orthonormalized,

V k = V m+1 Qm ( 29
)
where V m+1 has been computed at the end of the (i -1)st cycle and Qm stands for the unitary factor of the QR factorization of Hm P k = Qm R k . The relation that holds at this stage is:

AZ m P k R -1 k = V k (30)
From ( 16), we deduce that Z k = Z m P k R -1 k . Knowing V k and Z k , the solution and the residual can be updated with

x (i+1) 1 = x (i+1) 0 + Z k V T k r (i+1) 0 (31) r (i+1) 1 = (I -V k V T k )r (i+1) 0 ( 32
)
and we can initiate the Arnoldi process in (26) for the next system.

IV. ONERA-M6 wing aeroelastic analysis

In this study the numerical experiments have been performed with the well known ONERA-M6 fixed wing configuration which has been extensively used for CFD solvers validation in transonic flow conditions. In this work we use the RANS solver provided by elsA for the steady rigid and aeroelastic analyses [START_REF] Cambier | The Onera elsA CFD software: input from research and feedback from industry[END_REF][START_REF] Cambier | Verification and validation of the Onera elsA flow solver on RANS benchmarks[END_REF]. The elsA adjoint and coupled-adjoint solvers have also been used for the computation of rigid and flexible derivatives. The latest improvements to the Krylov solvers for the solution of the adjoint linear system are described in [START_REF] Jadoui | Comparative Study of Inner-Outer Krylov Solvers for Linear Systems in Structured and High-Order Unstructured CFD Problems[END_REF]. A multi-block structured mesh featuring a C-H topology is used (Fig. 1). It consists of 3.8 million grids divided into 42 blocks. The flight conditions are a free-stream Mach number of 0.84 at an angle of attack of 3.06 degrees. The convective fluxes are discretized by an upwind Roe scheme associated to a Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL) reconstruction and a Van Albada flux limiter. The one-equation Spalart-Allmaras turbulence model is selected. The surface contours in the bottom plot of Fig. 2 below show typical results for the ONERA-M6 wing. The pressure coefficient contours identify a lambda-shock along the mid-chord of the wing. For the aeroelastic analysis a simple but realistic finite element model has been designed (Fig. 1). The stiffness of this model can be easily tuned to get stronger or weaker fluid-structure interaction. The pressure coefficient contours at the aeroelastic equilibrium are plotted in the upper part of Fig. 2 and can be compared to the rigid contours. The maximum vertical displacement is 0.14 meters corresponding to 11.7 % of the wing span. To get a better insight of the effect of flexibility on the pressure distribution, we report in Fig. 3 the Cp distributions for two sections at y = 0.60 m and y = 1.12 m. The vertical displacement distributions associated to the front and rear spars as well as the twist increment distribution are plotted in Fig. 4. The rigid analysis results in a lift coefficient C L = 0.27 whereas the aeroelastic analysis, at the same angle of attack, results in a lower lift coefficient C L = 0.23.

V. Numerical experiments A. Preconditioning

We are interested in the solution of linear systems of adjoint equations for a finite volume discretization on multi-block meshes. As already mentioned, all computations are performed using a RANS flow model and a fully linearized one-equation Spalart-Allmaras turbulence model. The spatial discretization relies on a second-order upwind Roe scheme. The associated system of equations is usually very stiff as it exhibits terms of very different amplitudes, depending on the location in the fluid domain, e.g., in the flow wake, near the wall surface or at the far boundary. Some recent work by the authors [START_REF] Jadoui | Comparative Study of Inner-Outer Krylov Solvers for Linear Systems in Structured and High-Order Unstructured CFD Problems[END_REF] demonstrated the numerical robustness of the inner-outer GMRES solvers for such difficult problems. This class of solvers embeds two level of preconditioning. The outer preconditioner is a standard GMRES solver pre-conditioned itself by an inner stationary operator. The GMRES preconditioner offers the desirable property of addressing the global eigenspectrum of the entire domain. This also makes the inner-outer solver less sensitive to the quality of the second-level preconditioner.

Two stationary preconditioning strategies are considered for the inner GMRES. The first one consists in a block version of a standard block Lower-Upper Symmetric Gauss-Seidel (LU-SGS) iterative solver [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. LU-SGS is a particular case of LU-SSOR with a unity relaxation factor. The reader should refer to [START_REF] Peter | Large stencil viscous flux linearization for the simulation of 3D compressible turbulent flows with backward-Euler schemes[END_REF] for further details about LU-SGS as an iterative solver for the implicit phase of Backward-Euler Newton steady solvers. In the current implementation, LU-SGS is combined with a domain decomposition method. More specifically, the globalization of the preconditioner is achieved with a Restricted Additive Schwarz method [START_REF] Cai | A Restricted Additive Schwarz Preconditioner for General Sparse Linear Systems[END_REF].

LU-SGS is applied to a first-order diagonally dominant upwind approximation of the flux Jacobian matrix inspired by [START_REF] Yoon | Lower-upper symmetric-Gauss-Seidel method for the Euler and Navier-Stokes equations[END_REF]. This operator is based on a first-order spatial discretization of the convective and of the viscous fluxes using a simplifying thin layer assumption [START_REF] Peter | Large stencil viscous flux linearization for the simulation of 3D compressible turbulent flows with backward-Euler schemes[END_REF]. This strategy leads to a very compact stencil for the preconditioning matrix which will be denoted as J AP P O1 in the remainder of this paper. The second one is a Block Incomplete Lower Upper (BILU(k)) factorization applied to either an approximate or exact flux Jacobian matrix. For the so-called first-order exact Jacobian matrix J E X O1 , a first-order spatial Roe scheme is used for the discretization of the mean-flow convective fluxes and a 5-point corrected centered discretization scheme is used for the diffusive fluxes. The BILU(k) preconditioner will be applied either to the first-order approximate Jacobian matrix J AP P O1 , or to the first-order exact Jacobian matrix J E X O1 . About memory footprint, J E X O1 has a 9-point stencil in 2D whereas a 5-point stencil is associated with J AP P O1 . In 3D, we have a 7-point stencil for J AP P O1 and a stencil of 19 points for J E X O1 . Consequently, a better robustness is achieved but at the price of about twice the storage for J E X O1 compared to J AP P O1 . For a better understanding, we have reproduced the 3D stencils in Fig. 5. In order to tackle the different scales in the system Jacobian matrix, an additional level of preconditioning is considered in the present work. We perform a block-diagonal left-right scaling, i.e., row-column scaling, of the system matrix leading to a modified system of the form

D -1 L AD -1 U M -1 t = D -1 L b, (33)
where D -1 L and D -1 U are lower and upper inverses of the LU factors of the block-diagonal of matrix A. This leads to a scaled system matrix with unit diagonal. The preconditioner is then based on an approximate BILU factorization of

D -1 L AD -1 U .
Even if the scaling is performed externally, i.e., the approximate Jacobian matrix is left/right pre-multiplied before applying BILU, the condition number of the BILU factors is greatly reduced and the rounding errors during the factorization process are smaller. About products involving the system matrix A, the operation is performed on-the-fly as the matrix is not stored.

B. Performance of recycling Krylov solvers for aerodynamic adjoint systems

Having coded our new FGCRO-DR solver, the primal objective in this section is to assess its performances and compare to FGMRES-FR for the solution of a fluid adjoint problem. This exercise will also confirm the expected algebraic equivalence that exists between FGCRO-DR and FGMRES-DR and proved by Carvalho et al. in [START_REF] Carvalho | A flexible generalized conjugate residual method with inner orthogonalization and deflated restarting[END_REF]. Since we solve a single linear adjoint system, the subspace recycling step (step 2 to step 4 of Algorithm 2 in Appendix A) is skipped. Actually, Fig. 6 illustrates an essential numerical aspect pertaining to the FGCRO-DR method. More precisely, FGCRO-DR seems to be more sensitive than FGMRES-DR to the rounding errors during the orthogonalization process. Indeed, for FGCRO-DR we observe a stagnation of the relative residual after 10 -8 when a single step of orthogonalization is applied, even for a moderate size Krylov subspace (FGCRO-DR(30,10,10) with nMGS=1).

The best way to restore convergence in a reasonable number of iterations is to perform two steps of orthogonalization (FGCRO-DR(30,10,10) with nMGS=2). The main reason that FGCRO-DR suffers from rounding errors during the orthogonalization process comes essentially from the two-step Arnoldi process. Indeed, maintaining the orthogonalization condition between the recycled space V k and the inner Krylov space V m-k during the Arnoldi process is challenging for our stiff numerical problem. To point out the two-step orthogonalization process, (25) can be rewritten as

AZ m-k = V k B m-k + V m-k Hm-k (34)
In Fig. 6 we see that the theoretical numerical equivalence, in exact arithmetics, between FGCRO-DR and FGMRES-DR is demonstrated. As we are dealing with RANS flow equations associated to a fully linearized turbulence model, the Jacobian operators may be poorly conditioned in certain areas of the fluid domain. To mitigate this issue, we have considered a block diagonal scaling of the system matrix as defined in (33). Fig. 7 illustrates the benefit of this scaling. Indeed, scaling has a strong impact on the stagnation observed at the beginning of the convergence by shrinking the plateau.

Also, true and least-squares relative residuals match more accurately even at tiny convergence levels, i.e., 10 -9 in this case. Scaling seems then to be a good alternative to re-orthogonalization, and also a cheaper one. In Fig. 7 we notice a slightly lower convergence rate for the scaled problem so that both scaled and original problems converge in approximately the same number of iterations. However, in the context of a partitioned solution strategy, we seek for an early decrease of the fluid residual and, in this sense, scaling is very attractive.

C. Performance of the current partitioned solver for coupled-adjoint systems

This subsection illustrates the impact of preconditioning strategies applied to the fluid adjoint problem that improve the current partitioned solver in terms of robustness and computational time. We recall that an inner-outer Krylov solver is used in this work and then the selection of the preconditioner of the inner Krylov solver is of interest. The preconditioning strategies have been presented in section V.A.

Considering the relative size of the stencils associated to the exact first-order and the approximate first-order Jacobian matrices, one needs roughly three times more storage for J E X O1 compared to J AP P O1 but resulting in a better efficiency as we can see in Fig. 8. Indeed, BILU(0) applied to J E X O1 outperforms LU-SGS with a gain of 45% in terms of iterations. Moreover, the scaling has the ability to enhance the convergence of the coupled-adjoint system by shrinking the plateau that appears at the beginning of each fluid-structure coupling (see Fig. 9). From these numerical experiments, we conclude that the best option for performance of the partitioned solver is to associate a diagonal scaling to the BILU(0) factorization of the first-order exact Jacobian matrix. We also recall that the flow equations are systematically physically non-dimensionalized in the first place.

D. Performance of recycling strategies for aerostructural adjoint systems

The results presented in section are associated to an approximate linearization using the constant eddy viscosity assumption. However, we point out that the conclusions can be readily transposed to a fully linearized system of equations.

Our objective is to propose simple criteria to assess whether recycling can be activated. The first criterion consists in monitoring the norm of the right-hand side of the fluid block. Indeed, recycling is expected to be efficient for slowly varying systems. As we are solving a coupled problem, we know that after a given number of fluid-structure (F/S) couplings, the structural source term stabilizes. This stabilization of the structural part means that the aerostructural adjoint system solution reduces to an aerodynamic adjoint system only.

In Fig. 10, we plot the assembled right-hand side (red curve) which sums the structural one (RHS stru , blue curve) and the aerodynamic one (RHS aer o , green curve). We recall that RHS aer o is constant during the solution process and only RHS stru varies. We observe in Fig. 10 that RHS stru varies significantly during the first fluid-structure cycles and stabilizes after six couplings. Fig. 11 clearly shows the benefit of subspace recycling on the relative residual convergence of the partitioned solver. As just mentioned, the recycling strategy is relevant after six fluid-structure cycles. More precisely, we have considered two cases: the first one triggers the recycling strategy from the fifth cycle and the second one from the sixth cycle. Obviously, recycling is not relevant before six couplings.

Fig. 12 proposes another criterion for the activation of the recycling process. More precisely, at the beginning of the ith fluid-structure cycle, we compare the norms of the initial relative residual r (i) 0 and that of the relative residual projected orthogonally to the recycled subspace

V k = [v 1 , ..., v k ] denoted as r (i) 1 = (I -V k V T k )r (i) 0 .
The ratio r (i) 1 / r (i) 0 can be interpreted as a measure of the degree of colinearity between the two residual vectors. Thus, a value close to unity indicates that the previous approximation V k of the invariant subspace of the system matrix is relevant for recycling. By the way, as fluid-structure cycles proceed, we expect this approximation to get closer to the true invariant subspace. One interesting way to assess the approximation error is to monitor the distance between subspaces V (i) k and

V (i+1) k . This distance is expressed as dist(V (i) k , V (i+1) k ) = (I -V (i) k [V (i) k ] T )V (i+1) k [V (i+1) k
] T . Fig. 11 reflects these considerations. The convergence of the FGMRES-DR solver is reported for reference. As expected, it is directly comparable to the FGCRO-DR curve without recycling. Turning to the convergence of FGCRO-DR with recycling from the fifth F/S cycle, we observe a degradation of convergence. Clearly, the recycled space is not sufficiently informative to improve convergence during the subsequent cycles. However, recycling at the next cycle dramatically speeds up convergence.

E. Deflation strategies based on error approximation

As part of the GCRO-DR Niu et al. [START_REF] Niu | Accelerated GCRO-DR Method for Solving Sequences of Systems of Linear Equations[END_REF] proposed to recycle both the approximate error and the harmonic vectors corresponding to the smallest harmonic Ritz values in magnitude. Let first recall the definition of the error approximation. If we denote x s and x s-1 the approximate solutions obtained at the end of cycle (s) and cycle (s -1) of GMRES(m), and x the exact solution, then the approximate error is defined as z s = x sx s-1 . It represents the correction to x s-1 selected from K m (A, r s-1 ). The specific choice of z s for recycling is interesting since it represents the best available approximation of A -1 r s-1 . Indeed, De Sturler demonstrated that the true error vector e s = xx s is the optimal solution of the linear system Ae = r s-1 . Therefore, deflating the error approximation z s enables the GCRO-DR solver to compute a new search direction thus avoiding the z s -direction in the subsequent cycles. This strategy leads to a new solver which is denoted as LGCRO-DR(m,kp,p) where p stands for the number of error approximations that are recycled. Typically, p varies between 1 and 3. Note that we want to keep the size k of the recycled space unchanged, irrespective of the small number of recycled error, that is, p < k. Baker [START_REF] Baker | A technique for accelerating the convergence of restarted GMRES[END_REF] proposed a numerical criterion based on both sequential and skip angles to analyze the co-linearity of residuals between cycles. We recall a simple expression of the sequential and skip angles relying on the residual norm (see [START_REF] Niu | Accelerated GCRO-DR Method for Solving Sequences of Systems of Linear Equations[END_REF]). Let r s+1 and r s be the (s+1)st and the (s)th residual vectors from LGCRO-DR(m,kp,p). Then the sequential angle between two residual vectors is

cos∠(r s+1 , r s ) = r s+1 r s . (35)
Also, let r s+1 and r s-1 be the (s+1)st and the (s-1)st residual vectors from LGCRO-DR(m,kp,p). Then the skip angle between two residual vectors is

cos∠(r s+1 , r s-1 ) = r s+1 r s-1 . (36)
We point out that the expression of the skip angle (36) is only valid for the LGCRO-DR method. Otherwise, the skip angle is formulated for the non-flexible GCRO-DR as:

cos∠(r s+1 , r s-1 ) = r s+1 r s-1 - < Az s+1 , Az s > r s+1 r s-1 , (37)
where z s+1 = x s+1x s ∈ K m (A, r s ) stands for the error approximation of the (s+1)st cycle. Fig. 13 shows how the recycling of the approximation improves the robustness of the GCRO-DR solver. In order to highlight the interest of the error approximation strategy, we decided to solve the aerodynamic adjoint system for two different sizes of the Krylov space. We consider both GCRO-DR(60,20) and GCRO-DR(30,10). As we can see in Fig. 13, GCRO-DR(60,20) converges after 1900 iterations. Even though a plateau is observed at the beginning, the convergence of GCRO-DR(60,20) seems to follow the superlinear convergence property of GMRES without restarting. Fig. 14a depicts the variation of the sequential and skip angles for GCRO-DR(60,20) throughout the cycles. At the beginning, we observe for both sequential and skip angles low values which corroborates the plateau that we see in Fig. In contrast, sequential and skip angles are gradually improved throughout the cycles. More precisely, the skip angle reaches up to 80 degrees whereas the sequential angle reaches up to 60 degrees. The second test case consists in reducing the size of the Krylov space by considering GCRO-DR(30,10). In that case, the GCRO-DR solver stagnates. Indeed, this stagnation is clearly illustrated in Fig. 14b where the sequential and skip angles never exceed 2 degrees. Let see how the recycling of error approximations affects the convergence behaviour. More specifically, a maximum of two error approximations have been recycled. The red curve in Fig. 13 shows numerically the capability of one error approximation to retrieve the convergence after 9000 iterations. We point out that we have computed the error approximation in an economical way (see [START_REF] Niu | Accelerated GCRO-DR Method for Solving Sequences of Systems of Linear Equations[END_REF]). This economical method has the advantage of straightforward and easy to implement but is inexact. The LGCRO-DR(30,8,2) solver outperforms the LGCRO-DR(30,9,1) and improves even more the solution with a gain of 10 %. Also, Fig. 14c and14d shows the correlation between the convergence and the value of the angles. At first glance, there is no major difference between the two Figures. However, if we look a little closer, we see that both angle values in Fig. 14d are higher than those in Fig. 14c. Practically speaking, the lower the angle values the larger the plateau at the beginning of the cycle.

It reflects exactly the size of the plateau of red and black curves in the Fig. [START_REF] Jadoui | Comparative Study of Inner-Outer Krylov Solvers for Linear Systems in Structured and High-Order Unstructured CFD Problems[END_REF]. Apart from that, no difference appears.

VI. Conclusion

In this exploratory work, recycling Krylov subspace methods have been investigated for improving the partitioned solution of aerostructural adjoint systems.

Numerical experiments have been conducted on turbulent transonic flows over the three-dimensional ONERA M6 wing. In order to challenge our advanced Krylov solvers, a very flexible structural model of the wing was designed leading to a strongly coupled aeroelastic problem. Many of the improvements investigated in this study have proven beneficial to the convergence of Krylov solvers but often for problems of limited size or complexity.

At the beginning of this study, the current strategy was based on a FGMRES-DR solver for the fluid block with a varying right-hand side, i.e., the contribution of the structural source terms is updated at each fluid-structure coupling. However, this class of Krylov solvers is not able to take advantage of valuable spectral information produced during the previous fluid-structure cycles. This is mainly due to the fact that the deflation strategy is based on subspace augmentation. We then decided to implement a projection-based deflated Krylov solver which is able to virtually recycle any subspace at hand. The The current work particularly focused on the FGCRO-DR solver. More specifically, this class of solvers can recycle a spectral information from a linear system to another one.

A first step consisted in developing relevant numerical ingredients in order to improve robustness and efficiency of our inner-outer Krylov solvers. The combination of varying preconditioning strategies associated to proper diagonal scaling of the system matrix led to substantial improvement in the convergence rate, especially when the BILU(0) approximate factorization is applied to the exact first-order flux Jacobian matrix.

These numerical ingredients have been introduced in FGMRES-DR and FGCRO-DR and a numerical comparison showed the equivalence of these two solvers in terms of convergence rate and CPU cost. One important aspect which deserves to be mentioned is that scaling of the system matrix dramatically reduces the plateau of the residual convergence curve, which is a desirable feature in the framework of partitioned solvers when we want to achieve fast convergence between two fluid-structure couplings.

In the context of FGCRO-DR, we then explored subspace recycling from previous fluid-structure cycles. This showed very satisfactory convergence improvements, providing that the recycling process is activated after sufficient convergence of the fixed-point solver. Indeed, effective recycling of spectral information for the solution of a sequence of linear systems assumes that system matrix and right-hand side vary slowly. We proposed some distance criterion to monitor the beginning of recycling. We demonstrated a beneficial impact of recycling on the coupled adjoint system with gain of 20 % in terms of number of iterations.

Finally, we also investigated the recycling of error approximations for the aerodynamic adjoint system and for the standard (non flexible) GCRO-DR solver. Sequential and skip angles have been computed as a relevant criteria for appreciating the convergence. This loose formulation also brought significant improvements to the convergence rate of the Krylov solver. The next step will be to extend the error approximation recycling theory to the flexible class of GCRO solver and to demonstrate its added value in the contex of our partitioned coupled-adjoint solver.

  stru Approximate solution of the fluid adjoint problem 9: end if 10: for k ← 1, n cpl do 11:

Fig. 1 M6Fig. 2

 12 Fig. 1 M6 wing aeroelastic model: 42 block-structured RANS CFD mesh and FEM internal layout.

Fig. 3 Fig. 4

 34 Fig. 3 Comparison of rigid and aeroelastic pressure coefficient section plots at y=0.60m and y=1.12m.

a) Exact 2 Fig. 5

 25 Fig. 5 Stencils of the various Jacobian matrices.

Fig. 6

 6 Fig. 6 Comparison of the numerical performances between FGCRO-DR and FGMRES-DR. Impact of the orthogonalization step on the FGCRO-DR solver: one step orthogonalization (nMGS=1) deteriorates the convergence after 10 -8 while a re-orthogonalization restores convergence (nMGS=2). We recall the numerical parameters as follows: m = 30, m i = 10, k = 10.

Fig. 7

 7 Fig. 7 Adjoint relative residual norm convergence history of FGMRES-DR and FGCRO-DR. Impact of the scaling on the convergence. We recall the numerical parameters of solvers as follows: m = 30, m i = 10, = 10.

Fig. 8 Fig. 9

 89 Fig. 8 Adjoint relative residual norm convergence history of FGMRES-DR. Impact of various preconditioners.

Fig. 10

 10 Fig. 10 Variation of structural and aerodynamic right-hand sides with respect to F/S cycles. Residual norms are scaled by the norm of the fluid right-hand side.

Fig. 11 Fig. 12

 1112 Fig. 11 Relative residual convergence history for both FGCOR-DR FGMRES-DR. Impact of recycling on the convergence.

Fig. 13

 13 Fig. 13 Adjoint relative residual norm convergence history for GCRO-DR and LGCRO-DR. Impact of recycling the error approximations on the convergence. No scaling strategy is considered in this case.

13

 13 Fig.13shows how the recycling of the approximation improves the robustness of the GCRO-DR solver. In order to highlight the interest of the error approximation strategy, we decided to solve the aerodynamic adjoint system for two different sizes of the Krylov space. We consider both GCRO-DR(60,20) and GCRO-DR(30,10). As we can see in Fig.13, GCRO-DR(60,20) converges after 1900 iterations. Even though a plateau is observed at the beginning, the convergence of GCRO-DR(60,20) seems to follow the superlinear convergence property of GMRES without restarting. Fig.14adepicts the variation of the sequential and skip angles for GCRO-DR(60,20) throughout the cycles. At the beginning, we observe for both sequential and skip angles low values which corroborates the plateau that we see in Fig.13.

  Fig.13shows how the recycling of the approximation improves the robustness of the GCRO-DR solver. In order to highlight the interest of the error approximation strategy, we decided to solve the aerodynamic adjoint system for two different sizes of the Krylov space. We consider both GCRO-DR(60,20) and GCRO-DR(30,10). As we can see in Fig.13, GCRO-DR(60,20) converges after 1900 iterations. Even though a plateau is observed at the beginning, the convergence of GCRO-DR(60,20) seems to follow the superlinear convergence property of GMRES without restarting. Fig.14adepicts the variation of the sequential and skip angles for GCRO-DR(60,20) throughout the cycles. At the beginning, we observe for both sequential and skip angles low values which corroborates the plateau that we see in Fig.13.

  Fig. 14 Sequential and skip angles: (a) GCRO-DR(60,20), (b) GCRO-DR(30,10), (c) LGCRO-DR(30,9,1) and (d) LGCRO-DR(30,8,2).

Appendix A Flexible GCRO with deflated restarting (FGCRO-DR) a sequence of linear systems with varying right-hand sides.

Algorithm 2 FGCRO-DR(m,k)

1: Choose m, the maximum size of the subspace, and k, the desired number of approximate eigenvectors. let tol be the convergence tolerance. Choose an initial guess x 0 . Compute r 0 = b -Ax 0 , and set i = 1.

2: if C k and Z k are defined (from solving a previous linear system) then 3:

Perform m steps of FGMRES, solving min c -Hm y 2 for y and generating V m+1 , Z m and Hm . 9:

x 1 = x 0 + Z m y 10:

Compute the k eigenvectors gj of (H m + h 2 m+1,m H -T m e m e T m ) gj = θ j gj associated with the smallest magnitude eigenvalues θ j and store in P k .

13:

Let [Q,R] be the reduced QR-factorization of Hm P k .

14:

Solve min V T m+1 r i-1 -Ḡm y 2 for y.

26:

Compute the k eigenvectors gj of ḠT m Ḡm gj = θ j ḠT m V T m+1 W m gj associated with smallest magnitude eigenvalues θ j and store in P k .

29:

Let [Q,R] be the reduced QR factorization of Hm P k . 30:

33: end while 34: Recycling of Z k and C k (for the next system)