Laura Grigori
email: laura.grigori@inria.fr

Olivier Tissot
email: olivier.tissot@inria.fr.

Scalable Linear Solvers based on Enlarged Krylov subspaces with Dynamic Reduction of Search Directions

Keywords: Krylov subspace methods, Conjugate Gradient, Communication reducing algorithms

Krylov methods are widely used for solving large sparse linear systems of equations. On distributed architectures, their performance is limited by the communication needed at each iteration of the algorithm. In this paper, we study the use of so-called enlarged Krylov subspaces for reducing the number of iterations, and therefore the overall communication, of Krylov methods. In particular, we consider a reformulation of the Conjugate Gradient method using these enlarged Krylov subspaces: the enlarged Conjugate Gradient method. We present the parallel design of two variants of the enlarged Conjugate Gradient method as well as their corresponding dynamic versions where the number of search directions is dynamically reduced during the iterations. For a linear elasticity problem with heterogeneous coefficients, using a block Jacobi preconditioner, we show that this implementation scales up to 16, 384 cores, and is up to 5.7 times faster than the PETSc implementation of PCG.

Solveurs linéaires scalables basés sur des sous-espaces de Krylov Élargis avec réduction dynamique des directions de recherche

Résumé : Les méthodes de Krylov sont largement utilisées pour résoudre les systèmes linéaires creux et de grande taille. Sur une architecture distribuée, leur performance est limitée par les communications nécessaires à chaque itération de l'algorithme. Dans ce papier, on étudie l'usage de sous-espaces de Krylov élargis pour réduire le nombre d'itérations, et ainsi le total des communications, des méthodes de Krylov. En particulier, on considère une reformulation de la méthode du Gradient Conjugué qui utilise ces sous-espaces de Krylov élargis : le Gradient Conjugué élargi. On présente le design parallèle de deux variantes de la méthode du Gradient Conjugué élargi ainsi que les versions dynamiques associées, où le nombre de directions de recherche est réduit dynamiquement pendant les itérations. Pour un problème d'élasticité linéaire avec des coefficients hétérogènes, en utilisant un preconditioneur de type Jacobi par bloc, on montre que cette implémentation passe à l'echelle jusqu'à 16, 384 coeurs, et est jusqu'à 5, 7 fois plus rapide que l'implémentation de PCG présente dans PETSc.

Introduction

The discretization of partial differential equations, used to model physical phenomena, or optimization problems lead to linear systems of the form Ax = b where A is a sparse matrix. When A becomes very large, iterative methods based on Krylov subspaces are the method of choice [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. In this paper, we consider the case where A ∈ R n×n is symmetric (A = A) positive definite (x Ax > 0 for all x = 0). The Conjugate Gradient method [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF], and its preconditioned form, is a well-known method for solving such linear systems.

However, solving these linear systems efficiently on large scale computers remains a challenging problem. One difficulty is the high cost of communication compared to the computation on these machines [START_REF] Demmel | Minimizing communication in sparse matrix solvers[END_REF][START_REF] Demmel | Communicationoptimal parallel and sequential qr and lu factorizations[END_REF]. Recently, a lot of effort has been put in enhancing the performance of Krylov methods by avoiding global communication [START_REF] Chronopoulos | s-step iterative methods for (non)symmetric (in)definite linear systems[END_REF][START_REF] Carson | Avoiding communication in nonsymmetric lanczos-based krylov subspace methods[END_REF], overlapping communication with computation [START_REF] Ghysels | Hiding global synchronization latency in the preconditioned conjugate gradient algorithm[END_REF], or decreasing the number of iterations by searching in multiple directions at once [START_REF] Spillane | Adaptive Multi Preconditioned Conjugate Gradient: Algorithm, Theory and an Application to Domain Decomposition[END_REF][START_REF] Grigori | Enlarged Krylov Subspace Conjugate Gradient Methods for Reducing Communication[END_REF]. In this paper, we focus on the third approach, more precisely on the enlarged Conjugate Gradient method (ECG) [START_REF] Grigori | Enlarged Krylov Subspace Conjugate Gradient Methods for Reducing Communication[END_REF][START_REF] Grigori | Reducing the communication and computational costs of enlarged krylov subspaces conjugate gradient[END_REF].

After recalling Orthodir and Orthomin variants of ECG, we show the explicit link between the 2 methods. This gives a rigorous justification of an observation already made concerning the robustness of Orthodir compared to Orthomin in [START_REF] Grigori | Reducing the communication and computational costs of enlarged krylov subspaces conjugate gradient[END_REF]. Then we study theoretically the convergence behavior of ECG. We greatly improve the previous result in [START_REF] Grigori | Enlarged Krylov Subspace Conjugate Gradient Methods for Reducing Communication[END_REF], and show that ECG acts as if the smallest eigenvalues were somehow deflated. Then we present the parallel design of ECG. We consider both Orthodir and Orthomin variants, as well as dynamic versions of these variants that reduce dynamically the number of search directions in order to reduce the extra arithmetic cost in ECG compared to classical CG. In practice, we observe that enlarging the Krylov subspaces can drastically reduce the number of iterations. Indeed in the numerical experiments it is used with a block Jacobi preconditioner and acts as a second-level that, in a way, deflates the smallest eigenvalues; this is in accordance with the theory. This leads to a significant speed-up over the classical PCG. For instance for a 3D linear elasticity problem with heterogeneous coefficients with 4.5 millions of unknowns and 165 millions of nonzero entries, we observe that ECG is up to 5.7 times faster than the PETSc implementation of PCG, both using a block Jacobi preconditioner. This test case is known to be difficult because the classical onelevel preconditioners are not expected to be very effective [START_REF] Dolean | An introduction to domain decomposition methods[END_REF]. As it increases the arithmetic intensity and reduces the communication, it is well suited for modern and future architectures that exhibit massive parallelism. For the previous elasticity problem, we show that the method can scale up to 16, 384 threads, each one being bound to one physical core, which means that each core owns nearly 280 unknowns.

In summary, the contributions of the paper are the following. We provide a rigorous justification of the lack of robustness of Orthomin compared to Orthodir observed experimentally in [START_REF] Grigori | Reducing the communication and computational costs of enlarged krylov subspaces conjugate gradient[END_REF]. We give a proof of the speed of convergence of ECG -based on a direct extension of the proof of [6, Theorem 3.2] -which greatly improves the previous existing result presented in [START_REF] Grigori | Enlarged Krylov Subspace Conjugate Gradient Methods for Reducing Communication[END_REF]. This shows that ECG acts as a second level preconditioner that mitigates the effect of the smallest eigenvalues on the convergence of the iterative method. Hence it is sufficient to use as preconditioner a highly parallel method as Block Jacobi which bounds the largest eigenvalue of the preconditioned matrix. Finally, we introduce a parallel design embedding several variations of ECG whose scalability is assessed on different matrices and up to 16, 384 cores. We want to point out that our aim is not to design a specific solver for elliptic partial differential equations such as GenEO [START_REF] Spillane | Abstract robust coarse spaces for systems of pdes via generalized eigenproblems in the overlaps[END_REF] or multigrid preconditioners with some tuning. For a detailed comparison of such solvers we refer to Jolivet's thesis [START_REF] Jolivet | Domain decomposition methods. Application to highperformance computing[END_REF]. It is very likely that for these test cases, these solvers are more effective than ECG with a block Jacobi preconditioner. Nevertheless, unlike these methods ECG is an algebraic method. It does not require any information from the underlying PDE and does not rely on any assumption, except that the matrix is symmetric positive definite. Hence it can be seen as a black-box solver and integrated very easily in any existing code.

Enlarged Krylov Conjugate Gradients

Block Krylov methods

In 1980 O'Leary introduced the block Conjugate Gradient method [START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF] for solving SPD systems with several right-hand sides. In this seminal paper, she proved that block CG can converge significantly faster than CG. This idea was then generalized and extended to other classical Krylov methods as GMRES [START_REF] Robbé | Exact and inexact breakdowns in the block GMRES method[END_REF][START_REF] Langou | Iterative methods for solving linear systems with multiple righthand sides[END_REF] or BiCGSTAB [START_REF] El Guennouni | A block version of bicgstab for linear systems with multiple right-hand sides[END_REF]. Later Gutknecht [START_REF] Gutknecht | Block Krylov space methods for linear systems with multiple right-hand sides: an introduction[END_REF] introduced a general framework for defining Block Krylov subspaces.

Recently, Block Krylov methods are receiving an increasing attention in the HPC field [START_REF] Calandra | A modified block flexible GMRES method with deflation at each iteration for the solution of non-hermitian linear systems with multiple right-hand sides[END_REF][START_REF] Agullo | Block gmres method with inexact breakdowns and deflated restarting[END_REF][START_REF] Jolivet | Block iterative methods and recycling for improved scalability of linear solvers[END_REF][START_REF] Spillane | An adaptive multipreconditioned conjugate gradient algorithm[END_REF][START_REF] Kreutzer | GHOST: Building blocks for high performance sparse linear algebra on heterogeneous systems[END_REF]. They appear to be well suited for modern computers' architectures with a high level of parallelism because they allow to reduce the number of global synchronizations, while also featuring a higher arithmetic intensity at the cost of some extra computations.

Enlarged Krylov subspaces

In [START_REF] Grigori | Enlarged Krylov Subspace Conjugate Gradient Methods for Reducing Communication[END_REF], the authors define so-called enlarged Krylov subspaces. First, the matrix A is reordered by partitioning its graph into N subdomains (using METIS [START_REF] Karypis | Metis -unstructured graph partitioning and sparse matrix ordering system[END_REF] for example). Then, the initial residual r 0 is split into t vectors denoted R e 0 (i) , 1 ≤ i ≤ t. In the original paper the authors use t = N . It is important to note that the case t < N can be dealt with many ways as long as r 0 = t i=1 R e 0 (i) (Fig. 1). This is of particular interest in practice because typically N will correspond to the number of MPI processes. The parameter t is called the enlarging factor. In practice for a given t the splitting of r 0 does not have a high impact on the convergence of the method. In the numerical experiments we construct the initial enlarged residual R e 0 = [R e 0 (1) , . . . , R e 0 (t) ] as the leftmost example in Fig. 1.

RR n°9190

Figure 1: Illustration of the ordering of A into 8 subdomains obtained with METIS [START_REF] Karypis | Metis -unstructured graph partitioning and sparse matrix ordering system[END_REF] and several admissible splittings of r 0 into 3 vectors.

Then, the enlarged Krylov subspace of order k denoted K k,t (A, r 0 ) is defined as the block Krylov subspace of order k associated to A and the enlarged residual R e 0 . More precisely, and following the notation introduced in [START_REF] Gutknecht | Block Krylov space methods for linear systems with multiple right-hand sides: an introduction[END_REF],

K k,t (A, r 0 ) : = K k (A, R e 0 ) (1)
= span R e 0 , AR e 0 , . . . , A k-1 R e 0 .

(

) 2
Using this definition and following [START_REF] Grigori | Reducing the communication and computational costs of enlarged krylov subspaces conjugate gradient[END_REF] it is possible to derive two variants (Orthomin and Orthodir) of the enlarged Conjugate Gradient (ECG) algorithm (Fig. 1). More precisely, the enlarged approximate solution is a matrix of size n × t denoted X k , and the sum of its columns gives the approximate solution of the original system. We denote R k the enlarged approximate residual, and similarly we obtain the approximate residual of the original system by summing its columns. P k is a matrix of size n×t called search directions, it corresponds to the A-orthonormalization of Z k . We denote α k the optimal step, unlike in CG algorithm it is not a scalar but a matrix of size t × t. Depending on the method for constructing Z k+1 , it is possible to derive two variants of ECG: Orthomin and Orthodir.

Orthomin (Omin) corresponds to Block CG [START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF]:

β k = (AP k ) R k , (3)
Z k+1 = R k -P k β k . ( 4
)
This method is very similar to the one originally proposed by Hestenes and Stiefel [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF] because it constructs the new descent directions Z k+1 using R k and P k . Orthodir (Odir) corresponds to the Block Lanczos algorithm but with the inner product induced by A:

γ k = (AP k ) (AP k ), (5)
ρ k = (AP k-1 ) (AP k ), (6)
Z k+1 = AP k -P k γ k -P k-1 ρ k . (7)
It is the block equivalent of the homonym method defined in [START_REF] Ashby | A taxonomy for conjugate gradient methods[END_REF]. Unlike the previous variant, Z k+1 is constructed using P k and P k-1 .

Both Orthodir and Orthomin produce Z k+1 that is A-orthogonal to P i for i ≤ k. Then the search directions P k+1 are defined as

P k+1 = Z k+1 (Z k+1 AZ k+1 ) -1/2 . ( 8
)
Unlike CG algorithm a breakdown would occur if Z k+1 AZ k+1 is singular, i.e., Z k+1 is not full rank. Although rare this situation can happen in practice and several variants have been developed in order to handle this case [START_REF] Ji | A breakdown-free block conjugate gradient method[END_REF][START_REF] Grigori | Reducing the communication and computational costs of enlarged krylov subspaces conjugate gradient[END_REF][START_REF] Dubrulle | Retooling the method of block conjugate gradients[END_REF][START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF]. Overall, both Orthomin and Orthodir generate P k+1 such that

P k+1 AP i = 0, ∀i ≤ k (9)
P k+1 AP k+1 = I. (10)
Consequently, the ECG method can be summarized in Algorithm 1. Another difference with the original block CG algorithm is that the search directions are A-orthonormalized at each iteration: P k is used as search directions instead of Z k . It has be shown numerically that using this variant can increase the numerical stability of the method [START_REF] Dubrulle | Retooling the method of block conjugate gradients[END_REF].

Given a preconditioner M -1 , the idea for applying left preconditioning to the (block) Conjugate Gradient method is to remark that M -1 A is self-adjoint with respect to the M -inner product [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. Then by replacing A by M -1 A, and the transpose by M in the algorithm (Fig. 1), it follows the preconditioned enlarged Conjugate Gradient method. In fact, some simplifications occur and the algorithm remains exactly the same except the definition of Z k that is slightly different. More precisely, it follows that the preconditioned Orthomin method corresponds to, Z k+1 = (I -

P k P k A)M -1 R k , (11)
and the preconditioned Orthodir method corresponds to,

Z k+1 = (I -P k P k A -P k-1 P k-1 A)M -1 AP k . (12)
In both cases, the initialization also slightly differs because

Z 1 = M -1 R e 0 .
Overall, the preconditioner is applied once per iteration, as in the classical CG method.

Equivalence between Orthodir and Orthomin

In what follows, we assume exact arithmetic and we study the connection between these two methods with the aim to derive formulas that links the approximate quantities of both variants. Indeed, by construction the approximate

RR n°9190

Algorithm 1 ECG algorithm.

1:

P 0 = 0 2: Z 1 = R e
P k = Z k (Z k AZ k ) -1/2
6: construct Z k+1 using (3)-(4) (Orthomin) or ( 5)-( 7) (Orthodir)

α k = P k R k-1 7: X k = X k-1 + P k α k 8: R k = R k-1 -AP k α k 9: if || t i=1 R (i) k || 2 < ε then 10:
13:

k = k + 1 14: end for 15: x k = t i=1 X (i) k
solutions computed by Orthodir and Orthomin are equal. Hence, the approximate residuals are also equal. But this does not imply that the search directions generated are equal even if they belong to the same space. We denote with a tilde the variables related to Orthomin and with a hat the variables related to Orthodir, e.g., P k are the A-orthonormalized search directions generated during Orthodir. In order to simplify the presentation we consider that no breakdowns have occurred, i.e., Z k , P k and Z k , P k are all well-defined.

Since Orthomin and Orthodir rely on the same projection process [START_REF] Grigori | Reducing the communication and computational costs of enlarged krylov subspaces conjugate gradient[END_REF] (they both search an approximate solution in K t,k , such that the corresponding residual is orthogonal to K t,k ), we know that X k = X k . It follows that:

P k α k = P k α k , (13)
P k P k = P k P k . (14)
Hence, there exists δ k ∈ R t×t orthogonal and such that P k = P k δ k . By definition we also have,

R k -R k-1 = -A P k α k . (15)
A simple computation using the previous relationships gives,

-Z k+1 δ k α k = -A P k α k + P k γ k α k + P k-1 ρ k α k , (16)
= R k -P k β k -Z k + P k P k AR k-1 , (17)
and,

P k P k AR k-1 = Z k ( Z k A Z k ) -1 Z k AR k-1 , (18)
= Z k ( Z k A Z k ) -1 Z k A Z k , (19)
= Z k . (20)
Thus,

Z k+1 = -Z k+1 δ k α k . ( 21
)
This result is a generalization of a previous result presented by Manteuffel et al. [2, p. 1550] for the classical CG. In fact, the authors show that z k = Π k i=0 (-α k ) z k but they never consider explicitly the A-orthonormalized search directions. In particular, they define z k+1 using z k (for Omin) and z k-1 (for Odir). This explains the slight difference between our generalization and their result.

When k becomes large,

α k = P k R k-1 and || α k || 2
is more likely to be low because R k-1 is supposed to converge to 0 and P k is A-orthonormalized -the same reasoning applies for α k . This result is very interesting because it shows that, since δ k is an orthogonal matrix, when k becomes large || Z k+1 || 2 can be significantly lower than || Z k+1 || 2 . Hence, the conditioning of Z k+1 A Z k+1 could be much worse than that of Z k+1 A Z k+1 , possibly leading to a breakdown when computing its Cholesky factorization (line 5 in Fig. 1). It is remarkable to notice that even for classical CG method, this has already been noticed by Manteuffel et al. in [2, p. 1551-1552]: "If BCA is indefinite, Omin may still be used, but the previous direction vector [...] should be stored. Then, if α i = 0, control can switch to the 3-term recursion of Odir to get p i+1 ". In practice, this phenomenon is indeed observed: there are cases where Orthomin breaks down while Orthodir does not [START_REF] Grigori | Reducing the communication and computational costs of enlarged krylov subspaces conjugate gradient[END_REF]. In conclusion, Orthodir is expected to be more reliable than Orthomin. However, Orthodir is also more costly than Orthomin: the construction of Z k+1 requires twice as many flops and memory as for Orthomin.

Convergence study

As previously mentioned, O'Leary [START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF] proved that block CG can converge significantly faster than the classical CG. In [START_REF] Grigori | Enlarged Krylov Subspace Conjugate Gradient Methods for Reducing Communication[END_REF], it is proved that ECG converges faster than CG but there is no further information on the speed of convergence of ECG. This section is dedicated to the proof of the following result.

Theorem 1 Let x k be the approximate solution given by the Enlarged Conjugate Gradient method with an enlarging factor t at step k. Then we have:

||x k -x * || 2 A ≤ C √ κ t -1 √ κ t + 1 2k ( 22
)
where κ t = λ n λ t and C is a constant independent of k whose exact expression will be given in the course of the proof.

Before starting the proof, we want to point out that this result is a big improvement of the theorem stated in [START_REF] Grigori | Enlarged Krylov Subspace Conjugate Gradient Methods for Reducing Communication[END_REF] because it explains that ECG's convergence is, to some extent, closer to that of Deflated-CG [START_REF] Dostál | Conjugate gradient method with preconditioning by projector[END_REF] rather than that of the classical CG.

Proof 1

The key idea of the proof is to remark the close link between ECG and block CG. In particular, we follow the proof of Theorem 3.2 by Jie Chen in [START_REF] Chen | A deflated version of the block conjugate gradient algorithm with an application to gaussian process maximum likelihood estimation[END_REF] for block CG, and we adapt it to our case of interest: ECG. The approach of the proof is similar to that of O'Leary [START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF] but the final result is slightly different.

First, we need to write the error at iteration k as a polynomial of A evaluated in the initial error. Indeed,

x k = x 0 + t j=1 q kj (A)Ad (j) 0 , (23)
where d

(j) 0

is the i-th component of the enlarged initial error [START_REF] Chen | A deflated version of the block conjugate gradient algorithm with an application to gaussian process maximum likelihood estimation[END_REF], and q kj is a polynomial of degree not exceeding k -1.

Hence,

e k = x k -x * = t j=1 (1 + q kj (A)A)d (j) 0 (24) = t j=1 p kj (A)d (j) 0 , (25)
where p kj (X) = 1 + q kj (X)X is a polynomial of degree not exceeding k and such that p kj (0) = 1, ∀k, j.

Let A = Φ ΛΦ be the spectral decomposition of A. Let us rewrite e k according to this decomposition,

e k = t j=1 p kj (A)d (j) 0 = t j=1 Φp kj (Λ)Φ d (j) 0 . (26)
Thus,

||e k || 2 A = e k Ae k (27) =   t j=1 Φp kj (Λ)Φ d (j) 0   A   t j=1 Φp kj (Λ)Φ d (j) 0   (28) =   t j=1 d (j) 0 Φp kj (Λ)   A   t j=1 p kj (Λ)Φ d (j) 0   . ( 29
)
This final expression is a generalization of the expression that occurs in the proof of convergence of the CG algorithm in [START_REF] Saad | Iterative methods for sparse linear systems[END_REF].

Let j ∈ {1, . . . , t}. We denote p := p k ĵ in order to simplify the notations. Following Chen [START_REF] Chen | A deflated version of the block conjugate gradient algorithm with an application to gaussian process maximum likelihood estimation[END_REF], let us define

-p(Λ) -1 (2I -p(Λ))Φ d (j) 0 = F 1 F 2 , ∀j = j, (30)
Φ d

( j) 0 = f 1 f 2 . ( 31
)

RR n°9190

The other polynomials p kj are chosen such that p kj = τ j (2 -p) where the τ j are defined as the components of the solution of

F 1 τ = f 1 .
From those definitions, it follows that for j = j p kj (Λ)Φ d

(j) 0 = (2I -p(Λ))Φ d (j) 0 τ j (32) = -p(Λ)(-p(Λ) -1 (2I -p(Λ))Φ d (j) 0 )τ j (33) = -p(Λ) F 1 F 2 τ j . (34)
Hence, when summing up those terms

j = j p kj (Λ)Φ d (j) 0 = -p(Λ). (35)
Thus,

t j=1 p kj (Λ)Φ d (j) 0 = p(Λ) f 1 f 2 -p(Λ) f 1 F 2 F -1 1 f 1 (36) = 0 0 -p(Λ 2 )F 2 F -1 1 p(Λ 2 ) Φ d ( j) 0 (37)
:= 0 0 E P Φ d ( j) 0 . (38)
Finally,

||e k || 2 A = d ( j) 0 Φ 0 E 0 P A 0 0 E P Φ d ( j) 0 (39)
≤ ||d

( j) 0 || 2 A E E E P P E P 2 . ( 40
)
Let us focus on the right term in this last inequality,

E E E P P E P 2 ≤ ||P || 2 0 0 F 2 F -1 1 I 2 (41) ≤ ||P || 2 a 2 + 1 ( 42
)
where a is the largest singular value of F 2 F -1 1 . Finally, we have

||e k || 2 A ≤ ||d ( j) 0 || 2 A a 2 + 1 ||P || 2 . ( 43
)
If we replace p (and therefore P ) by the optimal choice we can rewrite the bound as the following min-max problem,

||e k || A ≤ ||d ( j) 0 || A (a 2 + 1) min p∈P k 1 max t≤i≤n |p(λ i )|. ( 44
) RR n°9190
Finally, it is possible to use Chebyshev polynomials to estimate this min-max quantity [START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF][START_REF] Saad | Iterative methods for sparse linear systems[END_REF],

min p∈P k 1 max t≤i≤n |p(λ i )| ≤ 2 √ κ t -1 √ κ t + 1 k . ( 45
)

Dynamic reduction of the search directions

In what follows, we recall an approach for reducing the block size in the Orthodir method during the iterations presented in [START_REF] Grigori | Reducing the communication and computational costs of enlarged krylov subspaces conjugate gradient[END_REF]. The idea is to reduce the arithmetic and memory extra costs of Orthodir while maintaining its good convergence behavior. As explained in the survey [START_REF] Gutknecht | Block Krylov space methods for linear systems with multiple right-hand sides: an introduction[END_REF] the key idea to reduce the block size is to monitor the rank of R k . Once R k becomes rank deficient, it means that a part of the approximate solution has already converged at iteration k. More precisely for i ≥ k -1, there exists a linear combination (independent of i) of columns of X i that remains constant. As a consequence, there exists a linear combination of search directions that is not useful anymore to compute the approximate solution. The idea is to remove these search directions in the next iterations. As R k-1 is an n × t matrix with n large, it is preferable to avoid computing the rank of R k-1 directly. In [START_REF] Grigori | Reducing the communication and computational costs of enlarged krylov subspaces conjugate gradient[END_REF], it is shown that the rank of

α k = P k R k-1 can be computed instead.
The method presented in [START_REF] Grigori | Reducing the communication and computational costs of enlarged krylov subspaces conjugate gradient[END_REF] can be divided into two parts. At each iteration of the algorithm (Fig. 2) a Singular Value Decomposition of α k is computed (line 8). If the numerical rank of α k is below a given tolerance then the search directions are reduced accordingly (line 12) and some of them are kept in order to keep the A-orthogonality property (line 14, 24 and 25). Although computing the SVD of α k at each iterations induces an extra cost compared to Orthodir, this operation does not involve any communications and it is negligible because α k is a small matrix of size t × t. Furthermore, as the search directions P k are reduced, the dominant operation of Krylov iterations in terms of flops; the matrix product (AP k ), and the application of the preconditioner (M -1 AP k ), are cheaper.

Curing breakdowns in Orthomin

As explained previously Orthomin version of ECG can break down. There exist several methods to overcome this issue and in the following we recall the Breakdown-Free block CG method defined in [START_REF] Ji | A breakdown-free block conjugate gradient method[END_REF]. Starting from the original algorithm of O'Leary [START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF], the authors propose to perform a rank-revealing QR decomposition of Z k+1 and then drop its null part before A-orthonormalizing it. They show that in exact arithmetic this allows to continue the algorithm with nearly no further modification. The resulting algorithm is given in Algorithm 3.

From a practical point of view the size of P k+1 can be reduced, but at each iteration Z k+1 is of size n × t because the size of R k remains constant. Hence the matrix product AP k is cheaper but the application of the preconditioner

RR n°9190

Algorithm 2 ECG D-Odir algorithm.

1:

P 0 = 0 2: Z 1 = R e 0 3: H = ∅ 4: k = 1 5: for k = 1, . . . ,
P k = Z k (Z k AZ k ) -1/2
7:

α k = P k R k-1 8: α k = U k Σ k V k 9:
let s k be the number of singular values of α k bigger than ε def 10:

if s k < s k-1 then 11: α k = U k α k 12: P k = P k U k 13:
α k = α k (1 : s k , :)

14:

H = [H, P (:, s k : s k-1 )]

15:

P k = P k (:, 1 : s k ) 16:
end if 17:

X k = X k-1 + P k α k 18: R k = R k-1 -AP k α k 19: if || t i=1 R (i) k || 2 < ε then 20:
γ k = (AP k ) (AP k ) 23: ρ k = (AP k-1 ) (AP k ) 24: δ k = (AH) (AP k ) 25: Z k+1 = AP k -P k γ k -P k-1 ρ k -Hδ k 26: k = k + 1 27: end for 28: x k = t i=1 X (i) k RR n°9190
M -1 R k is not. Furthermore computing a rank-revealing QR factorization of Z k+1 cannot be neglected because Z k+1 is of size n × t. In summary, as this method was not meant for efficiency, but rather for improving the stability of Orthomin, it does not allow to save as many computations as in the dynamic variant of Orthodir.

Algorithm 3 BF-ECG algorithm.

1:

P 0 = 0 2: Z 1 = R e
P k = Z k (Z k AZ k ) -1/2
6:

α k = P k R k-1 7: X k = X k-1 + P k α k 8: R k = R k-1 -AP k α k 9: if || t i=1 R (i) k || 2 < ε then 10:
β k = (AP k ) R k 13: Z k+1 = R k -P k β k 14:
Z k+1 = ortho(Z k+1 )

15:

k = k + 1 16: end for 17: x k = t i=1 X (i) k 3 Parallel design 3.

Data distribution

As it is usually the case in parallel implementations of Krylov methods, we assume that the unknowns are distributed among the processors. We also assume that each processor owns different unknowns. Thus all the variables whose size scales as the size of the linear system (X k , R k , P k , AP k , Z k ) are distributed row wise among the processors according to the distribution of the unknowns. All variables whose size scales as the enlarging factor t (α k , β k , γ k , ρ k ) are replicated on all the processors. Locally they are stored contiguously and column by column. There is no allocation or deallocation of memory during the iterations. In particular, when using Dynamic Orthodir or Breakdown-Free Orthomin the memory is not freed when the block size is reduced. The local memory consumption of preconditioned Orthodir and Orthomin on P processors is summarized in Table 1. For completeness, we also add the local memory consumption of the classical CG algorithm, described in [START_REF] Saad | Iterative methods for sparse linear systems[END_REF] for instance, where only 5 vectors and 2 scalars are needed.

[P k ; P k-1 ] [AP k ; AP k-1 ] Z k R k X k α k [γ k ; ρ k ]/β k

Cost analysis of ECG

Our implementation of ECG is based on Reverse Communication Interface [START_REF] Kalhan | Reverse communication interface for linear algebra templates for iterative methods[END_REF]. For one iteration of ECG, it requires external routines to apply the sparse matrix product and the preconditioner to a set of vectors. Indeed the implementation of these routines highly depends on the linear system to be solved. This is why we do not take into account these operations in our cost analysis.

Given n, t such that t n, we denote V, W tall and skinny matrices of size n × t whose rows are distributed among the processors, and α is a matrix of size t × t replicated on the P processors. Following [START_REF] Kreutzer | GHOST: Building blocks for high performance sparse linear algebra on heterogeneous systems[END_REF], it is possible to decompose the iterations of ECG (and more generally block CG) into the following kernels: [START_REF] Kreutzer | GHOST: Building blocks for high performance sparse linear algebra on heterogeneous systems[END_REF]),

• V ← V + W α (tsmm in
• α ← V W (tsmtsm in [START_REF] Kreutzer | GHOST: Building blocks for high performance sparse linear algebra on heterogeneous systems[END_REF]),

• Cholesky factorization of α (potrf),

• triangular solve of α with several right-hand sides (trsm).

Following ECG algorithm (fig. 1), each iteration of Orthodir and Orthomin consists of 3 tsmm (lines 7, 8, and 12), 4 tsmtsm (lines 5, 6, 9, and 12), 1 potrf (line 5) and 2 trsm (line 5). Indeed, the line 5 of the algorithm (fig. 1) can be decomposed as, Doing so allows us to avoid calling the sparse matrix set of vectors product for computing AP k at the price of an extra trsm. Hence the difference between the two algorithms is the construction of Z k+1 (line 12). The tsmtsm and tsmm for constructing Z k+1 in Orthodir (equations ( 5)-( 7)) are twice as much costly as for Orthomin (equations ( 3)-( 4)).

AZ k ← A * Z k C ← tsmtsm(Z k , AZ k ) C ← potrf(C) P k ← trsm(Z k , C) AP k ← trsm(AZ k , C)
As matrices of size t × t are replicated among the processors, we notice that tsmm, Cholesky factorization of α and triangular solve of α are local operations without any communication. Hence we use the corresponding LAPACK routines: gemm, potrf (dense Cholesky factorization) and trsm (dense triangular solve with several right-hand sides). However V and W are distributed and tsmtsm is not a local operation. The LAPACK routine gemm is called to compute the local product V i W i followed by a call to MPI_Allreduce.

Thus, the only kernel operation that requires a communication is tsmtsm and 4 calls to MPI_Allreduce are done per iteration. It is usually assumed that during a call to MPI_Allreduce the number of messages sent and received on the network is equal to log 2 (P ) -although the exact number depends on the MPI implementation [START_REF] Thakur | Optimization of collective communication operations in mpich[END_REF]. Moreover it is a blocking operation: when completed all the processors are synchronized. This is why in practice, as in the classical CG, the communication cost is dominated by 2 calls to MPI_Allreduce: the one after the sparse matrix set of vectors product (line 5) and the one after the preconditioner (line 12), because they occur after operations with a potential load imbalance between processors.

In summary, the detailed costs of one iteration of Orthodir and Orthomin in terms of flops, words, and messages are indicated in Table 1. For the sake of comparison, we recall the complexity of the CG algorithm described in [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. We also report the number of MPI_Allreduce in parenthesis, in addition to the order of magnitude of the number of messages. In summary, one iteration of ECG is approximately t 2 times more costly in terms of flops than one iteration of CG. While the number of messages is of the same order, the number of words is also t 2 times larger. Indeed there is a trade-off between these extra costs and the reduction of the number of iterations due to the enlargement of the search spaces.

Cost of dynamic reduction of ECG

The implementation of the dynamic reduction of the search directions within Orthodir follows Algorithm 2. In practice, we use LAPACK routine gesvd and only compute the right singular vectors of α k denoted U k . We check the singular values obtained. If there are some lower than ε √ t , which is the criterion proposed in [START_REF] Grigori | Reducing the communication and computational costs of enlarged krylov subspaces conjugate gradient[END_REF], we call geqrf on U in order to perform the updates P U k , AP U k and U k α k in-place with ormqr. Since P k and AP k are stored in a column major fashion the selection of the columns is done at no cost. Similarly H is not explicitly defined. However the selection of the first rows of α k implies an inplace memory rearrangement.

The implementation of Breakdown-Free Orthomin is similar to Orthomin except the computation of a rank-revealing QR decomposition of Z k+1 . As Z k+1 is distributed, it is not reasonable to use a LAPACK kernel to compute it. Instead we use a modification of Chol-QR algorithm [START_REF] Yamazaki | Mixed-precision cholesky qr factorization and its case studies on multicore cpu with multiple gpus[END_REF] which is a cheaper but less stable alternative to TS-RRQR [START_REF] Demmel | Communication avoiding rank revealing qr factorization with column pivoting[END_REF][START_REF] Demmel | Communicationoptimal parallel and sequential qr and lu factorizations[END_REF]. Its implementation is very easy using the LAPACK routine pstrf (Cholesky with pivoting) for computing (R, π) at line 2. Following [START_REF] Ji | A breakdown-free block conjugate gradient method[END_REF] we use the default tolerance of pstrf for detecting exact rank deficiency of Z k+1 .

Algorithm 4 Chol-RRQR Input: P , ε Output: Q 1 orthogonal such that

P π = Q 1 Q 2 R 11 R 12 0 R 22
where π is a permutation and all the diagonal elements of R 11 are larger than ε 1: µ ← P P

2: Compute (R, π) such that π µπ = R R with R = R 11 R 12 0 R 22
and all the diagonal elements of R 11 are larger than ε 2 3: P 1 ← P π(:, 1 : size(R 11 ))

4: Q 1 ← P 1 R -1 11 4 Numerical experiments 4.

Description of the parallel environment

In the experiments we use a block Jacobi preconditioner, associating at each block a MPI process. Before calling ECG, each MPI process factorizes the diagonal block of A corresponding to the local row panel that it owns. At each iteration of ECG, each MPI process performs a backward and forward solve locally in order to apply the preconditioner. Hence the application of the preconditioner does not need any communication. It is likely that there exists better preconditioners than block Jacobi for our test cases, however we are interested in the iterative method rather than in the preconditioner. In particular, we do not want to target specific applications and aim at being as generic as possible. Although in theory it is possible to apply any preconditioner within this implementation, in practice it is essential that applying this preconditioner to several vectors at the same time is not too costly, e.g, a sublinear complexity with respect to the number of vectors.

The following experiments are performed on a machine located at Umeå University as part of High Performance Computing Center North (HPC2N), called Kebnekaise. It is a heterogeneous machine formed by a mix of Intel Xeon E5-2690v4 (Broadwell) with 2x14 cores (and E7-8860v4 for large memory computations), Nvidia K80 GPU and Intel Xeon Phi 7250 (Knight's Landing) with 68 cores. In our experiments, we use the so-called compute nodes, which are formed by Intel Xeon E5-2690v4 (Broadwell) with 2x14 cores. For a detailed description of the machine, we refer to the online documentation 1 .

We compile the code (and its dependencies) using Intel toolchain installed on the machine: mpiicc (based on icc version 18.0.1 20171018) and MKL [START_REF] Wang | Intel math kernel library[END_REF] version 2018.1.163. We use PETSc [START_REF] Balay | PETSc users manual[END_REF] in order to compare ECG implementation to PETSc PCG implementation. In particular, PETSc is configured to use MKL-PARDISO as exact solver for sparse matrices in the block Jacobi preconditioner. For partitioning the matrix we are using the METIS library downloaded and installed by PETSc.

Test cases

The Ela matrices arise from the linear elasticity problem with Dirichlet and Neumann boundary conditions defined as follows div(σ(u)) + f = 0 on Ω (46)

u = 0 on ∂Ω D (47) σ(u) • n = 0 on ∂Ω N ( 48
)
where Ω is a unit cube. The matrices Ela N correspond to this equation discretized with FreeFem++ [START_REF] Hecht | New development in freefem++[END_REF] using a triangular mesh with 1600 × N × N points on the corresponding vertices and P1 finite elements scheme. ∂Ω D is the Dirichlet boundary, ∂Ω N is the Neumann boundary, f is some body force, u is the unknown displacement field. σ(.) is the Cauchy stress tensor given by Hooke's law: it can be expressed in terms of Young's Modulus E and Poisson's ratio ν. For a more detailed description of the problem see [START_REF] Jolivet | Scalable domain decomposition preconditioners for heterogeneous elliptic problems[END_REF][START_REF] Grigori | Robust algebraic Schur complement preconditioners based on low rank corrections[END_REF][START_REF] Spillane | An adaptive multipreconditioned conjugate gradient algorithm[END_REF]. We consider a heterogeneous beam made of several layers of a hard material (E 1 , ν 1 ) = (2 × 10 11 , 0.25) and a soft material (E 2 , ν 2 ) = (10 7 , 0.45), i.e., discontinuous E and ν (Figure 3). This test case is known to be difficult because the matrix is ill conditioned. In particular, the classical one-level preconditioners are not expected to be very effective [START_REF] Dolean | An introduction to domain decomposition methods[END_REF]. As previously pointed out, ECG is an algebraic method that does not rely on any particular assumption on the matrix, except that it is symmetric positive definite. As an illustration, we also test the implementation on two SPD matrices coming from the Sparse Matrix Collection of Tim Davis [START_REF] Davis | The university of florida sparse matrix collection[END_REF]: Flan 1564 and Bump 2911. Numerical properties of the test matrices are summarized in Table 2.

Results

In all the experiments the tolerance is set as the default tolerance of PETSc, i.e., 10 -5 and the maximum number of iterations is set to 5000. The right-hand side is chosen uniformly random and normalized and the initial guess is set to 0. We do not use any kind of threading and use 28 MPI processes per node. Unless otherwise stated, we use one OpenMP thread per MPI process -we also perform numerical experiments to observe the effect of threading in the last section.

Impact of the enlarging factor

First we study the impact of the enlarging factor t on the methods. We fix the number of processors and we vary the value of t for the 4 methods: Orthodir (Odir), Orthodir with dynamic reduction of the search directions (D-Odir), Orthomin (Omin) and Breakdown-Free Orthomin (BF-Omin). The results obtained are summarized in Table 3 For Flan 1565 the number of MPI processes is fixed to 56. We remark that the runtime is decreasing until t = 12 and then it increases slightly. When t is relatively small the 4 methods are comparable. However as t increases the effect of dynamic reduction becomes more visible. With t = 28, D-Odir is almost 10% faster than Odir. On the other hand, as we are detecting exact rank deficiency of Z k+1 , BF-Omin did not reduce the size of the block and as no breakdown occurs it is slightly more costly than Omin. We also tested ε √ t as the tolerance for detecting breakdowns but this does not allow the method to converge. Overall, for this matrix, the best method is D-Odir with t = 12.

For Bump 2911 we fix the number of MPI processes to 112. For this matrix the reduction of the number of iterations is not balancing the increase in flops. For instance, the number of iterations for D-Odir(8) is 695, and for D-Odir [START_REF] Dolean | An introduction to domain decomposition methods[END_REF] it is 665, which represents a decrease of only 4% in the number of iterations. According to Theorem 1, it is very likely that in this case the preconditioned matrix does not have a cluster of small eigenvalues, hence the convergence of the method is not significantly improved when enlarging the Krylov subspace. However, we also notice that using the dynamic Orthodir variant (D-Odir) allows to reduce significantly the runtime when t is large: D-Odir is around 20% faster than Odir.

We also perform this study for Ela 20 with 112 MPI processes. First we remark that a breakdown occurs with Orthomin for all the values of t that we tested. This behavior of elasticity matrices had also been reported in [START_REF] Grigori | Reducing the communication and computational costs of enlarged krylov subspaces conjugate gradient[END_REF]. Using BF-Omin effectively cures the breakdowns but does not allow the method to converge within the prescribed maximum number of iterations. Similarly, D-Odir does not converge when t = 4, but performs very well when t is larger. On the contrary, Odir is very stable and converges for all the values of t tested. As for Bump 2911 we observe that D-Odir is around 20% faster than Odir when t = 24. Overall, for this matrix, the best method is D-Odir with t = 24.

In order to better understand how the dynamic reduction of the search directions affects the convergence we plot the normalized residual and the block size (dash line) as a function of the iteration count for Flan 1565 (Fig. 4a) and Ela 20 (Fig. 4b). We notice that the convergence is not really affected by the reduction of the search directions because the number of iterations remains almost the same. However the block size is effectively reduced as soon as the method starts to converge. We note that the Ela 20 test case is very favorable: the block size is reduced even a bit before the convergence and the number of iterations is lower than when the search reductions are not reduced.

In conclusion, D-Odir is the best method over the different variants of ECG that we tested: it is a good compromise between the stability of Odir and the efficiency of the classical CG. Nevertheless, there exists matrices such as Bump 2911 for which the reduction of the number of iterations does not compensate the extra cost of ECG compared to the classical CG, even when using the dynamic reduction of the search directions. These results support the theoretical convergence study that was done in the previous section. ECG(t) is acting as if the t smallest eigenvalues of the matrix were deflated. Finally, we notice that values of t between 8 to 24 are good default parameters. Indeed, such values allow to effectively reduce the number of iterations while maintaining an affordable cost per iteration.

Strong scaling study

Following the parameter study, we perform a strong scaling study on Flan 1565 and Ela 30. As Bump 2911 does not seem particularly well suited for the method we do not perform the strong scaling study on this matrix.

For Flan 1565 we compare PETSc PCG and D-Odir with t = 12, the best choice over the parameters we tested. The resulting runtimes are presented in Table 4. When the number of MPI processes is relatively low ECG scales as well as PETSc, i.e., almost linearly. As the number of iterations is significantly reduced with D-Odir [START_REF] Kalhan | Reverse communication interface for linear algebra templates for iterative methods[END_REF], there is about 20% speed-up compared to PETSc at such scales. Nevertheless, we notice that for 2,016 MPI processes PETSc is significantly faster than ECG. This is likely because the number of iterations with PETSc is reduced with respect to 1,008 MPI processes. This behavior is not expected because it is known that block Jacobi preconditioners are not scalable (see [START_REF] Dolean | An introduction to domain decomposition methods[END_REF] for instance). Indeed, we observe that the number of iterations is D-Odir [START_REF] Dolean | An introduction to domain decomposition methods[END_REF] PETSc effectively increasing both for ECG and CG when the number of MPI processes increases.

Then we make the same comparison on Ela 30 test case for which we use D-Odir and t = 24, as discussed previously. The resulting runtimes are summarized in Table 5. We observe that both PETSc and ECG are scaling very well. Enlarging the Krylov subspaces allows us to reduce drastically the number of iterations: D-Odir(24) performs around 25 times less iterations than CG. As a consequence, D-Odir [START_REF] Kalhan | Reverse communication interface for linear algebra templates for iterative methods[END_REF] is more than 5 times faster than PETSc PCG at small scale and around 2.5 times faster at large scale. We believe that this relatively poor scaling of D-Odir(24) compared to PETSc PCG at large scale is due to the implementation that is not as optimized as PETSc which has been developed for many years. For instance, the routine we use for computing the sparse matrixset of vectors multiplication is certainly not as optimized as that of PETSc for computing the sparse matrix-vector multiplication. Also, we mentioned that we are currently performing 4 calls to MPI_Allreduce per iteration, but that could be reduced to 2 by fusing them. Furthermore, we could use Pipelining [START_REF] Ghysels | Hiding global synchronization latency in the preconditioned conjugate gradient algorithm[END_REF] or Communication-Avoiding based on s-step methods [START_REF] Hoemmen | Communication-avoiding Krylov subspace methods[END_REF][START_REF] Carson | Avoiding communication in nonsymmetric lanczos-based krylov subspace methods[END_REF]

Dependence on the mesh size

Given the importance of the parameter t regarding the efficiency of the method, we perform a study of the convergence of the method with respect to the mesh size for the elasticity test case. More precisely, we consider the mesh used for generating the Ela 30 matrix, then we coarsen it by dividing the number of points in each dimension by 2 1/3 . Thus, we generate 3 additional elasticity matrices on which we perform a weak scaling experiment. Our major focus is not the weak scaling of ECG, but rather the comparison between PETSc's CG and D-Odir in terms of runtime.

The results are summarized in Table 6. First, we fix t = 24 and we perform a weak scaling study. We observe that D-Odir [START_REF] Kalhan | Reverse communication interface for linear algebra templates for iterative methods[END_REF] is always between 2.5 to 2.9 times faster than PETSc PCG, but the gap tends to slightly decrease when the number of MPI processes increases. As ECG(t) is acting as if the t smallest eigenvalues of the matrix were deflated, it seems natural to use smaller values of t for the smaller matrices. Indeed, we perform another set of experiments where we vary t as well as the size of the problem. We observe that this slightly improves the results for the smaller matrices, we obtain for example an overall speed-up of 3.1 compared to 2.9 for the smallest matrix. As expected, this does not have a significant effect on the larger matrices. Another interesting observation is that even if the number of iterations is almost constant when increasing both t and the number of MPI processes, it does not improve the scaling because the cost of one iteration also increases. Overall, D-Odir is still at least 2.5 times faster than PETSc CG.

Impact of threads on performance

One motivation for enlarging the Krylov subspaces is to increase the arithmetic intensity of the resulting methods. This is particularly interesting to take advantage of the so-called manycore architecture as Nvidia GPUs, Intel Xeon Phi, Omin [START_REF] Agullo | Block gmres method with inexact breakdowns and deflated restarting[END_REF] Odir [START_REF] Kalhan | Reverse communication interface for linear algebra templates for iterative methods[END_REF] # omp time (s) speed-up time (s) speed-up or Sunway SW26010 used in the Sunway TaihuLight supercomputer. As the implementation relies on the MKL library which is multi-threaded [START_REF] Wang | Intel math kernel library[END_REF], it is straightforward to assess its efficiency on the Xeon Phi processors.

In order to do so, we perform the following experiments on NERSC's supercomputer Cori. It consists in two partitions, one with Intel Haswell processors and another one with the last generation of Intel Xeon Phi processors: Knights Landing (KNL). More precisely, the second partition consists in 9,688 singlesocket Intel Xeon Phi 7250 (KNL) processors with 68 cores each. For a detailed description of the machine, we refer to the online documentation2 . We compile the code (and its dependencies) using the default compilers and libraries installed on the machine: icc version 18.0.1, cray-mpich version 7.6.2, MKL version 2018.1.163 and METIS version 5.1.0. We consider Ela 30 test case and we study the impact of threads on the strong scaling of Odir [START_REF] Kalhan | Reverse communication interface for linear algebra templates for iterative methods[END_REF]. We do not use the dynamic reduction of the search directions in order to keep the cost of one iteration constant during the solve to better understand the effect of threading. We both increase the number of MPI processes from 64 to 2048 and the number of threads from 1 to 8 -this means at most 2, 048 × 8 = 16, 384 threads, each one being bound to one physical core.

The results obtained are summarized in Fig. 5a. First of all, we notice that there is a trade-off between using threads or MPI processes because the number of MPI processes dictates the preconditioner. Indeed, there are as many blocks in the block Jacobi preconditioner as the number of MPI processes, thus increasing the number of MPI processes deteriorates the quality of the preconditioner. For instance, using 512 MPI processes takes 123s, and using 64 MPI processes with 8 threads each takes 179s: it is an increase of 50% compared to 512 MPI processes. Nevertheless, we observe that using more than 2 threads, and up to 8, has always a significant effect on the speed-up, even when the number of MPI processes is high. For instance, as shown in Table 5b, increasing the number of threads from 1 to 8 with a fixed number of 2, 048 MPI processes leads to an decrease in runtime of nearly 3. Of course, we are not close to full efficiency when using multiple threads, but we are still taking advantage of the BLAS 3 routines. This is illustrated by the Table 5b where we compare the speed-up obtained by using threads for Omin [START_REF] Agullo | Block gmres method with inexact breakdowns and deflated restarting[END_REF], which corresponds to the classical PCG, and Odir [START_REF] Kalhan | Reverse communication interface for linear algebra templates for iterative methods[END_REF]. We observe that using more than 2 threads is not effective at all with Omin(1) whereas it always significantly increases the speed-up with Odir [START_REF] Kalhan | Reverse communication interface for linear algebra templates for iterative methods[END_REF].

Finally, we are able to obtain an overall speed-up of 50 when using 16, 384 cores with respect to 64 cores. Compared to an ideal speed-up of 256, it may seem that this result is not very good (around 20% of efficiency), however it is well-known that Krylov methods may face efficiency issues at very large scale 3 -in practice, such difficulties are overcome by using preconditioning strategies well adapted to the underlying problem. Furthermore, it is important to notice that the matrices tested are relatively small, but they allow us to simulate extreme scale situation: with 16, 384 cores the average number of unknowns per core is around 280. We have shown that in such cases, using enlarged Krylov subspaces allows to increase arithmetic intensity while decreasing the communication by drastically decreasing the overall iterations. Thus it takes advantage of the current trend in hardware architecture for reaching exascale.

Conclusion

In this paper, we have studied the Enlarged Conjugate Gradient method. It relies on so-called enlarged Krylov subspaces which can be seen as particular cases of block Krylov subspaces. The parallel efficiency of the approach has been assessed, and we have shown that this method is scalable up to 16, 384 cores and it is up to 5.7 times faster than PETSc's implementation of PCG.

First, we have thoroughly studied the method from a theoretical point of view. In particular, we have described the link between the two variants of the method (Orthodir and Orthomin), and thus explained their differences in terms of robustness. We also have studied its convergence rate and we have showed that the Enlarged Conjugate Gradient method is acting as if the smallest eigenvalues of the matrix were somehow deflated. Moreover, we have presented dynamic versions of Orthodir and Orthomin where the number of search directions is adaptively reduced during the iterations. Starting from the theory, we have described the link between the two variants of the method, and thus explained their differences in terms of robustness. Then, we studied its convergence rate and we have showed that ECG is acting as if the smallest eigenvalues of the matrix were somehow deflated. Then, we have described the parallel design of the method, including the two variants as well as their dynamic versions. Numerical experiments show that enlarging the Krylov subspaces allow to reduce significantly the number of iterations with respect to the classical PCG method. Furthermore, the reduction of the search directions allows to reduce the cost of the extra arithmetic operations induced by the method. Overall, the proposed solver is up to 5.7 times faster than PETSc PCG for an elasticity matrix. Also, as our implementation is based on BLAS 3 kernels only, it offers a good scaling when increasing the number of threads per MPI process. Thus, it is scaling up to 16, 384 threads, each one binded to a physical core, and it seems well adapted for so-called manycore architectures.

Throughout this work, the only assumption we make is that the matrix is symmetric positive definite. Hence the resulting methods are very generic and completely algebraic. For instance, it is straightforward to use D-Odir for solving linear systems with several right-hand sides. Similarly, ECG can be used with any preconditioner that could be used with CG. Thus, it can be used as a blackbox solver that can be integrated very easily in any existing code. As it increases the arithmetic intensity and reduces the communication, it is well suited for modern and future architectures that exhibit massive parallelism. According to the theoretical study and the numerical experiments, it is particularly well adapted for matrices with a small number of low eigenvalues.

0 3 : k = 1 4 :

 34 for k = 1, . . . , k max do 5:

0 3 : k = 1 4 :

 34 for k = 1, . . . , k max do 5:

Figure 2 :

 2 Figure 2: Local distribution of the data: Orthodir on the left and Orthomin on the right.

Figure 3 :

 3 Figure 3: Heterogeneity pattern of the Young's modulus and Poisson's ratio for elasticity matrices.

Figure 4 :

 4 Figure 4: Convergence of D-Odir (odir-1) compared to Odir (odir-0). The dash line represents the number of search directions for D-Odir. In parenthesis the difference of iteration count to reach convergence between D-Odir and Odir (+ means that D-Odir took more iterations to converge).

  (w.r.t. 64 MPI/1 OpenMP) (a) The bars represent the runtime (left) and the lines represent the corresponding speed-up with respect to 64 MPI with 1 thread each (right).

Figure 5 :

 5 Figure 5: Strong scaling study for Ela 30 matrix on Cori (# omp stands for the number of threads assigned to each MPI processes).

Table 1 :

 1 Complexity and memory consumption of Orthodir, Orthomin and CG where t is the enlarging factor, n is the number of rows of A and P is the number of processors. In parenthesis is indicated the number of calls to MPI_Allreduce.

	sparse matrix set of vectors
	form Z k AZ k
	Cholesky factorization
	update P k and AP k

Table 2 :

 2 Test matrices.

	Name	Size	Nonzeros	Problem
	Flan 1565	1,564,794 117,406,044	Structural problem
	Bump 2911 2,911,419 130,378,257 Reservoir simulation
	Ela 20	2,118,123 74,735,397	Linear elasticity
	Ela 30	4,615,683 165,388,197	Linear elasticity

Table 3 :

 3 Runtime results (in seconds) for Flan 1565 (P = 56), Ela 20 (P = 112) and Bump 2911 (P = 112). The ++ means that the maximum number of iterations (5000) was reached and the -means that a breakdown occurred.

	RR n°9190

Table 4 :

 4 Strong scaling study for Flan 1565. The speed-up is the ratio between PETSc runtime and ECG runtime.

	CG

Table 5 :

 5 Strong scaling study for Ela 30. The speed-up is the ratio between PETSc runtime and ECG runtime.

Table 6 :

 6 Weak scaling study. The dimension of the matrix is denoted n, and t denotes the enlarging factor. The speed-up (sp.-up) is the ratio between PETSc runtime and ECG runtime.

	on top of ECG

RR n°9190

https://www.hpc2n.umu.se/resources/hardware/kebnekaise RR n°9190

http://www.nersc.gov/users/computational-systems/cori/configuration/

This is well illustrated by HPCG benchmark: http://www.hpcg-benchmark.org/custom/ index.html?lid=155&slid=293 RR n°9190

RR n°9190 RESEARCH CENTRE PARIS -ROCQUENCOURT Domaine de Voluceau, -Rocquencourt B.P. 105 -78153 Le Chesnay Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgments

The computational resources were provided by the High Performance Computing Center North (HPC2N) at Umeå (Sweden), and the National Energy Research Scientific Computing Center (NERSC) at Berkeley (US), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

The work of this author is funded by the NLAFET project as part of European Union's Horizon 2020 research and innovation programme under grant agreement No 671634).

Mots-clés : Méthodes des sous-espace de Krylov, Gradient Conjugué, Algorithmes qui réduisent les communications


Romaric Guillerm
email: guillerm@laas.fr

Hamid Demmou
email: demmou@laas.fr

Safety Management Method in Complex System Engineering

Keywords: Safety requirement, Requirement engineering, Complex system I

The main objective of System Engineering is the successful development of complex system. It is based on the application of iterative and recursive processes on each phase or step of the system development. One critical process is the requirement management, particularly when it deals with the safety requirements. These one are non-functional requirements and are related to emergent properties, which come from the integration of the different system components. They must be identified as soon as possible, because they are guards to validate or not the system, which can require changes in system architecture. Moreover, they are formulated at system level and need to be derived at sub-system level.

The objective of this paper is to propose a safety management method based on well-known safety methods, in order to organize the different tasks to make the system safe. The method focuses mainly on the definition of the system safety requirements following risk and hazard analysis, and also on their derivation according to a top-down approach. It is based on the well-known Failure Mode, Effects, and Criticality Analysis (FMECA) and the use of Fault Trees and Event Trees.

INTRODUCTION

Modern systems are increasingly complex [START_REF] Chavalarias | French roadmap for complex systems 2008-2009[END_REF]. Indeed, they integrate more and more different technologies, offer more functions, and have complex interactions between their components. The processes and the design methods must evolve to reflect this growing complexity [START_REF] Juristo | Is the european industry moving toward solving requirements engineering problems[END_REF], [START_REF] Komi-Sirvio | Great challenges and opportunities of distributed software development -an industrial survey[END_REF]. In particular, for our purposes, the management of properties such as reliability or security [START_REF] Avizienis | Basic concepts and taxonomy of dependable and secure computing[END_REF] must evolve accordingly, to ensure and enable the necessary level of confidence [START_REF] Rasmussen | Risk management in a dynamic society: A modelling problem[END_REF]. For an effective consideration of safety in the design process, it is necessary to consider safety in overall studies by using the engineering system processes. The safety properties must be defined globally; that is to say elicited [START_REF] Goguen | Techniques for requirements elicitation[END_REF]. Once these safety properties are identified, they must be derived locally to be actually realized by the system. The local properties associated with subsystems must be satisfied to ensure the global properties, reaching issues of traceability [START_REF] Gotel | An analysis of the requirements traceability problem[END_REF], [START_REF] Sahraoui | Requirements traceability issues: Generic model, methodology and formal basis[END_REF] and requirements engineering [START_REF] Sommerville | Software engineering (update)[END_REF].

Requirements Engineering (RE) is one of the System Engineering (SE) processes. RE is a crucial process within the development of complex system. Safety requirements are classified as non-functional requirements and are related to emergent system properties. They cannot be attributed to a single system component. Furthermost, non-functional requirements are fundamental to determine the success of a system. Two activities are defined in RE. The first one concerns requirements development including the processes of elicitation, documentation, analysis and validation of requirements. The second one concerns requirement management which includes the processes of maintainability management, changes management and requirements traceability.

The work presented in this paper concerns a part of our approach for the integration of safety in system engineering processes [START_REF] Guillerm | System engineering approach for safety management of complex systems[END_REF]. It is an improvement and extension of the method presented in [START_REF] Guillerm | Combining fmeca and fault trees for declining safety requirements of complex systems[END_REF], that was inspired from [START_REF] Lindsay | Derivation of safety requirements for an embedded control system[END_REF] with an engineering process and requirements point of view. The originality of the approach is to take into account the safety requirements in system engineering process to facilitate traceability of these requirements throughout the life cycle of the system. It concerns the two activities of RE: the development and the management activities. The paper presents a method that allows to define, derive and derive system safety requirements, with the combination of several FMECA (Failure Mode, Effects, and Criticality Analysis) [START_REF] Buzzatto | Failure mode, effects and criticality analysis (fmeca) use in the federal aviation administration (faa) reusable launch vehicle (rlv) licensing process[END_REF], Fault Trees analysis [START_REF] Lee | Fault tree analysis, methods, and applications -a review[END_REF] and Event Trees analysis [START_REF] Villemeur | Sûreté de fonctionnement des systèmes industriels[END_REF]. This paper contains five parts. The second one presents the system engineering framework of the method. The third one exposes the method for safety requirement definition and derivation, with its different steps. The fourth one presents briefly the result of one first case study. Finally, the last section concludes the paper and presents some perspectives.

II. CONTEXT

In this part, the context of our work is exposed. The first section presents the System Engineering notion. Then, the SE standard that we adopt is presented with its useful concept of building block that split the design into different system layers. To finish, a focus is done on the safety requirement management.

A. System Engineering

System Engineering (SE) is an interdisciplinary approach, whose objective is to assist the development of new systems. It contains collaborative and interdisciplinary processes of resolution of problems, supporting knowledge, methods and techniques resulting from the sciences and experiment to define a system, which satisfies an identified need, and is acceptable for the environment, while seeking to balance the total economy of the solution, on all the aspects of the problem in all the phases of the development and the life of the system. SE concepts are adequate specifically for complex problems [START_REF] Sahraoui | Issues in systems engineering research[END_REF]. SE is the application of scientific and engineering efforts to:

-Transform an operational need into a description of system performance parameters and a system configuration, through an iterative process of definition, synthesis, analysis, design, test and evaluation.

-Integrate reliability, safety, maintainability, expandability, survivability, human engineering and other factors into the total engineering effort to meet cost, schedule, supportability and technical performance objectives.

SE is the global framework of the approach proposed in this paper.

B. EIA-632

A standard currently used in the industrial and military fields is the EIA-632 standard [START_REF]Processes for engineering systems. Electronic Industries Alliance standard[END_REF]. Our work is also based on it. Briefly, this standard covers the product life cycle from the needs capture to the transfer to the user. It is constituted by 13 processes grouped into 5 sets (see Figure 1):

1. Technical management processes (three processes): these processes monitor the whole process ranging from the initial idea to building a system until the delivery of the system.

2. Acquisition and supply processes (two processes): these processes ensure the supply and acquisition (and are very close to logistics).

3. System design processes (two processes): these processes deal with the elicitation and the acquisition of requirements and their modeling, the definition of the logical design and its physical solution.

4. Product realization processes (two processes): these processes deal with the implementation is-sues of the system design and its use.

5. Technical evaluation processes (four processes): these processes deal with verification, validation and testing issues.

C. Building Block Concept

The EIA-632 standard adopts an original and interesting system decomposition based on the concept of "building block". A building block is the association between one (or several) final product and a set of enabling products, as shown in Figure 2. In fact, the system is seen as a hierarchy of building blocks. The solutions defined in the upper layer (level) blocks, described by a set of specified requirements, are allocated as input requirements for the lower layer blocks (see Figure 3). Finally, the building block decomposition is stopped when blocks correspond to on-the-shelf components or when their realization can be subtracted. With this description, we identified the need of deriving the safety requirements through the hierarchical decomposition.

D. Safety Requirement Management

To situate the position of our method for deriving safety requirements, the Figure 4 gives an overview of the involved EIA-632 system engineering processes.

Among the different possible sources of safety requirement we can find the requirement provided by some dependability analysis as shown in the Figure 4. In this paper we consider this source of requirements. The proposed approach is used to define and derive safety requirement with different safety analysis.

III. THE SAFETY MANAGEMENT METHOD

In this second part, the safety management method is presented. First, an overview of the method is given, followed by an explanation about the different kinds of safety requirements that are taken into account in the current version of the method. Afterwards, the 9 steps of the method are explained in details.

A. Overview

The method assumes that a complex system is composed of some subsystems (the principle of Building Block of the EIA-632 standard). It combines FMECA, fault trees and event trees, and has the objective to define all safety requirements at system level and to derive them locally at subsystems level with a goal of traceability. The Figure 5 summarizes the process associated to the method and illustrates how the different steps are integrated together.

B. Classification of the considered Safety Requirements

The method enables to identify and deals with several kinds of safety requirements. We classify these requirements into subcategories, which are:

-Reliability requirement, that claims a quantitative objective in term of reliability properties.

-Architectural requirement, that defines an architectural design to deal with safety (like redundancies).

-Active functional security requirement, that is related to additional security equipment (protective barrier) that can participate to reduce the probability of an accident.

-Passive functional security requirement, that is related to an additional protective or mitigation equipment that can reduce the severity of an accident.

C. Step 1: Risks Identification

The first step is to identify and classify all the system risks. These can be human actions, external failures, internal failures (those of the system) or environmental conditions. The classification must be done in two groups: the risks representing internal failure modes and the other risks.

D. Step 2: System Failure Modes Analysis

The second step is to begin the analysis of risks that correspond to system failure modes. The recommended method is the FMECA [START_REF] Buzzatto | Failure mode, effects and criticality analysis (fmeca) use in the federal aviation administration (faa) reusable launch vehicle (rlv) licensing process[END_REF], that is a technique used to identify, prioritize, and eliminate potential failures from a system, a design or a process. Concretely, this step is to complete few columns of the FMECA table (others than severity, probability, criticality and corrective action) (see Table I). For each system function, we identify failure modes, causes of these modes and effects on the system (possibly depending on the phase, state or mode). For the identification of failure modes, lists of generic modes have been defined in some standards like CEI 60812 [START_REF]Techniques d'analyse de la fiabilité des systèmes[END_REF]. The effects are here the potential accidents. In fact, we also propose some changes in the classical FMECA to clarify the method, (see Table I). A distinction is made between the probability and the detectability of the failure modes and those of the effects. Indeed, between a failure mode and an effect (accident), there is a set of involved "cofactors" (protection barrier, environmental condition ...), recorded in the "condition" column of the table. These conditions will be identified during the step of consequences analysis.

The assessment of the probability of the risk and the assessment of the severity, probability and criticality of the effect will be done during the step of risk assessment. The corrective actions will be proposed at the risk mitigation step.

E. Step 3: Other Risks Analysis

This step is similar to the previous step of system failure modes analysis, but focuses on the other risks (external). The recommended method is to use the principle of an FMECA, that we can call here RECA (Risks, Effects, and Criticality Analysis) (see Table II). This step is to complete few columns of the RECA table (others than severity, probability, criticality and corrective action). The effects are also the potential accidents. The same remarks as for the FMECA remain true concerning the probability, the detectability and the severity of the failure modes and the effects, and the conditions.

F. Step 4: Consequences Analysis

In this step, the consequences of all the identified risks (system failure modes and others) must be analyzed. It is to identify how the risks contribute to an accident. It can be done using event trees [START_REF] Villemeur | Sûreté de fonctionnement des systèmes industriels[END_REF] to visualize the possible chains of events that led from the risk to the accident, through branching points representing protective measures or interventions (cofactors) (see Figure 6). The minimal cuts associated with the various accidents are also identified. A generic example is given in Figure 6. The effects correspond to accidents, whereas the term consequence refers to events or factors involved in the causes to consequences relationship, starting from the analyzed risk.

G. Step 5: Causes Analysis

The fifth step is to conduct an internal analysis of the system by identifying the causes of system failure modes. These causes analysis must lead to subsystems failure modes. For this step, the use of fault trees [START_REF] Lee | Fault tree analysis, methods, and applications -a review[END_REF] is recommended. Indeed, a fault tree provides a simple modeling way to represent the interactions between components from the point of view of reliability. Static fault trees use traditional Boolean logic functions to represent the combination of component failures (events) that cause system failure.

So, the top event of each tree corresponds to a system cause. The objective is to determine the causes of the top event (using logical operators such as AND and OR) in the subsystems. The leaves of the fault tree correspond to sub-systems failure modes (see Figure 7). In fact, the analysed system failure modes correspond either directly to a system risk (defined in the first step), or to a cofactor of an event tree which is a system failure mode.

H. Step 6: Sub-systems Failure Modes Analysis

An analysis of the subsystems failure modes should be done in parallel, using FMECA. The subsystems failure modes used in step 5 re-appears (the principle of the FMECA analysis). This FMECA will define the corrective actions at the sub-systems level that are representative of subsystem reliability requirements.

I. Step 7: Risks Assessment and Mitigation

The seventh step is the central one. It deals with risks assessment and risks mitigation with definitions of corrective actions.

1) Assessment

The risk assessment consists in defining the severity and the probability of the identified accidents, in order to evaluate the criticality. This information must be recorded in the various FMECA and the RECA tables. Concerning the probability, this one must be evaluated based on the fault trees and the event trees. Finally, it must be decided whether the risks are acceptable or not.

2) Mitigation

The risk mitigation consists in advocating corrective actions (to be filled in the FMECA and the RECA) for the risks qualified as "non-acceptable" during the risk assessment step, in order to make them become "acceptable".

The corrective actions can:

-Reduce the probability of the accident, by:

• Fixing an objective of reliability with a reliability requirement (at system or subsystem level).

• Modifying the system architecture for a better reliability (with redundancies for example) with an architectural requirement, that derives from the reliability requirement of the objective.

• Adding additional security equipment (protective barrier) with an active functional security requirement, that derives from the reliability requirement of the objective. During the next iteration of the method, reliability requirements will be defined for this security equipment based on the analysis of the failure modes in which it participates.

-Try to satisfy a criterion, for example:

• A single failure criterion, adding a security equipment (barrier) to increase the necessary number of failures that lead to the occurrence of an accident, with an active functional security requirement.

• A spatial dispersion criterion, with an architectural requirement.

• A redundancy with separate development criterion, with an architectural requirement.

-Reduce the severity of the accident, by:

• Adding a protection or mitigation equipment, with a passive functional security requirement.

Note: This is not the only possible corrective actions (preventive maintenance for example). Other types of corrective actions will be incorporated in future work to improve the process.

J. Step 8: Safety Requirements Synthesis

Before eventually transferring the change requests to modify the system, this step will summarize the results in terms of requirements, derivation of requirements, and traceability links between requirements and between requirements and accidents. As in the first version of the method [START_REF] Guillerm | Combining fmeca and fault trees for declining safety requirements of complex systems[END_REF], the derivation part is based on the following 3 types of relations:

-System causes and system corrective actions, -System causes and sub-systems failure modes, -Sub-systems failures modes and sub-systems corrective actions.

A generic example of this synthesis is given in Figure 8.

K. Step 9: Stop Criterion

The process is finished once all the risks are considered as "acceptable". If this is not the case, a change request must order to modify the system and the method must be reapplied from the beginning by updating the different analysis.

IV. APPLICATION ON AN EXAMPLE

In order to validate the potential benefit of the proposed method, we have applied it on a short running example. It is the "lift" system with the following initial considered subsystems: "cable" and "motion controller". The main considered function is to "bring up or down the passenger of the cabin".

A. Step 1: Risks identification

Two risks are identified and considered:

-"Cabin fall", as a system failure mode, -"Fire", as an environmental condition.

B. Step 2: System Failure Modes Analysis

The "cabin fall" risk appears in the FMECA of the system because it corresponds to an internal failure mode. The final version of the FMECA table is given in Figure 9, where all columns are already fielded (which is normally not the case at the beginning of the study). We can notice the difference that we do in our new version of FMECA between the probability of the failure mode and the probability of the accident (which is the effect). This difference can be justify with the event tree of the step 4.

C. Step 3: Other Risks Analysis

The second risk, "fire", appears in the RECA of the system. One more time, it is the final version of the RECA after several iterations of the methodology that is given in Figure 10, where all columns are already fielded.

D. Step 4: Consequences Analysis

Then, the event tree of the Figure 11 exposes the result of the consequences analysis. The "cabin fall" failure mode can have the "crash of the cabin" accident as consequences, if the emergency brake is out of service. This emergency brake was added after an iteration of the methodology in the step 7 -risks assessment and mitigation -in order to reduce the probability of the accident caused by the "cabin fall" risk. About the consequences of the "fire" risk, a second event tree, Figure 12, shows that it brings to an accident, but with different severity depending if the fire extinguisher is operational or not. This fire extinguisher was added after an iteration of the methodology in the step 7 -risks assessment and mitigation -in order to reduce the severity of the accident caused by the "fire" risk.

E. Step 5: Causes Analysis

In the Figure 13, the fault tree of the "cabin fall" failure mode shows the result of the causes analysis. The "cabin fall" happens if there is an undesirable descent order from the motion controller or if the two cables are broken together. The Figure 14, 15 and 16 show respectively the FMECA of the "emergency brake", "motion controller" and "cable" sub-systems. These sub-systems are the one that are involved in the causes of the "cabin fall" failure mode or in the consequences of the two risks, "cabin fall" and "fire". The failure modes that appears in these sub-system FMECA corresponds to leaves in the fault tree or to cofactors in the event trees. Fig. 14. FMECA of the "emergency brake" sub-system Fig. 15. FMECA of the "motion controller" sub-system Fig. 16. FMECA of the "cable" sub-system

G. Step 7: Risks Assessment and Mitigation

Finally, we succeed to make these two risks acceptable after two iterations of the method, by: -Modifying the system architecture with a redundancy of two cables, in order to reach the reliability for the "cabin fall" system failure mode.

-Adding an "emergency brakes" sub-system as one active security equipment, in order to satisfy the single failure criterion.

-Adding a "fire extinguisher" sub-system as one passive security equipment, in order to reduce the severity of the "fire" environmental condition.

H. Step 8: Safety Requirements Synthesis

The final safety requirements synthesis for this case study is presented in figure 17, where appears all the necessary requirements and there traceability links with the risks. The "cabin fall" risk is treated by two main system safety requirements which are "crash of the cabin must have a frequency < 1E-9/h" and "must have 1 emergency brake".

The first one reduces the probability of this risk and is derived into three safety requirements at subsystem level:

-"Cable break must have a frequency < 1E-7/h", associated to the "cable" subsystem, -"Controller order stuck in the down state must have a frequency < 1E-7/h", associated to the "motion controller" sub-system, -"Inoperative emergency brake must have a frequency < 5E-3/h", associated to "cable" sub-system.

Another safety requirement derive from this first system safety requirement, which is "must have 2 cable in redundancy", in order to reach the reliability level of the failure mode. It corresponds to a modification of the architecture of the system with a redundancy.

The second one -"must have 1 emergency brake"participates also to the reduction of probability of the accident caused by the failure mode. But mainly, it was added to satisfy the single failure criterion, that's way it is linked directly with the risk, and not as a derivation of the first system safety requirement.

Concerning the "fire" risk, in this example, it is treated by only one system safety requirement: "must have 1 fire extinguisher". The objective of this requirement is to reduce the severity of the effect of this risk.

V. CONCLUSION

The method provides a support framework to define system safety requirements with an objective of traceability and requirements management. The interest is multiple for the safety field: the method deals with the safety elements (failure modes, safety requirements...) and it is done with a comprehensive system engineering (with traceability and requirements derivation) which is a factor contributing to safer systems. This method is compatible with the standard EIA-632 [START_REF]Processes for engineering systems. Electronic Industries Alliance standard[END_REF], and it extends the principle and strengthens the links between failure modes researches and analysis (FMECA), causes analysis and effects analysis.

In this work, several safety attributes are taken into account, like reliability, passive security and active security. They correspond to the given classification of safety requirements, which are themselves defined from the corrective actions. Other requirements concerning maintainability or availability should also be considered in further study. The probability, the severity and the criticality was treated through the FMECA.

However, the work still doesn't consider the detectability

Fig. 1 .

 1 Fig. 1. EIA-632 Standard System Engineering Processes

Fig. 2 .

 2 Fig. 2. One Building Block

Fig. 3 .

 3 Fig. 3. Multilayer Building Block

Fig. 4 .

 4 Fig. 4. Dependability analysis as a source of requirements

Fig. 6 .

 6 Fig. 6. Event tree

Fig. 7 .

 7 Fig. 7. Fault tree

Fig. 8 .

 8 Fig. 8. Requirements traceability synthesis

Fig. 9 .

 9 Fig. 9. FMECA of the "lift" system

Fig. 10 .

 10 Fig. 10. RECA of the "lift" system

Fig. 11 .

 11 Fig. 11. Event tree of the "cabin fall" failure mode

Fig. 12 .

 12 Fig. 12. Event tree of the "fire" risk

Fig. 13 .

 13 Fig. 13. Fault tree of the "cabin fall" failure mode F. Step 6: Sub-systems Failure Modes Analysis

Fig. 17 .

 17 Fig. 17. Requirements traceability synthesis for the "lift" example

TABLE I .

 I FMECA

TABLE Fig. 5. Overview of the safety management method

TABLE II .

 II RECA TABLE

aspect. We also should update the tool "SafetyLab" (http://homepages.laas.fr/rguiller/) that implements the first version of the method presented in [11].