Keywords: numerical analysis, domain decomposition methods, algebra, scientific computation . . . .

 plus methods are used to build the new Hybrid-FETI. The basic idea is to develop a new algorithm that can use both methods at the same time by choosing in each interface the most suited condition depending on the characteristics of the problem. By doing this we search to have a faster and more robust code that can work with configurations that the base methods will not handle it optimally by himself. The performance is tested on a contact problem. The following part involves the development of a new implementation for the S-FETI method [39], the idea is to reduce the memory usage of this method, to make it able to work in larger problem. Different variation for this method are also proposed, all searching the reduction of directions stored each iteration of the iterative method. Finally, an extension of the FETI-2LM method to his block version as in S-FETI, is developed. Numerical results for the different algorithms are presented.

Z or G owned by subdomain Ω (s) , in black the interfaces where F times the red modes is also non null. Interface point of view (Right):

In doted lines the coarse modes describing FZ or FAG owned by local interface Γ (sq) between subdomains Ω (s) and Ω (q) , in black the modes where, due to Γ (sq) , Z T (FZ) is non null. . . . . . . . . 2.6 One-way split of six subdomains . . . . . . . . . . . . . . . . .

Introduction

In the last decades and thanks to the increasing computational power, faster robust and accurate algorithms had been developed to solve numerically a large variety of problems modeled by Partial Differential Equations (PDE). The use of multiple processors to increase the speed of calculations leads to the search of strategies in parallelism that allows to profit of this new computers architectures. Different parallel iterative and direct methods for solving linear systems had been developed [START_REF] Sameh | On Stable Parallel Linear System Solvers[END_REF], [START_REF] Duff | Parallel implementation of multifrontal schemes[END_REF], [3], both with positive and negative features in terms of speed, memory and accuracy. Results of iterative methods based on Krylov spaces usually depend on the condition number of the matrix representing the system, and the memory requirements are usually not an issue, on the other hand direct methods are more robust but the use of memory can be a problem for some large systems.

For problems coming from discretization of Finite Element Methods we have properties that allow us to use a different approach, that is the Domain Decomposition Methods (DDM), this methods can be considered as an hybrid between iterative and direct methods [START_REF] Quarteroni | Domain Decomposition Methods for Partial Differential Equations[END_REF], [START_REF] Dolean | An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation: Other Titles in Applied Mathematics[END_REF]. They are based on the partition of the domain of the problem into subdomains, where smaller system of equations are defined. From this division, this methods can be categorized into two large groups, the overlapping and non overlapping methods.

In this work we will focus in some DDM with non overlapping interfaces, mainly the method of Finite Elements Tearing and Interconnecting or FETI and other related methods [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF], [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method-part I: A faster alternative to the two-level FETI method[END_REF], [START_REF] Farhat | Two-level domain decomposition methods with Lagrange multipliers for the fast iterative solution of acoustic scattering problems[END_REF].

The objective of our work is to develop new FETI methods, that apply to particular cases, improving the results shown by the existing methods. Also we will expand the results on one of the existing FETI methods to extend his 1 applications to problems in which the current formulations does not allows it.

Contributions of this thesis

The following work is based in one of the most known of the non overlapping domain decomposition methods, the Finite Element Tearing and Interconnecting method (FETI) [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF], [START_REF] Farhat | Optimal convergence properties of the FETI domain decomposition method[END_REF].

The works on this method since its first formulation have produced several improvements to the same [START_REF] Farhat | The two-level FETI method Part II: Extension to shell problems, parallel implementation and performance results[END_REF], [START_REF] Daniel | A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems[END_REF], but also it has permitted the formulation of new FETI methods [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method-part I: A faster alternative to the two-level FETI method[END_REF], [START_REF] Farhat | Two-level domain decomposition methods with Lagrange multipliers for the fast iterative solution of acoustic scattering problems[END_REF]. Within this context, and using the similarities on the constructions of two of the most used FETI methods, namely the original FETI-1LM and the later developed FETI-2LM. We will formulate a new algorithm based in this two methods, that tries to take advantage of the good properties of one and the other. After the development of this new method, and after a better understanding of the same we can exhibit his advantages in contact cases where it outperforms the FETI-1LM and FETI-2LM methods.

In a different line of work, this time following the development of the new S-FETI method [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF], we will continue the analysis within, in order to find new, faster or more robust forms of the same. Different variations of it will be formulated and tested, all of them trying to improve the existing results shown by the method. We will also extend the application of S-FETI to a larger class of problems, with a new implementation, base on a sparse storage, that reduces the memory limitations of this method.

Next we will try to use the ideas presented in the formulation of S-FETI to develop new FETI algorithms with certain improvements, however since they are in a basic stage they present several issues that leads to new sources of research.

Finally, we use the same idea that led to the S-FETI formulation but this time is applied to the FETI-2LM method, in order to develop a new Block version of this last one.

Before giving more details about FETI and the other FETI-like methods, we want to recall some of the basic linear algebra tools needed in this thesis to understand this type of algorithms, we refer to the iterative solvers for linear Lemma 0.2. If A p g 0 ∈ K p then A p+q g 0 ∈ K p for every q > 0 Proof. By induction. For q ≥ 0, if we have A p+q g 0 ∈ K p then A p+q g 0 = p-1 k=0 α k A k g 0 and therfore

A p+q+1 g 0 = p-2 k=0 α k A k+1 g 0 + α p-1 A p g 0 = p-2 k=0 α k A k+1 g 0 + α p-1 p-1 k=0 β k A k g 0 = p-1 k=0 γ k A k g 0 (3)
Lemma 0.3. The Krylov space succession its strictly increasing from 1 to p max then it stagnates from p = p max Proof. If p is the smallest integer that makes A p g 0 dependent of previous vectors, then the vectors (g 0 , Ag 0 , A 2 g 0 , . . . , A p-1 g 0 ) are linearly independent and K q has a dimension of q, for every q ≤ p. In particular K p has a dimension p.

Furthermore, A p g 0 ∈ K p and, from Lemma 0.2, every vector A p+q g 0 is in K p , for every q > 0, which implies that K p+q = K p for every q > 0.

We then have: K 1 ⊂ • • • ⊂ K p = K p+q for every q > 0. And by definition of p max , we have that p = p max Theorem 0.4. The solution of the linear system Ax = b is in the affine space x 0 +K p max Proof. From Lemma 0.2 and Lemma 0.3 the vectors (g 0 , Ag 0 , A 2 g 0 , . . . , A p max -1 g 0 ) are linearly independent and

A p max g 0 = p max -1 k=0 α k A k g 0 ( 4
)
In this equation, the coefficient α 0 is non null, from which, multiplying doth terms by A -1 we obtain A p max -1 g 0 = p max -1 k=0 α k A k-1 g 0 [START_REF] Axelsson | A Black Box Generalized Conjugate Gradient Solver with Inner Iterations and Variable-Step Preconditioning[END_REF] which is contradictory with the linear dependency of the vectors.

If we divide both terms in Equation 4 by α 0 and passing all terms to one side, we have

g 0 + p max -1 k=1 α k α 0 A k g 0 - 1 α 0 A p max g 0 = 0 ⇔ Ax 0 -b + p max -1 k=1 α k α 0 A k g 0 - 1 α 0 A p max g 0 = 0 ⇔ A         x 0 + p max -1 k=1 α k α 0 A k-1 g 0 - 1 α 0 A p max -1 g 0         = b ⇔ (6)
In practice to build this spaces all we have to do is to compute the basis of the space, but we will never use the "natural" base because it degenerates numerically as it grows. In practice if we use the regular double precision in a standard machine, after about 16 iteration the new values of the succession A p g 0 start to be linearly dependent, and depending on the matrix A some values are way too small or too big to be represented. for i = 1 to j do 7:

α i = (w, v i ) 8: w = w -α i v i 9:
end for 10:

v j+1 = w w 11: end loop solution can be written as

x p = x 0 + V p z p [START_REF] Axelsson | Generalized Conjugate Gradient Methods[END_REF] where z p is a p dimension vector. This approximation allows the writing of the error and residual vectors as e p := x px = x 0x + V p z p = e 0 + V p z p g p := Ax pb = Ae p = Ae 0 + AV p z p = g 0 + AV p z p [START_REF] Barka | Scalability of FETI-2LM methods on HPC clusters for antenna RCS applications[END_REF] A Krylov method consist, in one hand, of an algorithm to compute the base of the Krylov Space and on the other hand an optimal criteria to determine the approximate solution x p . This is made, a priori, by minimizing the error or residual using an adapted norm.

Lanczos method

We will now explain how to build the solution in the special case when A is a symmetric matrix. If we define h ij to be the coefficient of orthogonalization of Av j against v i , and also h j+1,j to be the norm of the w vector we get after the orthogonalization and V p the matrix of the first p vectors of the Arnoldi's base, then we have that

V t p V p = I p AV p = V p+1 H p+1,p (9)
where

H p+1,p :=                             
                             (10)
and so

H p := H pp = V t p AV p (11)
For the case of a symmetric matrix the matrix H p is also symmetric an thus tridiagonal, so the algorithm for the construction of the Base of Arnoldi is simplified. h j,j-1 = h j-1,j 7:

w = wh j-1,j v j-1 8:

h jj = (w • v j ) 9:
w = wh jj v j for symmetric matrices, now called Basis of Lanczos. This algorithm has the property of only use short recurrences in the computation so his cost is constant in every iteration.

If the matrix is also positive definite, the Lanczos method consist on minimize the error in the norm defined by A.

E(x p ) = x p -x 2 A = (A(x p -x) • (x p -x)) = (g p • e p ) (12)
This approximate solution has the next properties Theorem 0.5. The approximate solution x p of the Lanczos method is the projection of x in x 0 + K p for the inner product derived from A Proof. From 12 x p is the element from x 0 + K p which has a distance to x minimal in the A-norm.

Corollary 0.6. The residual vector g p = Ax pb of the Lanczos method is orthogonal to K p Proof. Direct of the properties of projections in affine spaces.

All it is missing is the practical computation of x p . To do so, of 12 and 8 we have that E(x p ) = (A(e 0 + V p z p ) • (e 0 + V p z p )) = (AV p z p • V p z p ) + 2(g 0 • V p z p ) + (g 0 • e 0 ) [START_REF] Chan | Collective Communication: Theory, Practice, and Experience: Research Articles[END_REF] To minimize this we only need to use the part that depends on z p and so the problem is reduce to the minimization of the functional

J p (z p ) = 1 2 (V t p AV p z p • z p ) + (V t p g 0 • z p ) = 1 2 (T p z p • z p ) + (y p • z p ) (14)
where T p = H p is the matrix of the orthonormalization coefficients. This is a classical finite dimension minimization problem from were we have the next results Lemma 0.7. If A is a symmetric positive definite matrix then J (x) = 1 2 (Ax

• x) -(b • x) is strictly convex Proof. =J (αx + (1 -α)y) = = 1 2 α 2 (Ax • x) + α(1 -α)(Ax • y) + 1 2 (1 -α) 2 (Ay • y) -α(b • x) -(1 -α)(b • y) = αJ (x) + (1 -α)J (y) + 1 2 (α 2 -α)(Ax • x) + 2α(1 -α)(Ax • y) + ((1 -α) 2 -(1 -α))(Ay • y) = αJ (x) + (1 -α)J (y) + 1 2 α(α -1) [(Ax • x) -2(Ax • y) + (Ay • y)] ( 15
)
Since A is positive definite, we have [START_REF] Ibrahima | Finite element approximations of the Lamé system with penalized ideal contact boundary conditions[END_REF] This functional has an absolute minimum in the unique point where his differential is zero, which is the point where Axb = 0

Corollary 0.9. The minimum of E(x p ) defined in 13 is the point x p = x 0 + V p z p , z p being the solution of the system T p z p = -y p [START_REF] Dolean | An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation: Other Titles in Applied Mathematics[END_REF] Proof. From previous theorem and lemma all we have to prove is that T p is positive definite which comes directly from the fact that the matrix A is also positive definite and the vectors V p are linearly independent.

The Lanczos method consist on building T p and y p then find z p and replace it in 7 to find the approach solution x p who minimizes the error in the A-norm.

One of the good things of this method is that the vectors of the base are calculated using a short recurrence, but the main drawback is that for the computation of x p we need to solve a bigger system every step, so the cost grows as the number of iteration.

Conjugate gradient

The Lanczos method will be a much better algorithm if it could use a short recurrence for the calculation of the approach solution x p . This can be done if the first component of z p are the ones of z p-1 which will give a formula of the type x p = x p-1 + α p v p [START_REF] Dongarra | A Proposal for a Set of Level 3 Basic Linear Algebra Subprograms[END_REF] In order to do so, the base of the Krylov space must be one in which the projection matrix W t p AW p is a diagonal one. But we will not be able to compute the error E(x p ) because e 0 will be unknown, the only way to test the convergence of the method is by using the dimensionless residual

Ax p -b b < ( 21
)
So we are forced to compute the successive gradients in order to control the method. Since g p ∈ K p+1 ∩ K ⊥ p we have that g p = ρv p+1 , so rather than using the orthonormal base of vectors v p we can use the orthogonal base of the gradients, even if it goes to zero, because we will stop before have any representation problem.

Let G p = (g 0 , g 1 , . . . , g p-1 ) as we just saw, G p = V p ∆ p where ∆ p is a diagonal matrix. The projection matrix is also tridiagonal symmetric positive definite G t p AG p = ∆ t p V t p AV p ∆ p = ∆ p T p ∆ p = Tp [START_REF] Dostál | FETI based algorithms for contact problems: scalability, large displacements and 3D Coulomb friction[END_REF] it admits the factorization Tp = Lp Dp Lt p and so

G t p AG p = Lp Dp Lt p ⇔ L-1 p G t p AG p L-t p = Dp ( 23
)
The previous shows that the matrix W p = G p L-t p made of linear combinations of G p is a A-orthogonal base of K p .

Since the projection matrix W t p AW p is diagonal, this base is ideal to use a short recurrence in the computation of the solution of the optimization problem 12, as it can be built using the relation

W p Lt p = G p (24)
Let (w 0 , w 1 , . . . , w p-1 ) be column vectors of W p and (γ 0 , γ 1 , . . . ) the sub-diagonal elements of L-t p , the equation 24 implies w 0 = g 0 and γ j-1 w j-1 + w j = g j ∀j > 0 [START_REF] Duff | Parallel implementation of multifrontal schemes[END_REF] With the different relations between x p , g p and w p we can formulate a new method without only using short recurrences in the construction of every new vector. Actually, from the previous equation we have g 0 = Ax 0b and w 0 = g 0 [START_REF] Dureisseix | A numerically scalable domain decomposition method for the solution of frictionless contact problems[END_REF] From the properties of the base W p we have

x p = x p-1 + ρ p-1 w p-1 ⇔ g p = g p-1 + ρ p-1 Aw p-1 [START_REF] Farhat | A scalable dualprimal domain decomposition method[END_REF] We know that g p is orthogonal to K p and so to w p-1 , this allows us to obtain ρ p-1 (g p ẇp-1 ) = (g p-1 • w p-1 ) + ρ p-1 (Aw p-1 • w p-1 ) = 0 ⇔ ρ p-1 = -(g p-1 • w p-1 )

(Aw p-1 • w p-1 ) (28)
From the equation 25 we can build the new w p with the previous w p-1 and g p w p = g pγ p-1 w p-1 [START_REF] Farhat | Optimal convergence properties of the FETI domain decomposition method[END_REF] The coefficient γ p-1 is also computed using the A-orthogonality relation between w p and w p-1

(w p • Aw p-1 ) = (g p • Aw p-1 ) -γ p-1 (w p-1 • Aw p-1 ) = 0 ⇔ γ p-1 = (g p • Aw p-1 ) (Aw p-1 • w p-1 ) (30)
The method defined this way i called the Conjugate gradient method and it can be resumed in 3. w p = g pγ p-1 w p-1 [START_REF] Chan | Collective Communication: Theory, Practice, and Experience: Research Articles[END_REF]: end loop From a theoretical point of view, we do not need any orthogonalization to build the new w p descent direction, but in limited arithmetic, errors are transmitted from each actualization, so in practice, to apply this method in the domain decomposition framework, we will need to reconjugate the vectors of the base, and so the storage of this vectors for a robust method is mandatory [START_REF] Roux | Acceleration of the Outer Conjugate Gradient by Reorthogonalization for a Domain Decomposition Method for Structural Analysis Problems[END_REF].

Considering this, the lines 11, 12 from the algorithm 3 are now replaced by the loop necessary to build

γ i = (g p • Aw i ) (Aw i • Aw i ) , i = 0, . . . , p -1 (31)
and

w p = g p - p-1 i=0 γ i w i ( 32
)
Proof. From the properties of projection in an affine space (A t A(x px) • w p ) = (A(x px) • Aw p ) = (g p • Aw p ) = 0, ∀w p ∈ K p [START_REF] Farhat | The two-level FETI method Part II: Extension to shell problems, parallel implementation and performance results[END_REF] We have naturally introduced the scalar product defined by A t A that is symmetric and positive definite if A is invertible. We could think that would be appropriate to use the Conjugate Gradient method to the system

A t Ax = A t b ( 35
)
This equation is call the "N ormal equation and it does not have a very good conditioning since it can be as much as the square of the original conditioning of A.

The best is to compute the solution using a short recurrence, so from the theorem 0.10 we know that the any A t A-orthogonal base of K p , W p implies that W t p A t AW p should be diagonal.

If we look at the structure of H p+1,p the matrix H t p+1,p H p+1,p is full, in this case the previous basis have no use. To build a A t A-orthogonal base we only need to apply the modified Gram-Schmidt procedure to the vectors obtained by successive multiplication for the matrix, using the A t A-norm. With this considerations we have the next short recurrences x p = x p-1 + ρ p w p g p = g p-1 + ρ p Aw p [START_REF] Farhat | Implicit Parallel Processing in Structural Mechanics[END_REF] and from the A t A-orthonormal properties we have

(g p • Aw p ) = 0 ⇔ (g p-1 • Aw p ) + ρ p (Aw p • Aw p ) = 0 ⇔ ρ p = -(g p-1 • Aw p ) (37)
and so the algorithm is described in 4

This algorithm forces to save the vectors of the base, but is numerically stable, Algorithm 4 ORTHODIR method 1: Initialization 2: g 0 = Ax 0b 3: w 0 = g 0 4: loop Iterate p = 1, ..., until convergence 5:

ρ = -

(g p-1 •Aw p-1 ) (Aw p-1 •Aw p-1 )
6:

x p = x p-1 + ρw p-1 7:

g p = g p-1 + ρAw p-1 8:

w p = Aw p-1 [START_REF] Bernardi | Domain Decomposition by the Mortar Element Method[END_REF]:

for i = 0 to p -1 do

10: γ = - (Aw p •Aw i ) (Aw i •Aw i ) 11:
w p = w p + γw i 12:

end for 13: end loop in any case we will deal with this issue in the next chapters.

A small variation of this method, equivalent to the regular and with the same properties can be made if we change the line 8 and build the next direction using the gradient, instead of the previous orthonormalized direction, i.e w p = g p [START_REF] Gosselet | On the initial estimate of interface forces in FETI methods[END_REF] of how to compute this products will be given when presenting the first domain decomposition method used in this work.

2. Global reduction operations to add the contribution of each process when computing the different scalar products.

Using the MPI libraries, or in general any communication protocol, needed to do the exchanges in parallel algorithm produces synchronization points in the code that need to be reduced as much as possible to avoid major impact in the total computation time.

With the previous basic linear algebra and parallel computing considerations we are ready to presents the main work of this thesis.

Chapter1

Hybrid FETI method

In the domain decomposition framework, the Finite Element Tearing and Interconnecting (FETI) methods have proven to be very effective in the solution of real engineering applications during the last years. This is one of the reason why its current development keeps up to this day, always searching for faster, precise and robust new versions. In this context we have developed a new method built from the base of two existing FETI methods, namely the FETI-1LM and FETI-2LM. This new method tries to recover the good properties of each method, in configurations where the use of one or the other is not so clear. For the development of this new method, first we need to understand the basics of both basic FETI methods from a theoretical point of view, but also the implementation which will be crucial to show how the new method works. This chapter is presented first introducing the FETI method in his classic version, including the preconditioners used to achieve good performance and some implementation considerations. Then the FETI method with two lagrange multipliers is explained to finally show the new hybrid method that arises from mixing the two previous FETI methods. The chapter ends with numerical results for all the methods.

Ω 2 Γ 3 Ω 1 Figure 1.1 -Two subdomain splitting 1.1 Basic FETI method 1.1.

Model problem and discretization

Model Problem

To begin, we will show the development of the method for a simple model, with the most basic domain decomposition configuration. All the ideas will later be extended to different elliptic problems and configurations. Let us first consider the Poisson problem with Dirichlet boundary condition

       -∆u = f in Ω u = 0 on ∂Ω (1.1)
where Ω ∩ R d , d = 2, 3 is a bounded domain. To find his variational form, the Stokes formula is used in 1.1, so the problem is now:

Find u ∈ H 1 0 (Ω) such that Ω ∇u • ∇v = Ω f v ∀v ∈ H 1 0 (Ω) (1.2)
The domain is now divided into two smaller subdomains Ω (1) and Ω (2) . Let Γ (3) = ∂Ω (1) ∩ ∂Ω (2) be the interface between both subdomains as in Figure 1.1. This division allows the formulation of two new smaller problems in each subdomain that inherits the Dirichlet condition, or any other condition, in a part of the boundary. This problems are written as

       -∆u (s) = f (s) in Ω (s) u (s) = 0 on ∂Ω (s) \ Γ (3) (1.3) with s = 1, 2.
It is clear that the solution of 1.1 will satisfy this equations, but the contrary is not always true, due to the differences that occur in Γ (3) .

We continue by using again the Stokes formula to find the variational form of 1.3

Ω ∇u (s) ∇v (s) = Ω f (s) v (s) + Γ (3) ∂u (s) ∂n (s) v (s) , ∀v (s) ∈ H 1 0∂Ω (s) \Γ (3) (Ω (s) ) (1.4)
For a function v ∈ H 1 0 (Ω), in particular the solution of the problem, his restriction in the subdomains Ω (s) will be continuous on the interface Γ (3) . On the other hand two functions u (s) that satisfies the local Laplace equations 1.3 will not necessarily share the same values on Γ (3) . Instead they can be used to build a more general global solution, but not smooth enough as required.

To do this construction, but at the same time recover the unique solution of 1.1 the two variational equations 1.4 are added, giving the following variational equality

Ω ∇u∇v = Ω f v + Γ (3)
∂u (1) ∂n (1) + ∂u (2) ∂n (2) v (3) ∀v ∈ H 1 0 (Ω) (1.5) where v (3) = v (1)

| Γ (3) = v (2) | Γ (3) .
This formulation shows that a new condition is necessary to have an equivalence between the solution of the global problem and the solution of the local ones. For u to be in H 1 0 (Ω) this admissibility condition imposes the continuity on the interface u (1) = u (2) on Γ (3) (1.6)

and also for the same reason (explicit in equation 1.5) we need a condition on the flux ∂u (1) ∂n (1) + ∂u (2) ∂n (2) = 0 on Γ (3) (1.7)

In general, a non overlapping domain decomposition method consist in introducing boundary conditions on Γ (3) to complement the local equations 1.3 and to iteratively find the values of these boundary conditions for which both continuity 1.6 and equilibrium 1.7 interface conditions are satisfied, meaning that this solution will be the exact same as the global searched one.

Depending on the condition imposed, two basic methods can be derived, the Schur complement method and the FETI method, later we will show a third method, also of the FETI type, that comes from using a different condition on the interface.

The Schur complement method consists in enforcing consistent Dirichlet boundary conditions on Γ (3) so continuity condition 1.6 is automatically satisfied u (1) = u (2) = u 3 on Γ (3) (1.8)

The local Dirichlet problem to be solved in parallel for a given u 3 on each subdomain is

           -∆u (s) = f (s) in Ω (s) u (s) = 0 on ∂Ω (s) \ Γ (3)
u (s) = u 3 on Γ (3) (1.9)

We reduce the computation to find the value of u 3 for which the equilibrium interface condition 1.7 is satisfied. From equations 1.9, the functions ∂u (1) ∂n (1) and ∂u (2) ∂n (2) are continuous depending on u 3 . The Schur complement method consists in solving iteratively a condensed interface problem to find u 3 whose residual is equal to ∂u (1) ∂n (1) + ∂u (2) ∂n (2) . The FETI method is based in enforcing consistent Neumann boundary conditions on Γ (3) so now the equilibrium interface condition 1.7 is automatically satisfied:

∂u (1) ∂n (1) = -∂u (2) ∂n (2) = λ on Γ (3) (1.10) the local Neumann problem to be solved in parallel for a given λ on each

subdomain is              -∆u (s) = f (s) in Ω (s) u (s) = 0 on ∂Ω (s) \ Γ (3) ∂u (s) ∂n (s) = ±λ on Γ (3) (1.11)
Now we compute the value λ on the interface for which the continuity condition 1.6 is satisfied. From equations 1.11, u (1) | Γ (3) and u (2) | Γ (3) are continuous functions depending on λ. The FETI method consists in solving iteratively a condensed interface problem to find λ and whose residual is equal to u (1)

| Γ (3) -u (2) | Γ (3) .
A different interpretation of the FETI method can be considered if we see the unknown λ as the Lagrange multiplier of the continuity condition 1.6. The solution of the global variational problem 1.2 is the field u of H 1 0 (Ω) that minimizes the energy functional

J(v) = 1 2 Ω ∇v • ∇v - Ω f v (1.12)
This minimization problem is equivalent to finding the couple of fields (u (1) , u (2) ) of (2) ) that minimizes the sum of the local energy functionals

H 1 0∂Ω (1) \Γ (3) (Ω (1) ) × H 1 0∂Ω (2) \Γ (3) (Ω
J 1 (v (1) ) + J 2 (v (2) ) = 1 2 Ω (1)
∇v (1) • ∇v (1) -

Ω (2) f (2) v (2) + 1 2 Ω (2)
∇v (2) • ∇v (2) -

Ω (2)
f (2) v (2) (1.13) under the continuity constraint u (1) | Γ (3) = u (2) | Γ (3) . This condition can be written under the weak form

Γ (3) u (1) -u (2) µ = 0 ∀µ ∈ H -1 2 (Γ (3) ) (1.14)
Now, consider the Lagrangian L(v (1) , v (2) , µ) = 1 2 Ω (1)

∇v (1) • ∇v (1) -

Ω (1) f (1) v (1) + 1 2 Ω (2)
∇v (2) • ∇v (2) -

Ω (2) f (2) v (2) - Γ (3) v (1) -v (2) µ (1.15) Γ Ω Ω 1 2 3 Figure 1.
2 -Two subdomain divisions with duplicated nodes.

and we note that the saddle point (u (1) , u (2) , λ) of L in H 1 0∂Ω (1) \Γ (3) (Ω (1) )×H 1 0∂Ω (2) \Γ (3) (Ω (2) )× H -1 2 (Γ (3) ) is precisely the point where the variational equations 1.11 and 1.14 are satisfied.

Discretization

Lets consider a discretization of variational equation 1.2 using a finite element method. This process works for different elliptic partial differential equations and different finite element discretizations. So from now on we can consider this as a more general work, as long as we have this type of discretization that will lead to a system of the following form

Kx = f (1.16)
The global stiffness matrix of the discrete problem can be arranged to have the block structure showed in equation 1.17, where subscripts i denote the inner degrees of freedom of subdomains Ω (1) and Ω (2) and subscript b is used for the nodes on the interface Γ (3) = ∂Ω (1) ∩ Ω (2)

            K (1) ii 0 K (1) ib 0 K (2) ii K (2) ib K (1) bi K (2) bi K bb                         x (1) i x (2) i x b             =             f (1) i f (2) i f b             (1.17)
The formulation of each local discretization matrix is made considering that each subdomain has its own mesh and also that the nodes of the interface Γ (3) are shared by both meshes as in Figure 1.2. So there are two interface blocks, one in each local matrix, noted with superscripts (1) and (2) . The local stiffness matrices of the two subdomains are

K (1) =        K (1) ii K (1) ib K (1) bi K (1) bb        K (2) =        K (2) ii K (2) ib K (2) bi K (2) bb        (1.18)
where

K (1) bb + K (2) bb = K bb .
The discretization of variational formulation of equation 1.4 in subdomain Ω (s) leads to the following systems of equations

       K (s) ii K (s) ib K (s) bi K (s) bb               x (s) i x (s) b        =        f (s) i f (s) b + h (s) b        (1.19) where f (1) b + f (2) b = f b and h (s)
b is the vector representing the discretization of the flux ∂x (s) ∂n (s) on Γ (3) .

From this we have and explicit relation between the inner an the interface nodes

x (s) i = K (s) -1 ii f (s) i -K (s) -1 ii K (s) ib x (s) b (1.20)
from 1.20 and 1.19 the relation linking the trace and the flux of a vector satisfying the inner subset is derived

h (s) b = K (s) bi x (s) i + K (s) bb x (s) b -f (s) b = K (s) bi (K (s) -1 ii f (s) i -K (s) -1 ii K (s) ib x (s) b ) + K (s) bb x (s) b -f (s) b = (K (s) bb -K (s) bi K (s) -1 ii K ib )x (s) b -(f (s) b -K (s) bi K (s) -1 ii f (s) i ) = S (s) bb x (s) b -c (s) b (1.21) S (s) bb = K (s) bb -K (s) bi K (s) -1 ii K (s)
ib is the Schur complement matrix. It is the discretization of the Dirichlet to Neumann mapping that defines the bi-continuous one to one correspondence between the trace and the flux on the boundary (or interface in our case) of a field that satisfies the Laplace equation inside the subdomain. It is symmetric positive definite if the K (s) matrix is symmetric positive definite.

The discretization of continuity 1.6 and flux 1.7 are

x (1) b = x (2) b (1.22) h (1) b + h (2) b = 0 (1.23)
The last condition, also called equilibrium, combined with 1.19, gives the following interface equation

K (1) bi x (1) i + K (1) bb x b -f (1) b + K (2) bi x (2) i + K (2) bb x b -f (2) b = 0 ⇔ K (1) bi x (1) i + K (2) bi x (2) i + (K (1) bb + K (2) bb )x b = f (1) b + f (2) b (1.24)
Finally, for two vectors defined on subdomains Ω (1) and Ω (2) to be considered as the restrictions of the solution of the global discrete problem 1.17, they must meet

• the inner equations in each subdomain        K (1) ii x (1) i + K (1) ib x (1) b = f (1) i K (2) ii x (2) i + K (2) ib x (2) b = f (2) i (1.25)
• the interface equation

K (1) bi x (1) i + K (2) bi x (2) i + K (1) bb x (1) b + K (2) bb x (2) b = f (1) b + f (2) b (1.26)
• the continuity across the interface Γ (3) x are both equal to the restriction of the global solution on Γ (3) then the inner equations 1.25 are the first two rows of 1.17 and the interface equations 1.26 are the third row. Meaning that the methodology derived only from linear algebra is valid for any finite element discretization of PDEs.

(1) b = x
The inner equations 1.25 are common solution vectors of local problems for any kind of boundary conditions on Γ (3) . Equations 1.27 and 1.26 are the actual condensed interface problem since the inner equations 1.25 establish that x (1) i and x (2) i can be derived from x

FETI one lagrange multiplier

The previous ideas are now formalized for the extended case where we have N s > 2 subdomains and and interface Γ , defined as

Γ = 1≤s,q≤N s (∂Ω (s) ∩ ∂Ω (q) )
(1.28)

In the FETI method the discrete flux, noted λ, is the unknown defined in the nodes along the interface and the jump of the solutions of the local Neumann problems are the gradient of the condensed interface problem. The discretization of the local Neumann problems in some subdomain s, can also be written as

K (s) x (s) = f (s) + t (s) T B (s) T λ (1.29)
where K (s) is the local stiffness matrix, f (s) the right-hand side vector, t (s) ∈ M #(∂Ω (s) )×#(Ω (s) ) are trace operators which extracts boundary degrees of freedom from subdomain Ω (s) and B (s) ∈ M #(Γ )×#(∂Ω (s) ) are discrete assembling matrices which connects pairwise degrees of freedom on the interface. In a general case, the stiffness matrix K (s) , the local solution vector x (s) and the local right-hand side vector f (s) are defined as

K (s) =        K (s) ii K (s) ib K (s) bi K (s) bb       
(1.30)

x (s) =        x (s) i x (s) b        , f (s) =        f (s) i f (s) b        (1.31) Figure 1.3 -Multiple interface node.
The i and b subscripts also denotes the interior and interface nodes of the subdomain respectively. So trace operator applied to the solution x (s) is such that

x (s) b = t (s) x (s) (1.32)
The discrete operator B (s) is a defined as the mapping of a vector in local interface ∂Ω (s) on the complete interface, so applied to the solution on the local interface x

(s)

b we can write the continuity condition across the total interface, degree of freedom per degree of freedom s B (s) t (s) x (s) = 0 (1.33) the restriction of B (s) on interface Γ ij , noted B (ij) , is defined as a signed boolean operator such that B (ij) and B (ji) have opposite signs, providing the continuity needed. Any node that it's in more than two subdomains, Figure 1.3, will generate the same number of continuity conditions and flux as the number of interfaces who shares it.

With the definition of B (s) and t (s) , the solutions x (s) of 1.29 and 1.33 are the searched restrictions in every subdomain of the global discrete solution of 1.17. The vectors x (s) are actually continuous and the definition of B (s) is that t (s) T B (s) T λ is zero for the inner nodes of Ω (s) and in the interface we have B (ij) T λ+B (ji) T λ = 0 thanks to the opposite sign. So, again, the assembly of local discrete Neumann equations 1.29 gives exactly global discrete equation 1.16.

The gradient of the condensed interface problem is defined as

g = s B (s) t (s) x (s) (1.

"Floating" subdomains

For most subdomains, we face the common case of finding that ∂Ω (s) ∩ ∂Ω is void, this means that the Dirichlet condition of the problem is not in Ω (s) , so the local discrete Neumann equations 1.29 are ill posed. If K (s) comes from the Laplace equation, its kernel are the constant fields in the subdomain, if it comes from three-dimensional linear elasticity, then the kernel is the subspace of rigid body motions, of dimension 6 in the case of subdomains simply connected.

The pseudo-inverse K (s) + is now needed and the Cholesky factorization with partial pivoting is used in the matrix K (s) to achieve this, and also because it allows to compute a generator of the kernel R (s) and a factorization of a maximal full rank sub-block. Given the pseudo-inverse K (s) + and the kernel generator R (s) , the solution x (s) of the discrete system of equation 1.29 can be written as a particular solution plus an undefined element of the kernel of K (s)

x (s) = K (s) + (f (s) + t (s) T B (s) T λ) + R (s) α (s) (1.35)
From equation 1.29 we see that the right-hand side belongs to the range space of matrix K (s) , and so is orthogonal to the kernel. This orthogonality constraint can be written

R (s) T f (s) + t (s) T B (s) T λ = 0 (1.36)
This last equation is the admissibility condition for the forces of a floating subdomain. Its interpretation is that fields belonging to the kernel must have zero energy.

Condensed interface problem

Replacing x (s) from equation 1.35 in the continuity condition 1.33 we have

s B (s) t (s) K (s) + t (s) T B (s) T λ + s B (s) t (s) R (s) α (s) = - s B (s) t (s) K (s) + f (s) (1.37)
To build the condensed interface problem this previous equation is used, and the equation 1.36, leading to the problem to be satisfied by λ and the vector of coefficients of the kernel components α

      F G G T 0             λ α       =       d c       (1.38)
Where: s) , jump of zero energy fields defined by α (s) in Ω (s)

• F = s B (s) t (s) K (s) + t (s) T B (s) T = s B (s) S (s) + bb B (s) T dual Schur complement matrix • Gα = s B (s) t (s) R (s) α (
• G T λ = . . . , B (s) t (s) R (s) , . . . T λ • d = - s B (s) t (s) K (s) + f (s)
• c = . . . , -b (s) T R (s) , . . .

T

The condensed interface system 1.38 is an hybrid system. It's solution λ satisfies the following orthogonality condition

µ T Fλ = µ T d, ∀µ/G T µ = 0 (1.39)
Now consider λ 0 , for example

λ 0 = AG(G T AG) -1 c (1.40)
where A is a matrix symmetric positive definite that is usually taken as being the identity, but it can also be defined as the preconditioner, later to be defined, or some scaling matrix. For details, see [START_REF] Daniel | Theoretical comparison of the FETI and algebraically partitioned FETI methods, and performance comparisons with a direct sparse solver[END_REF].

Then this λ 0 satisfies the admissibility constraint of equation 1.36 and G T (λλ 0 ) = 0. So, if P is any projector in the kernel of G T , then from equation 1.39 we have that λ is the solution of the following projected problem

P T FP (λ -λ 0 ) = P T (d -Fλ 0 ) (1.41)
The FETI method solves iteratively via a conjugate gradient algorithm the previous projected condensed interface problem 1.41, using the orthogonal projector in the kernel of G T .

Interpretation of projector P

The orthogonal projection in the kernel of G T can be written algebraically

P = I -AG G T AG -1 G T (1.42)
To compute the projection of a given vector g we mainly solve the problem

G T AG α = -G T g (1.43)
which is a global coarse grid problem whose unknowns are the coefficients of zero energy components of the solutions in the floating subdomains.

Now for for a given approximation λ p of the flux on the interface, the residual of the condensed interface problem is

g p = s B (s) S (s) + bb B (s) T λ p + s B (s) t (s) K (s) + f (s) = s B (s) t (s) x (s) p+ (1.44)
where u p+ i is the solution of the local Neumann problems, computed using the pseudo-inverse matrices

x (s) p+ = K (s) + f (s) + t (s) T B (s) T λ p (1.45)
so the gradient is equal to the jump of this particular solutions. From equations 1.42 and 1.43, the projected gradient P g p is

P g p = g p + AGα p = s B (s) t (s) x (s) p+ + s B (s) t (s) R (s) α (s) p (1.46)
So the projected gradient P g p is equal to the jump of the local particular solutions of Neumann problems x (s) p+ plus the term of zero energy fields with coefficients

α (s) p x (s) p = x (s) p+ + R (s) α (s) p (1.47)
The definition of the constraint 1.36 associated with the orthogonal projector 1.43 entails that the zero energy components of the jump of the local solution fields x (s) p are minimal in the sense that this jump is orthogonal to all the traces of zero energy fields

G T P g p = 0 ⇔ B (s) t (s) R (s) T P g p = 0 ∀s (1.48)
Computing the projected gradient P g p consists in fact in computing the coefficients α p of optimal local zero energy fields. For the linear elasticity problem, the zero energy fields are the rigid body motions. The underlying process is a kind of coarse grid smoothing of approximate solution that ensures a convergence rate for the overall FETI process asymptotically independent upon the number of subdomains. Hence, the FETI method with floating subdomains is a kind of two-level solver that is numerically scalable.

Local preconditioner

The coarse grid smoothing performed by the zero energy fields projector P gives a convergence rate independent upon the number of subdomains, but this is not enough to have a convergence rate that is also independent upon the mesh size. It's mandatory the use of a preconditioner, which for the FETI method is one of the "Dirichlet" type.

Consider t (s) the trace or restriction operator on the local interface of subdomain Ω (s) . Then the contribution of subdomain s to the condensed interface operator, defined in1.37, is

B (s) t (s) K (s) + t (s) T B (s) T (1.49)
and it just depends on the restriction in the interface Γ (s) = ∂Ω (s) of the pseudoinverse of K (s) , meaning

t (s) K (s) + t (s) T = (K (s) bb -K (s) bi K (s) -1 ii K (s) ib ) + =: S (s) bb + (1.50)
As this matrix is the pseudo-inverse of the Schur complement matrix on the interface Γ (s) , a preconditioner based on local contributions for FETI is

D -1 = s B (s) S (s) bb B (s) T (1.51)
where again S (s) bb is the Schur complement. This preconditioner tries to approximate the global inverse of local sums by the sum of the local inverse, meaning

      s B (s) S (s) + bb B (s) T       + s B (s) S (s) bb B (s) T (1.52)
The computation of this preconditioner applied to an interface vector w is done by solving the following local problem with Dirichlet boundary conditions on the interface Γ (s) defined by the assembled local vector t (s) T B (s) T w

      K (s) ii K (s) ib 0 I              w(s) i w(s) b        = t (s) T B (s) T w =       0 w (s) b       (1.53) whose solution is w(s) i = -K (s) -1 ii K (s) ib w (s) b (1.54)
and multiplying by the stiffness matrix to obtain

       K (s) ii K (s) ib K (s) bi K (s) bb               w(s) i w(s) b        =        0 (K (s) bb -K (s) bi K (s) -1 ii K (s) ib ) w(s) b        =       0 S (s) bb w (s) b       (1.55)
With this preconditioner it has been proved [START_REF] Mandel | Convergence of a substructuring methodwith Lagrange multipliers[END_REF], [START_REF] Farhat | Optimal convergence properties of the FETI domain decomposition method[END_REF], that the FETI method is asymptotically independent upon the mesh size. With a condition number for the projected condensed operator bounded by

C(1 + log( H h )) 2 (1.56)
where h is the mesh size and H is the characteristic subdomain size. Meaning, that when decreasing H the number of subdomains increases and when decreasing h the finite element mesh is refined (more elements).

So, with both local and global preconditioners, the number of iterations does not depend anymore on both the number of subdomain nor the mesh size. A second preconditioner, not mathematically optimal can be introduced, as an alternative to the cost of implementation of the Dirichlet one, the so called "Lumped" preconditioner that can de defined as

L -1 = s B (s) K (s) bb B (s) T (1.57)
where

K (s)
bb is the finite element discretization matrix on the interface nodes. This preconditioner works as an approximation of the local Schur complements, with a much more economical implementation because it does not require any additional storage and involves only matrix-vector products of sizes equal to the subdomain interfaces.

Both preconditioner were generalized to treat heterogeneous problems [START_REF] Daniel | A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems[END_REF] by just redefining the Boolean operator B (s) to a more general one B(s) := β (s) B (s) (1.58) such that s B(s) B(s) T = I, with β (s) a diagonal scaling matrix, usually based on the diagonal coefficients of the local stiffness matrix on the interface, the so called super-lumped scaling. This scaling is a mechanical consistent combination of the interface reaction forces from the Dirichlet problem in each subdomain. With this new scaled assembling operators, the preconditioners are written as

D -1 = s B(s) S (s) bb B(s) T (1.59) L -1 = s B(s) K (s) bb B(s) T (1.60)
The behaviour of both preconditioner can be seen [START_REF] Farhat | Optimal convergence properties of the FETI domain decomposition method[END_REF], [START_REF] Daniel | A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems[END_REF] and will be evocated in later sections. w p+1 = P D -1 P T g p+1 10:

for i = 0 to p do 11:

γ i = - (w i , Fw p+1 ) (w i , Fw i ) 12:
w p+1 = w p+1 + γ i w i 13:

Fw p+1 = Fw p+1 + γ i Fw i 14:
end for 15: end loop As we can see, a full reconjugation, instead of the classical CG update without any, is also used (see the line 10 of the algorithm), because in limited arithmetic the orthogonal properties of the CG method are lost, specially in the context of the FETI methods, where the multiplication by the operator is not totally accurate, making this part also crucial for good convergence rate [START_REF] Roux | Acceleration of the Outer Conjugate Gradient by Reorthogonalization for a Domain Decomposition Method for Structural Analysis Problems[END_REF]. Fw p+1 = P D -1 P Fw p+1 13: for i = 0 to p do 14:

Algorithm 6 FETI-1LM unsymmetric 1: Initialization 2: λ 0 = AG[G T AG] -1 (-R T c
γ i = - (Fw i ) T (Fw p+1 ) (Fw i ) T (Fw i ) 15:
w p+1 = w p+1 + γ i w i

16:

Fw p+1 = Fw p+1 + γ i Fw i 17:

end for 18:

end loop 19: end loop

Unsymmetric FETI algorithm

The previous method and algorithm only works when the matrix F is symmetric positive definite, for a more general case when our FETI operator F is no longer symmetric or non positive definite, as we will see in later sections, the CG algorithm is no longer appropriate, we use instead the ORTHODIR method with left preconditioner for unsymmetric matrices, mainly for its implementation, simplicity and good properties, as is equivalent to the GMRES algorithm [START_REF] Magoules | Parallel Scientific Computing[END_REF]Chapter 12]. One of the main theoretical differences between them two is the storage of previous computed directions which, in any case, it is done to perform a full reorthogonalization that allows the loss of orthogonality of the search directions and keeps a good convergence ratio.

The full description is done in Algorithm 6.

FETI with two lagrange multipliers 1.2.1 FETI-2LM

We now present the second FETI method that is used as a base to build our method. Originally introduced in [START_REF] Farhat | Two-level domain decomposition methods with Lagrange multipliers for the fast iterative solution of acoustic scattering problems[END_REF] as a solver for acoustic problems then extended in [START_REF] Roux | Optimization of Interface Operator Based on Algebraic Approach[END_REF] and [START_REF] Roux | Approximation of Optimal Interface Boundary Conditions for Two-Lagrange Multiplier FETI Method[END_REF] as a robust solver for more general problems. The basic idea of this method is to impose Robin boundary conditions in the interface to "glue" the local solutions in order for them to be the restrictions of the global solution. To have a method useful for different approximations of elliptic PDE, it will as in FETI, be derived from linear algebra.

To understand how this method works we start with a simple partition of the total domain Ω into 2 subdomains, as in Figure 1.1. From this splitting and some linear elliptic PDE problem defined in global domain Ω we use any finite elements discretization coming from this equations. Again, we use subscript i for internal nodes and b for the interface ones, the interface defined as Γ (3) = ∂Ω (1) ∩ Ω (2) . Lets consider the contribution of each subdomain Ω (s) , s = 1, 2 to the matrix and right-hand side the finite element discretization

K (s) =        K (s) ii K (s) ib K (s) bi K (s) bb        , f (s) =        f (s) i f (s) b        (1.61) where K (s)
bb represents the interaction matrices between the nodes on the interface and those from the interior of Ω (s) independent of each other, the same happens for f

            K (1) ii 0 K (1) ib 0 K (2) ii K (2) ib K (1) bi K (2) bi K bb                         x (1) i x (2) i x b             =             f (1) i f (2) i f b             (1.62) with K bb = K (1) bb + K (2) bb and f b = f (1) b + f (2)
b . Taking the second line to find x (2) i and replacing it in the first and third lines, the system that remains can be written in the following way

       K (1) ii K (1) ib K (1) bi K (1) bb + (K (2) bb -K (2) bi K (2) ii -1 K (2) ib )              x (1) i x b       =        f (1) i f (1) b + (f (2) b -K (2) bi K (2) ii -1 f (2) i )        (1.63)
The same treatment can be done, but eliminating the unknowns of subdomain 1. And so, for both subdomains we have

       K (s) ii K (s) ib K (s) bi K (s) bb + S (q) bb              x (s) i x b       =        f (s) i f (s) b + c (q) b        (1.64) with S (q) bb = K (q) bb -K (q) bi K (q) ii -1 K (q)
ib the Schur complement and c

(q) b := f (q) b -K (q) bi K (q) ii -1 f (q) i
the condensed right hand side in the subdomain q, opposite to s.

Previous equation means that the restriction of the global solution is, in each subdomain, a local solution of a problem with generalized Robin boundary condition in the interface. The operator of the generalized Robin condition is the Dirichlet to Neumann operator of the rest of the domain. In fact, the Dirichlet to Neumann operator describes exactly the behaviour of the boundary of the outer domain. Enforcing a local generalized Robin boundary condition using the Dirichlet to Neumann operator of the rest of the domain makes the interface to behave locally exactly as the rest of the domain forces it to do. Then, the local problem formulation takes exactly into account the coupling between the subdomain and the rest of the domain.

Local generalized Robin boundary conditions enforced on both sides of the same interfaced should be set up independently. This leads to introduce two independent interface variables, one for each side.

       K (s) ii K (s) ib K (s) bi K (s) bb + A (s) b               x (s) i x (s) b        =        f (s) i f (s) b + λ (s)        (1.65)
from this local equations, inner equations can be eliminated to get the equivalent local condensed problem

(K (s) bb -K (s) bi K (s) ii -1 K (s) ib + A (s) b ) + x (s) b = f (s) b -K (s) bi K (s) ii -1 f (s) i + λ (s) (S (s) bb + A (s) b )x (s) b = c (s) b + λ (s) (1.66)
And the explicit relation between the trace of the solution and the Lagrange multiplier is

x (s) b = (S (s) bb + A (s) b ) -1 λ (s) + (S (s) bb + A (s) b ) -1 c (s) b (1.67)
Since the boundary conditions on both sides are not consistent, the two interface conditions 1.27 and 1.26 must be enforced to make the solutions of the local problems 1.65 to be the restrictions of the solution of the global problem 1.17. Thanks to the second line of local equation 1.65, the interface equilibrium condition gives

A (1) b x (1) b + A (2) b x
(2) b = λ (1) + λ (2) (1.68)

The residual of the condensed interface problem is then

x (1) b -x (2) b = 0 λ (1) + λ (2) -A (1) b x (1) b -A (2) b x (2) b = 0 (1.69) where x (s)
b is a function of λ (s) according to equation (1.67). This two equations, combined with the local problem 1.65 are equivalent to the global block problem 1.62. In this form, equations 1.69 have a mixed form, to avoid this, a new combination of both equations is done, leading to a new condensed interface problem λ (1) + λ (2) b from 1.67 and using it in (1.70), the condensed interface problem used to compute λ (1) and λ (2) , can be explicitly defined

-(A (1) b + A (2) b )x (2) b = 0 λ (2) + λ (1) -(A (2) b + A (1) b )x
       I I -(A (1) b + A (2) b )(S (2) bb + A (2) b ) -1 I -(A (2) b + A (1) b )(S (1) bb + A (1) b ) -1 I              λ (1) λ (2)       =        (A (1) b + A (2) b )(S (2) bb + A (2) b ) -1 c (2) b (A (2) b + A (1) b )(S (1) bb + A (1) b ) -1 c (1) b        (1.71)
This is method is called FETI-2LM, like 2-Lagrange multiplier FETI method.

Arbitrary mesh partition

Lets show the extension of the previous method to the case of an arbitrary mesh with a total number of subdomains N s > 2. The idea is to take the 2 subdomain case and apply it to the different interfaces created between neighbours subdomains.

To define this we start by dividing the total interface boundary Γ into interface edges Γ j considering that an interface edge is a collection of connecting interface nodes. The total interface can be written as and we can define more precisely Γ j as

Γ = 1≤s,q≤N s (∂Ω (s) ∩ ∂Ω (q) ) (1.72) Ω 1 Ω 3 Ω 2 Ω 4 Γ 2 Γ 1 Γ 3 Γ 4 Γ 5 Γ 6
Γ j := Γ (sq) = ∂Ω (s) ∩ ∂Ω (q) , ∀s, q : ∂Ω (s) ∩ ∂Ω (q) ∅ (1.73)
The crosspoint are nodes where more than two subdomains are connected, see Figure 1.4 for an example of 4 subdomains and one crosspoint. We can see here, that two edges were created in the crosspoint, as diagonal subdomains are neighbours. For each interface edge Γ j and for each subdomain Ω (s) that intersects it, lets define a boolean matrix B (s)

Γ j ∈ M dim(Γ j )×dim(Ω (s) ) (R) by B (s) Γ j v (s) = v (s) b | Γ j (1.74)
where the vector v (s) defined in all the nodes of subdomain s can be written by separating the nodes of the complete local interface from the interior ones as

v (s) =        v (s) i v (s) b        .
Remark: This definition of interface edges can produce redundancies in the formulation of the problem. In practice this redundancies do not produce any negative impact in the FETI methods, on the contrary, for some cases such as the preconditioning in FETI-1LM they play a vital role, as we will see later in the text.

With previous definitions we can rewrite the local problem 1.65 and the continuity conditions 1.69 for the general case

        K (s) + Γ j ⊆Ω (s) B (s) T Γ j A (s) Γ j B (s) Γ j         x (s) = f (s) - Γ l ⊆Ω (s) B (s) T Γ l λ (s) Γ l , s = 1, . . . , N s B (s) Γ j x (s) -x b | Γ j = 0 ∀Γ j , Ω (s) ⊇ Γ j Ω (q) ⊇Γ j λ (q) Γ j -A (q) Γ j x b | Γ j = 0 ∀Γ j (1.75)
where, A (s)

Γ j is the augmentation matrix in subdomain Ω (s) for the interface Γ j , the same way λ (s) Γ j is the lagrange multipliers in the interface edge Γ j on the side of the subdomain s and the vector x b | Γ j is the restriction of the global unknown x to Γ j . Lets define also, for every Ω (s) the total augmentation matrix

A (s) := Γ j ⊆Ω (s) B (s) T Γ j A (s) Γ j B (s) Γ j (1.76)
Finally we can do the condensation of previous equations to obtain the interface problem

Ω (q) ⊇Γ j         λ (q) Γ j -A (q) Γ j B (s) Γ j K (s) + A (s) -1 Γ l ⊆Ω (s) B (s) T Γ l λ (s) Γ l         = - Ω (q) ⊇Γ j A (q) Γ j B (s) Γ j K (s) + A (s) -1 f (s) , ∀Γ j , ∀Ω (s) ⊇ Γ j (1.77)
This problem can also be written as Fλ = d where F is, in general, a non symmetric matrix, defined in the whole interface, so again, the ORTHODIR method with full reorthogonalization is used to solve the problem. One of the differences with FETI-1LM is that there are no efficient local preconditioners for this method, so for the cases were convergence is expected in both methods, the 2LM may be slower. On the contrary, there are cases (e.g. anisotropic materials, heterogeneities in the interface, etc.) when bad numerical features located across the interface will not allow the convergence of FETI-1LM, but the use of two independent lagrange multipliers grant a good handle of this problems and convergence is achieved, making the FETI-2LM a more robust method.

Ω S-1 Γ b-1 Ω S Ω S+1 Γ b

Optimal Interface Boundary Conditions

As seen in the definition of the FETI-2LM method, a big part of it, comes from the definition of the augmentation matrix A (s) . Different work as been done to achieve the best numerical choice, next we will show some of them.

From a theoretical point of view, in the case of 2 subdomain the optimal augmentation matrices are defined as bb , s = 1, 2, q = 1, 2, s q, the matrix that defines the method 1.71 is equal to the identity.

In a more general case of a one way division in N subdomains with no crosspoints, the optimal augmentation can be proved to be the Schur complement of the rest of the domain. Consider a one-way splitting as in Figure 1.5 and denote the nodes in the interface with the subscript b -1 for the left interface and b for the right. Then the local contributions can be written as

K (s) =             K (s) ii K (s) ib-1 K (s) ib K (s) b-1i K (s) b-1b-1 0 K (s) bi 0 K (s) bb             , f (s) =             f (s) i f (s) b-1 f (s) b             (1.79)
If we elimiinate the inner nodes, the contributions to the condensed matrix and right hand side of subdomain Ω (s) is

       S (s) b-1b-1 S (s) b-1b S (s) bb-1 S (s) bb        =         K (s) b-1b-1 -K (s) b-1i K (s) ii -1 K (s) ib-1 -K (s) b-1i K (s) ii -1 K (s) ib -K (s) bi K (s) ii -1 K (s) ib-1 K (s) bb -K (s) bi K (s) ii -1 K (s) ib         (1.80)        c (s) b-1 c (s) b        =         f (s) b-1 -K (s) b-1i K (s) ii -1 f (s) i f (s) b -K (s) bi K (s) ii -1 f (s) i         (1.81)
After assembling the local contributions, the global condensed problem is a block tridiagonal system

                 . . . . . . 0 0 S (s-1) b-1b-2 S (s-1) b-1b-1 + S (s) b-1b-1 S (s) b-1b 0 0 S (s) bb-1 S (s) bb + S (s+1) bb S (s+1) bb+1 0 0 . . . . . .                  ,                  . . . c (s-1) b-1 + c (s) b-1 c (s) b + c (s+1) b . . .                  (1.82)
If this system is factorized by successive condensations from both sides up to the subdomain Ω (s) , the following condensed problem is obtained in subdomain

s        S (-) b-1b-1 + S (s) b-1b-1 S (s) b-1b S (s) bb-1 S (s) bb + S (+) bb               x (s) b-1 x (s) b        =        c (-) b-1 + c (s) b-1 c (s) b + c (+) b        (1.83)
where the terms with (+) and (-) superscript are recursively denoted by

S (-) b-1b-1 = S (s) b-1b-1 -S (s) b-1b S (s) bb + S (+) bb -1 S (s) bb-1 S (+) bb = S (s) bb -S (s) bb-1 S (-) b-1b-1 + S (s) b-1b-1 -1 S (s) b-1b c (-) b-1 = c (s) b -S (s) bb-1 S (-) b-1b-1 + S (s) b-1b-1 -1 c (-) b-1 + c (s) b-1 c (+) b = c (s) b-1 -S (s) b-1b S (s) bb + S (+) bb -1 c (s) b + c (+) b (1.84)
Previous system suggest that the optimal term to add to the local matrix problem K (s) is S (-) on the left interface nodes and S (+) on the right ones, since Ω (s) is the only subdomain with right hand side non null, then c

            K (s) ii K (s) ib-1 K (s) ib K (s) b-1i K (s) b-1b-1 + S (-) b-1b-1 0 K (s) bi 0 K (s) bb + S (+) bb                         x (s) i x (s) b-1 x (s) b             =             f (s) i f (s) b-1 f (s) b             (1.85)
In this case, we have the following theorem Theorem 1.2. In a case of a one way splitting, the Jacobi iterative algorithm for the two-Lagrange multiplier FETI method with augmented term equal to the complete outer Schur complement (as in Equation 1.85) converges in a number of iteration equal to the number of subdomains minus one.

Proof. See Roux et al. [START_REF] Roux | Optimization of Interface Operator Based on Algebraic Approach[END_REF].

With this theorem we know that in every subdomain the optimal augmentation term, in a one way split, correspond to the Schur complement of the outer domain, which in practice is impossible to compute. Furthermore this matrix is a full matrix, even if it was available for free, using it in the generalized Robin boundary condition would make the local problem too expensive to solve so it is mandatory to find sparse approximations.

In this context, several approaches had been done [START_REF] Roux | Approximation of Optimal Interface Boundary Conditions for Two-Lagrange Multiplier FETI Method[END_REF], from which we can name the following:

Schur neighbour: In a physical context, the complete outer Schur complement matrix represent the interactions of all the degree of freedom of the subdomains, condensed on the interfaces. This condensation only with the neighboring subdomains, leads to approximate the complete outer Schur complement with the neighbour Schur complement. In this case the cost of computation and the exchange of data are reduced to the neighboring subdomains only. Even though this is a good theoretical approximation, in practice it still has a very expensive cost. From this first approach defines the next two approximations, that try to mimic the behavior of the neighbour Schur complement, but with an sparse and not expensive approximation of it, that gives good convergence for the FETI-2LM method.

Lumped approximation: It has been shown in [START_REF] Roux | Optimization of Interface Operator Based on Algebraic Approach[END_REF] that an approximation of the neighbor Schur complement matrix

K (s) bb -K (s) bi K (s) ii -1 K (s)
ib with its first term K (s) bb gives good results. This way to approximate, is very easy to implement since this matrix is computed by the neighboring subdomain during the assembly procedure and the integration of the contribution of the interface nodes. Only one exchange with the neighboring subdomain is required for this regularization procedure, the results are however not as good as we would hope, so a third and better approximation will be finally used.

Sparse approximation based on patches along the interface:

This sparse approximation for the neighboring Schur complement, leads to better results than the lumped approximation, as shown in [START_REF] Roux | Optimization of Interface Operator Based on Algebraic Approach[END_REF] and is the only one used later in the numerical results section. The goal is emulate the spectral density of the neighbour Schur complement matrix as close as possible.

A dense, but cheaper approximation can be made if the condensation process is done only considering a neighbouring area of the interface, that is, the Schur complement is computed by

S(s) = K (s) bb -K (s) bj K (s) -1 jj K (s) jb , s = 1, 2, ∀j ∈ V d (1.86)
where the nodes j are no longer all the interior nodes, but instead they belong to the set V d = { indexes of interior nodes, such that the minimum connectivity distance between each of these nodes and the nodes belonging to the interface is lower than d, d ∈ N}.

Even though this procedure reduces the computational cost of the Schur complement, it generates a dense augmentation matrix, which adds an extra cost to the forward-backward substitution in the method. However its from this idea that the sparse approximation is build.

We start again with the 2 subdomain case, as the general case is direct. The idea is to perform the condensation only in a small part of the interface, called patch [START_REF] Roux | Approximation of Optimal Interface Boundary Conditions for Two-Lagrange Multiplier FETI Method[END_REF]. Lets define the following subsets of indexes for the nodes of This approximation consists on defining an sparse augmentation matrix, in order to avoid the increase on the bandwidth of the local problem matrix. This augmentation matrix is build by local condensation along the interface and extraction of some coefficients. The algorithm to compute the augmentation matrix A (q) b , q = 1, 2 is described in 7.

Ω (s) , s = 1, 2: V Ω (s) = { indexes of nodes inside the subdomain Ω (s) } V Γ = { indexes of nodes inside the interface Γ } V j p = {
Algorithm 7 Sparse approximation of neighbour Schur complement 1: Initialization of p and d.

2: Construction of the sparse structure of the interface matrix

A (q) b ∈ R dimV Γ ×dimV Γ . 3:
Construction of the sparse structure of the subdomain matrix K (s) ∈ R dimV Ω (s) ×dimV Ω (s) . 4: Assembly of the matrix K (s) . 5: for j in V Γ do Computation of the dense matrix S(s) by condensation of the matrix K(s) on the patch V j p .

8:

Extraction of the coefficients of the line associated with the node j from the matrix S(s) and insertion inside the matrix A (q) b at the line associated with the node j. 9: end for 10: Construction of the symmetric matrix

A (q) b = 1 2 A (q) T b + A (q)
b .

The augmentation matrix build up to the last line of the algorithm is non symmetric, so we symmetrize the matrix to avoid this drawback.

The same algorithm is used in any other subdomain in the general case, considering one neighbor at the time to build one augmentation matrix A (s) Γ j for every neighbour that will be added to the corresponding local matrix K (s) . As an example the regular mesh with Q 1 -finite elements is presented in Figure 1.6. This numbering leads to the definition, for the node 13, with p = 1 and d = 1, 2, of the subsets of indexes V 13 1 = {7, 13, 19} V 13 1,1 With the previous computation of the augmentation matrix, the final algorithm is the one described theoretically in previous section Algorithm 6, but without any projection nor preconditioning, since the local problems are all well posed. The full orthogonalization process is kept as it is needed to avoid the loss of orthogonality produced when multiplying by the FETI operator.

In practice, even if the algorithms are similar between this method and the non symmetric version of FETI-1LM, we have to point out that most of the implementation differences are in the construction of local problems, since we need to add the computed A (s) matrix. The other big difference is in the computation by this new operator, this is done by changing the assembly of the residuals according to the definition in Equation 1.77. 1.3 New FETI as an hybrid between one and two lagrange methods

Development

Now we are ready to present a new method that arises "directly" if we think on the properties of each previous method and the fact that both share similarities in their implementations. The main idea is to mix the good convergence and speed of the FETI-1LM method with the robustness of the FETI-2LM for the more complex scenarios that usually we have to face in real life applications.

To show this we can see in Figure 1.8 a contact problem coming from the basic scheme model of a rolling bearing, this kind of problems can be easily extended to other contact schemes were we want to compute the interactions between objects of different materials, that can be as different as steel and rubber. The heterogeneous problem with contact boundary conditions imposed with the penalty method is a type of problem that can make the discretization matrix of the model and consequently the FETI operator, very ill conditioned.

In this configuration a natural interface will be the one that divides the two materials, if this is the case, then the FETI-1LM method with the extended preconditioner will be appropriate to handle the heterogenities. However, this is not always the case when you do an automatic partitioning of the global mesh,

P Figure 1.8 -A contact problem.
or if, for example, we use different discretizations that leads to different types of meshing, then we will have the case of a non conforming mesh. Other cases were anisotropic materials needs to be modeled, will also give bad numerical features on the interface. For those problems the FETI-1LM will not assure the usual fast convergence. The FETI-2LM method is in general a more robust method, due to the formulation of two independent lagrange multipliers, that can surpass this numerical issues.

If we think now in the interfaces that are completely contained in one of the two material, or in general do not present numerical issues, the regular FETI method will be more suited as we will get a faster convergence.

Using one or the other will gives us either fast or accurate results, but we can try to get both of these good features if we choose properly which method to apply in every interface. If we manage to make the choice that better suits every local interface, then we will most likely have a new hybrid method with both good qualities. This method we will call it Hybrid-FETI and it's formulation can be described mathematically as previous methods, but this time starting from a different basic configuration, as we need at least two different local interfaces. For simplicity reason the basic problem consist in three subdomains and two interfaces that are completely independent of one another and no crosspoints as shown in Figure 1.9. In any case the generalization is straightforward since the interfaces are always disconnected with FETI approaches.

Ω Γ Ω Γ Ω 1 2 3 1LM 2LM Figure 1.9 -Three subdomain divisions.
Using the subscript notation of previous sections, we have i for interior nodes and b for the ones on the interface Γ = Γ 1 ∪ Γ 2 , we have that for each subdomain Ω (s) , s = 1, 2, 3 the contribution to the matrix and right hand side of a finite element discretization.

K (s) =        K (s) ii K (s) ib K (s) bi K (s) bb        , f (s) =        f (s) i f (s) b        (1.89)
The global discretization matrix of the form Kx = f (as in 1.16) can be written from the assembling of the local problems

                  K (1) ii 0 0 K (1) ib 0 K (2) ii 0 K (2) ib 0 0 K (3) ii K (3) ib K (1) bi K (2) bi K (3) bi K bb                                     x (1) i x (2) i x (3) i x b                   =                   f (1) i f (2) i f (3) i f b                   (1.90)
up to here we have the same as previous methods, but now we need to look in detail the interactions that happen in this particular interface, so we can separate it and impose the different boundary conditions, to form the correspondent local problems for each method.

If we use the subscript b 1 for the interface Γ 1 in which we will impose the Neumann condition as in the 1LM method and the subscript b 2 for the 2LM interface where a Robin condition is imposed, then the previous local contributions are

K (1) bi =        K (1) b 1 i 0        , K (2) bi =        K (2) b 1 i K (2) b 2 i        , K (3) bi =        0 K (3) b 2 i        K (1) ib = K (1) ib 1 0 , K (2) ib = K (2) ib 1 K (2) ib 2 , K (3) ib = 0 K (3) ib 2 (1.91)
As for the interaction between the interfaces, we have

K bb = K (1) bb + K (2) bb + K (3) bb =        K (1) b 1 b 1 0 0 0        +        K (2) b 1 b 1 0 0 K (2) b 2 b 2        +        0 0 0 K (3) b 2 b 2        f b = f (1) b + f (2) b + f (3) b =        f (1) b 1 0        +        f (2) b 1 f (2) b 2        +        0 f (3) b 2        (1.92)
And the unknown displacements are

x b =       x b 1 x b 2       (1.93)
With the previous separation, we can form the local Neumann, Robin and hybrid problems to solve in each subdomain. For the first problem, we introduce a single langrage multiplier to have

       K (1) ii K (1) ib 1 K (1) b 1 i K (1) b 1 b 1               x (1) i x (1) b 1        =        f (1) i f (1) b 1 + λ (1) b 1        (1.94)
For the second subdomain, we form the hybrid local problem, where the lagrange multiplier introduced for the first subdomain is shared in the "left" interface, an we add 2 extra at the interface in the "right" and the augmentation matrix, as in the 2LM method, with this we have

             K (2) ii K (2) ib 1 K (2) ib 2 K (2) b 1 i K (2) b 1 b 1 0 K (2) b 2 i 0 K (2) b 2 b 2 + A (2) b 2                           x (2) i x (2) b 1 x (2) b 2              =              f (2) i f (2) b 1 + λ (2) b 1 f (2) b 2 + λ (2) b 2              (1.95)
And for the third subdomain,

       K (3) ii K (3) ib 2 K (3) b 2 i K (3) b 2 b 2 + A (3) b 2               x (3) i x (3) b 2        =        f (3) i f (3) b 2 + λ (3) b 2        (1.96)
Since in the interface between subdomains 1 and 2 we have one lagrange multiplier, with opposite values stored in each subdomain, we will denote it as

λ b1 = λ (1) b1 = -λ (2) b1
(1.97) Now, using the same condensation process as previous methods, we will use the explicit relation between the multipliers and the displacements to eliminate the unknowns of the inner nodes of each subdomain. Lets take the first subdomain, where we know from the first equation in 1.94 that

x (1) i = -K (1) -1 ii K (1) ib 1 x (1) b 1 + K (1) -1 ii f (1) i (1.98)
introducing it in the second equation we have

λ b 1 = K (1) b 1 i x (1) i + K (1) b 1 b 1 x (1) b 1 -b (1) f 1 = -K (1) b 1 i K (1) -1 ii K (1) ib 1 x (1) b 1 + K (1) b 1 b 1 x (1) b 1 -f (1) b 1 + K (1) b 1 i K (1) -1 ii f (1) i = S (1) b 1 b 1 x (1) b 1 -c (1) b 1 (1.99) here S (1) b 1 b 1 := K (1) b 1 b 1 -K (1) b 1 i K (1) -1 ii K (1) ib 1
is the Schur complement of the domain 1 in the interface 1LM and c

(1) b 1 := f (1) b 1 -K (1) b 1 i K (1) -1 ii f (1)
i . With this we can have the displacement as a function of λ b1

x

(1) b 1 = F (1) λ b 1 + F (1) c (1) b 1 (1.100)
where

F (1) := S (1) -1 b 1 b 1 (1.101)
is the inverse of the Schur complement. For the second subdomain, we have from the first line in 1.95 that

x (2) i = -K (2) -1 ii K (2) ib 1 x (2) b 1 -K (2) -1 ii K (2) ib 2 x (2) b 2 + K (2) -1 ii f (2) i (1.102)
replacing in the second and third equation, we have

       λ b 1 λ (2) b 2        =        K (2) b 1 i x (2) i + K (2) b 1 b 1 x (2) b 1 -b (2) b 1 K (2) b 2 i x (2) i + (K (2) b 2 b 2 + A (2) b 2 )x (2) b 2 -b (2) b 2        =         K (2) b 1 b 1 -K (2) b 1 i K (2) -1 ii K (2) ib 1 -K (2) b 1 i K (2) -1 ii K (2) ib 2 -K (2) b 2 i K (2) -1 ii K (2) ib 1 K (2) b 2 b 2 -K (2) b 2 i K (2) -1 ii K (2) ib 2 + A (2) b 2                 x (2) b 1 x (2) b 2         +         -f (2) b 1 + K (2) b 1 i K (2) -1 ii f (2) i -f (2) b 2 + K (2) b 1 i K (2) -1 ii f (2) i         =        S (2) b 1 b 1 S (2) b 1 b 2 S (2) b 2 b 1 S (2) b 2 b 2 + A (2) b 2               x (2) b 1 x (2) b 2        -        c (2) b 1 c (2) b 2        (1.103)
where again the S terms denotes the corresponding Schur complements of the subdomain 2 in the interfaces 1LM and 2LM and c

(2) b 1 := f (2) b 1 -K (2) b 1 i K (2) -1 ii f (2) i , c (2) b 2 = f (2) b 2 -K (2) b 1 i K (2) -1 ii f (2)
i . So both the displacements of both local interfaces are

       x (2) b 1 x (2) b 2        = F (2)        λ b 1 λ (2) b 2        + F (2)        c (2) b 1 c (2) b 2        (1.104)
here we have

F (2) :=        S (2) b 1 b 1 S (2) b 1 b 2 S (2) b 2 b 1 S (2) b 2 b 2 + A (2) b 2        -1 (1.105)
We can divide F (2) in blocks, such that

F (2) :=        F (2) b 1 b 1 F (2) b 1 b 2 F (2) b 2 b 1 F (2) b 2 b 2        (1.106)
then we can also build a separate explicit relation for both displacements in this subdomain

       x (2) b 1 x (2) b 2        =        F (2) b 1 b 1 λ b 1 + F (2) b 1 b 2 λ (2) b 2 + F (2) b 1 b 1 c (2) b 1 + F (2) b 1 b 2 c (2) b 2 F (2) b 2 b 1 λ b 1 + F (2) b 2 b 2 λ (2) b 2 + F (2) b 2 b 1 c (2) b 1 + F (2) b 2 b 2 c (2) b 2        (1.107)
For the third subdomain, we do the same previous treatment to obtain first

x (3) i = -K (3) -1 ii K (3) ib 2 x (3) b 2 + K (3) -1 ii f (3) i (1.108)
and then replace it in the second equation of 1.96

λ (3) b 2 = K (3) b 2 i x (3) i + (K (3) b 2 b 2 + A (3) b 2 )x (3) b 2 -f (3) b 2 = -K (3) b 2 i K (3) -1 ii K (3) ib 2 x (3) b 2 + (K (3) b 2 b 2 + A (3) b 2 )x (3) b 2 -f (3) b 2 + K (3) b 2 i K (3) -1 ii f (3) i = (S (3) b 2 b 2 + A (3) b 2 )x (3) b 2 -c (3) b 2 (1.109)
With analogous definitions of each term, as before. The equivalence in this case is given by the next equation

x (3) b 2 = F (3) λ (3) b 2 + F (3) c (3) b 2 (1.110)
where

F (3) := (S (3) b 2 b 2 + A (3) b 2 ) -1 (1.111)
is the inverse of the Schur complements plus the augmentation term in this subdomain.

The continuity conditions needed in every case will vary depending on the interface. For the interface 1LM, eliminating the jump of the interface solution is enough to achieve this, meaning

x (1) b 1 -x (2) b 1 = 0 (1.112)
but for the interface 2LM, since the interface condition is not consistent, as seen in 1.69 we need to add another constraint that gives us a residual of the form

x (2) b 2 -x (3) b 2 = 0 λ (2) b 2 + λ (3) b 2 -A (2) b 2 x (2) b 2 -A (3) b 2 x (3) b 2 = 0 (1.113)
this mixed form is recombined as in the 2LM method so together with Equation 1.112 we obtain the conditions needed to impose continuity across complete the interface

x (1) b 1 -x (2) b 1 = 0 λ (2) b 2 + λ (3) b 2 -(A (2) b 2 + A (3) b 2 )x (2) b 2 = 0 λ (2) b 2 + λ (3) b 2 -(A (2) b 2 + A (3) b 2 )x (3) b 2 = 0 (1.114)
With this equations plus the relations between the lagrange multipliers and the different displacement, coming from the solution of local problems 1.100,1.107,1.110, we can do the condensation of previous problems on to the interface for the hybrid method

                F (1) + F (2) b 1 b 1 -F (2) b 1 b 2 0 A (2) b 2 + A (3) b 2 F (2) b 2 b 1 I -A (2) b 2 + A (3) b 2 F (2) b 2 b 2 I 0 I I -A (2) b 2 + A (3) b 2 F (3)                                λ b 1 λ (2) b 2 λ (3) b 2                =                 -F (1) c (1) b 1 + F (2) b 1 b 1 c (2) b 1 + F (2) b 1 b 2 c (2) b 2 A (2) b 2 + A (3) b 2 F (2) b 2 b 1 c (2) b 1 + F (2) b 2 b 2 c (2) b 1 A (2) b 2 + A (3) b 2 F (3) c (3) b 2                 (1.115)
This matrix is not symmetric, so an ORTHODIR method with full reconjugation is used to solve it. We name the previous method Hybrid-FETI as in a hybridization between FETI-1LM and FETI-2LM.

Extension to a general problem

For the case of a general problem with N s > 2, we proceed as in the 2LM method, and we divide the total interface into neighbouring interface edges so that each one is treated depending on the chosen method.

We split the total interface nodes into this interface edges, that correspond to a collection of connecting interface nodes, denoted Γ j as we did for the 2LM case.

We also add a unique marker to this edges in order to know a priori which one correspond to each method, we denote the marked edges as Γ j 1lm and Γ j 2lm for both FETI-1LM and FETI-2LM methods respectively. Crosspoints in particular are edges of one node, were we can arbitrarily marked them 1LM or 2LM, we will see the impact this choice in the numerical experiments.

We also define the boolean matrix B (s) Γ j for subdomain Ω (s) as the restriction of a vector defined in Ω(s) to its values on the interface Γ j , i.e

B (s) Γ j v (s) = v (s) b | Γ j (1.116) again v (s) =        v (s) i v (s) b      
 with the notation for interior and interfaces nodes being i and b respectively.

With this definitions we can write the local problems for subdomain Ω (s) as well as the conditions to get the continuity, for the general case

           K (s) + Γ j 2lm ⊆Ω (s) B (s) T Γ j 2lm A (s) Γ j 2lm B (s) Γ j 2lm            x (s) = f (s) - Γ l ⊆Ω (s) B (s) T Γ l λ (s) Γ l , s = 1, . . . , N s B (s) Γ j x (s) -x b | Γ j = 0 ∀Γ j , Ω (s) ⊇ Γ j Ω (q) ⊇Γ j 2lm λ (q) Γ j 2lm -A (q) Γ j 2lm x b | Γ j 2lm = 0 ∀Γ j 2lm (1.117) where A (s) Γ j 2lm
is the augmentation matrix in subdomain Ω(s) for the interface Γ j 2lm , the other terms are defined analogous as for the case in FETI-2LM. Lets also define the total augmentation matrix for a subdomain that contains an edge "2lm", where if not the case, its value is zero

A (s) 2lm := Γ j 2lm ⊆Ω (s) B (s) T Γ j 2lm A (s) Γ j 2lm B (s) Γ j 2lm (1.118)
With this equation we can do the condensation and obtain the problem to solve in the interface. Since, in general, two methods are used, then equations in the interface will also need to be separated, finally we have the interface problems

Ω (q) ⊇Γ j 2lm         λ (q) Γ j 2lm -A (q) Γ j 2lm B (s) Γ j 2lm K (s) + A (s) 2lm -1 Γ l ⊆Ω (s) B (s) T Γ l λ (s) Γ l         = - Ω (q) ⊇Γ j 2lm A (q) Γ j 2lm B (s) Γ j 2lm K (s) + A (s) 2lm -1 f (s) ∀Γ j 2lm , ∀Ω (s) ⊇ Γ j 2lm (1.119) B (s) Γ j 1lm K (s) + A (s) 2lm +           λ Γ j 1lm + Γ l 2lm ⊆Ω (s) B (s) T Γ l 2lm λ (s) Γ l 2lm           + B (s) Γ j 1lm R (s) α (s) + B (q) Γ j 1lm K (q) + A (q) 2lm +           λ Γ j 1lm - Γ l 2lm ⊆Ω (q) B (q) T Γ l 2lm λ (q) Γ l 2lm           -B (q) Γ j 1lm R (q) α (q) = -B (s) Γ j 1lm K (s) + A (s) 2lm + f (s) + B (q) Γ j 1lm K (q) + A (q) 2lm + f (q) ∀Γ j 1lm = Ω (s) ∩ Ω (q) (1.120)
where the terms R (s) (or q) is a kernel generator for the matrix K (s) with some coefficients α (s) such that the admissibility condition 1.36 for this "floating domain" is fulfilled. We need then to use the same projection strategy explained in section 1.1.2, in order to treat this subdomains. In the case that the subdomain contains a 2LM interface large enough (for example, two or more nodes for a 2D elasticity problem, three for 3D elasticity, etc) the term A (s) 2lm makes that there is no longer a kernel, so R (s) is zero.

This interface problem can also be written as Fλ = d with F being non symmetric, so we use ORTHODIR version with left preconditioner and full reconjugation. We will talk about this preconditioner in the next section.

Preconditioner

If we want this method to be scalable for the mesh size, we need to find a preconditioner as in FETI-1LM. Up to this day no optimal preconditioner as been found for the FETI-2LM, but this method has been tested and shown numerically to be scalable for the mesh size on his own [START_REF] Farhat | Two-level domain decomposition methods with Lagrange multipliers for the fast iterative solution of acoustic scattering problems[END_REF]. No extra computation needs to be done for the interface edges where the Robin condition is imposed. On the other hand, for those where one lagrange multiplier is used, a preconditioner is mandatory in order to keep the convergence independent of the mesh size.

We will base the preconditioner on local information, just as the FETI-1LM does. If we look at the Dirichlet preconditioner from a mechanical point of view, we know it consist on solving a local problem whose boundary condition is the jump in the displacements fields along this subdomain interfaces. The solution of this problem allows to compute the corresponding interface traction forces and therefore a correction to the FETI operator.

From this process we built the Dirichlet preconditioner, as already detailed in subsection 1.1.3, by applying it to all the interface edges and then adding this local contributions.

The local problem to solve, then comes from solving the following Dirichlet problem in

Ω (s)       K (s) ii K (s) ib 0 I              w (s) i w (s) b        =       0 {w (s) b }       (1.121)
the subscript i is for interior nodes and b for the interface, so the jumps are imposed along the total local interfaces.

We compute the multiplication of a vector w restricted to the interface edges Γ j 1lm by our preconditioner as

D -1 Γ j 1lm w| Γ j 1lm = B (s) Γ j 1lm S (s) bb B (s) T Γ j 1lm + B (q) Γ j 1lm S (q) bb B (q) T Γ j 1lm w| Γ j 1lm , ∀Γ j 1lm : Γ j 1lm = ∂Ω (s) ∩ ∂Ω (q) (1.122)
where S (s) bb (or q) is the Schur complement of the subdomains interior nodes into the local interface and the matrices B (s) are the signed assembly matrices, but can also be taken as their generalisation for heterogeneous problems, i.e B(s) = B (s) β (s) with β (s) some scaling matrix.

The multiplication by this operator are made in an analogous way to the case 1LM, first computing the Schur complements in each side of the interface, then making the exchange with the neighbour subdomain. We can write the previous preconditioner as a global matrix of the form

D -1 w =            Γ j 1lm ⊆Γ D -1 Γ j 1lm            w (1.123)
with w any vector defined in the interface Γ .

One problem that we encounter in this preconditioner, happens when in some subdomain the local interface is composed of both 1LM and 2LM interface edges. In this case when we apply the previous preconditioner, we are also imposing a jump in the 2LM edges, but not of the same type, as a Robin condition is used here. In Figure 1.10 we can see the two different boundary conditions that were used, the upper one corresponding to the usual Dirichlet preconditioner, and the bottom one represent an option to build a new preconditioner.

The advantage of this boundary conditions is that now we obtain the exact local inverse of the Hybrid-FETI operator. The preconditioner is then changed to solve this problem in the subdomains with shared edges. It involves the computation of a different Schur complement, this time with the augmentation matrix that allows to mimic the behaviour of the local interface marked as 2LM and therefore of the whole subdomain. subdomain with at least one edge marked as 1LM

            K (s) ii K (s) ib 1lm K (s) ib 2lm 0 I 0 0 0 I                          w (s) i w (s) b 1lm w (s) b 2lm              =             0 {w (s) b 1lm } 0             (1.124)
where b 1lm , b 2lm are subscript for the union of local edges marked as 1LM and 2LM respectively. Here we impose the computed jump only in the 1LM edges and zero in the 2LM. Then we compute the respective traction forces by multiplying against the matrix of the local problem. This leads to the computation of the following Schur complement

S (s) bb =        K (s) b 1lm b 1lm 0 0 K (s) b 2lm b 2lm + A (s) b 2lm        -K (s) b 1lm i K (s) b 2lm i K (s) -1 ii        K (s) ib 1lm K (s) ib 2lm        =        S(s) b 1lm b 1lm S (s) b 1lm b 2lm S (s) b 2lm b 1lm S (s) b 2lm b 2lm + A (s) b 2lm        (1.125)
The structure of this matrix is the same as usual and correspond to the local inverse of the condensed interface problem (as seen in 1.105).

With this computation of the Schur complement, the preconditioner for every interface edge Γ j with one lagrange multiplier is defined the same way as previous, but remembering that the Schur complement now includes the augmentation matrices (whenever they exist)

R -1 Γ j 1lm w| Γ j 1lm = B (s) Γ j 1lm S(s) bb B (s) T Γ j 1lm + B (q) Γ j 1lm S(q) bb B (q) T Γ j 1lm w| Γ j 1lm , ∀Γ j 1lm : Γ j 1lm = ∂Ω (s) ∩ ∂Ω (q) (1.126)
the definition of every term, other than the Schur complement, are the same as previous (including the consideration of the B matrices as scaling ones). This preconditioner can also be written in the form of a global matrix denoted R -1 as in "Robin" preconditioner

R -1 w =            Γ j 1lm ⊆Γ R -1 Γ j 1lm            w (1.127)

Implementation

The parallel computation of the FETI methods is done by subdomains structures that are used to store the residuals of the FETI methods in general. This structures for both FETI-1LM and FETI-2LM are basically the same, so in practical terms even if we use 1 lagrange multiplier (gradient) for the interface, the values of this vector will be stored in both subdomains sharing the local edge (actually one saves λ| Γ j and the other -λ| Γ j ).

Depending on the method this computation of the residuals involves the resolution of the corresponding local Neumann and Robin problems in order to compute x (i) b . In general for our method there will be at least one subdomain where the local problem solved is the one that comes from imposing both conditions at the same time.

Lets recall how the gradients are computed for each method. To do so, lets take as example the division of a domain into two subdomains Ω (1) , Ω (2) with an interface Γ b . We take this simple case because the procedure is the same for every interface edge in a more general partition.

The preconditioned gradient to compute for the 1LM method, is then

g = D -1 (Fλ -d) = D -1 (x (2) b -x (1) b ) (1.128) where x (s)
b , s = 1, 2 are the displacement of the interface nodes of each subdomain. No projection will be shown here as the implementation is analogous to the 1LM method and applies only to subdomains where the local interface is completely marked as 1LM.

If we look at previous computation from an implementation point of view, we have that for each process that computes each subdomain, the gradient is obtained as

g =        g (1) b g (2) b        = D -1        x (2) b -x (1) b x (1) b -x (2) b        (1.129)
where the process s computes and stores g (s) b , s = 1, 2. We also note that g

(1) b = -g (2)
b as in this case the unknown is a single lagrange multiplier along the interface. Now for the 2LM method, the gradient is computed from equations 1.70, meaning that for each process we have

g =        g (1) b g (2) b        =        λ (1) b + λ (2) b -(A (1) b + A (2) b )x (2) b λ (2) b + λ (1) b -(A (1) b + A (2) b )x (1) b        (1.130) the values of g (1)
b and g (2) b are different, since 2 independent lagrange multipliers are used for this method, but we can see that the computation structure is the same as the FETI-1LM method. The difference here lies in the exchanges needed to compute both gradients, in the case of the 2LM method more than one exchange is needed.

The computation of the gradients is the base to the computation by the F operator, needed in the use of the iterative method. For the case of an arbitrary number of subdomains, the previous computations is extended by defining a division of the total interface into interface edges as previously done. So, the multiplication of a vector v defined in the interface by the FETI operator is done by:

• Solving the local problems K (s) + A (s) 2lm with the right hand side equal to

t (s) T B (s) v =       0 v (s) b       and computing x (s) b .
• Assembly of the corresponding gradient, passing the information from one subdomain to its neighbours The first part is similar for the FETI-1LM, FETI-2LM and our method, only changing the existence of the augmentation term when some 2LM marked edge is present. Is in the assembly part of the gradient where we need to do the mix of both method.

To explain this in detail, lets take the same case of a three subdomain division as Figure 1.9 with one lagrange multiplier to the "left" and two to the "right". If we want to solve it using the FETI-1LM method we will have the following preconditioned gradient at each iteration

g =                    g (1) b 1 g (2) b 1 g (2) b 2 g (3) b 2                    = D -1                    x (2) b 1 -x (1) b 1 x (1) b 1 -x (2) b 1 x (3) b 2 -x (2) b 2 x (2) b 2 -x (3) b 2                    = D -1 Γ b 1 + D -1 Γ b 2                    x (2) b 1 -x (1) b 1 x (1) b 1 -x (2) b 1 x (3) b 2 -x (2) b 2 x (2) b 2 -x (3) b 2                    (1.131)
where the multiplication by the preconditioner was explained in the subsection 1.1.3.

If we think on solving it using only the FETI-2LM method, then the gradient will be

g =                    g (1) b 1 g (2) b 1 g (2) b 2 g (3) b 2                    =                    λ (1) b 1 + λ (2) b 1 -(A (1) b 1 + A (2) b 1 )x (2) b 1 λ (2) b 1 + λ (1) b 1 -(A (2) b 1 + A (1) b 1 )x (1) b 1 λ (2) b 2 + λ (3) b 2 -(A (2) b 2 + A (3) b 2 )x (2) b 2 λ (3) b 2 + λ (2) b 2 -(A (3) b 2 + A (2) b 2 )x (1) b 2                    (1.132)
As for the preconditioned gradient of our method the computation is done by assembling

g =                    g (1) b 1 g (2) b 1 g (2) b 2 g (3) b 2                    = D -1 Γ b 1 1lm                    x (2) b 1 -x (1) b 1 x (1) b 1 -x (2) b 1 λ (2) b 2 + λ (3) b 2 -(A (1) b 2 + A (2) b 2 )x (2) b 2 λ (3) b 2 + λ (2) b 2 -(A (2) b 2 + A (1) b 2 )x (1) b 2                    (1.133) in here the preconditioner D -1 is reduced to D -1 Γ b 1 1lm
and affects the gradient only in the interface edge "1LM" leaving the last two lines untouched.

The previous computation of the residuals allows us to obtain in an analog way the multiplications by the FETI operator for each of the three methods, by changing only the corresponding gradient.

With this, the new method is complete and a formal algorithm can be described and implemented. As already said, the operator is non symmetric due to the fact that it shares blocks of the also non symmetric 2LM operator, so we use the iterative ORTHODIR method with left preconditioner and full reconjugation, for the same reasons as we use it in the unsymmetric and 2LM methods.

The FETI-1LM method for unsymmetric matrices described in Algorithm 6 is the base of our implementation. Both of them share theoretical resemblance, changing the definition of our operator F and also modifying D to D 1lm because when using the preconditioner, it is used only in the 1LM interface edges.

In practice the modifications are more complex so we are listing them

• Marking of the interface edges.

• Correct construction of the local problems, meaning that we add the augmentation matrix in the subdomains with an edge 2LM.

• Change in the computation of the residuals for the multiplication by our new operator.

• Modify the local Dirichlet problem to apply the different preconditioners and use them only in the 1LM interface edges.

• Project only in the subdomains with a complete local interface 1LM (this is done automatically, after computing the corresponding kernel in each subdomain).

In the next section we will test numerically our method to show the advantages of it in some kind of difficult problems.

Numerical results

In this section we will show a comparative performance of the methods presented in this chapter.

To begin with, we present the formulation of the Poisson problem in 3D with Dirichlet condition in a part of the boundary. Let Ω ⊆ R 3 be a bounded domain and let ∂Ω D ⊆ ∂Ω be a part of the boundary, where null Dirichlet conditions will be imposed. Given f ∈ L 2 (Ω), the problem is

Find u ∈ H 1 (Ω) such that              -ν∆u = f in Ω u = 0 on ∂Ω D ∂u ∂n = 0 on ∂Ω \ Ω D (1.134)
with ν ∈ R + , the parameter that will define the difference between materials and therefore in the local stiffness matrices of any possible divisions of Ω.

Let V := {v ∈ H 1 (Ω) : v| ∂Ω D = 0}. Then, the weak or variational formulation of previous problem is

Find u ∈ V such that ν Ω ∇u • ∇vdx = Ω f v, ∀v ∈ V (1.135)
This problem is discretized using tri-linear Q 1 finite elements functions, leading to already known global system of equations

Kx = f (1.136)
Figure 1.11 -Two material bar.

Two material bar

In our first test we will model a bar of two different materials under the gravity, see Figure 1.11, in this case

Ω = [0, 10] × [0, 1] × [0, 1], f = (0, 0, g) T . The boundary fixed is ∂Ω D = {(x, y, z) ∈ ∂Ω : x = 0 ∨ x = 10}.
The problem will consist on a bar with two blocks of a different material across the x-axis, they are located at the same distance of each boundary. The idea is to represent a localized heterogeneous zone in a general homogeneous domain, leaving only a few interfaces with a bad local conditioning. For each material we have the following values of their parameters ν 1 = 10 5 , "blue" blocks

ν 2 = 10 0 , otherwise (1.137)
Depending on the subdomain partition that we do, we change the size of the blocks to match the interfaces of the subdomains in the different divisions. Since the total domain is divided into slices, we assign one of values of ν to each subdomains, leaving only two in each case different from the rest, this way we only have heterogenities across four of the interfaces.

Before doing any comparative test against other FETI methods, we need to check the performance of both possible preconditioners for this method, For all cases in this chapter we use the following global stopping criterion

Kx p -f 2 f 2 < 10 -4 (1.138)
but we also use a second criterion, that is, the relative norm of the jump of the solution between subdomains [START_REF] Farhat | Implicit Parallel Processing in Structural Mechanics[END_REF].

x p -x p-1 2 x p 2 < 10 -4 (1.139)
The results are shown in Table 1.1. In this case there is no significant difference between both of them, even in some cases the Dirichlet preconditioner performs better than the Robin. We think that this is due to the fact that the small number of interfaces where the preconditioner is different is not high enough to let the differences to be clear. In examples where more subdomains share the two types of interface the difference is more evident (See the contact case 1.4.2).

We also note, that almost no augmentation in the number of iterations is present, even if we have increased the number of subdomains, this good quality is explained because the most part of the interfaces were marked as 1LM, meaning that we are very close to the regular FETI-1LM method. We hope in any case to keep this property even in the harder problems.

We will see now, how the Hybrid method performs when increasing the number of elements per subdomains, we hope that also some of the other properties of the methods FETI-1LM and FETI-2LM are hold in the Hybrid case, in particular the property of having a small increase in the number of iterations when the number of elements in a subdomain is augmented. Hence, we test different values of h (mesh size) for a constant of 64 subdomains, the results are summarized in Figure 1.12.

From Table 1.1 and Figure 1.12 we can see that, as expected, the number of iterations after augmenting the local size of the problems or the number of subdomains does not change considerably. This is due to the fact that this method is build from a combination of two numerically scalable methods.

One of the main issues of the new method, lies in the problem of which interfaces should be marked as 1LM or 2LM. The intuition tells that the interfaces where the materials change should be marked as 2LM and the rest as 1LM. We will now test this choice against a different one, where two extra interfaces are marked as 2LM, as seen in Figure 1.13. The idea is to "cover" the problematic interfaces, in this case is the jump of materials, but in other cases this interfaces will not be so clear, so extending the 2LM marked interfaces will result mainly in a more robust method, as the 1LM interfaces will not always achieve convergence, but it will be probably slower.

1LM 1LM 1LM 1LM 2LM 2LM 1LM 1LM a) b) 1LM 2LM 2LM 1LM 2LM 2LM 1LM 1LM
For the rest of this basic tests, no difference will be done between preconditioners as they are minimal.

The results can be seen in the Table 1.2 and, in terms of speed, keeping the number of 2LM interfaces at a minimum will improve the convergence rate. Given that this it not completely clear for every case, the extension of 2LM interfaces, in order to "cover" the difficult interfaces, will give us a more robust method.

In order to achieve a better performance, in the following for the Hybrid-FETI case we will mark as 2LM the interface that connects different materials and all the rest will be treated as 1LM. With the defaults parameters for the Hybrid method, now we can perform a comparative against the two base method. In general the problem we are solving in this section cause no performance impact to the FETI-1LM with the basic Dirichlet preconditioner, however what we hope to obtain in this case is an improvement against the behaviour of the FETI-2LM. The idea is to show later some examples where we only have a good convergence for the FETI-2LM method, in which we will also hope to be a better alternative.

We modify the configuration of previous case to solve something closer to the type of problems where we want to apply this method. Hence, only a different block is considered, located in the center of the structure, with different parameters of value

ν 1 = 10 5 , center block ν 2 = 10 0 , otherwise (1.140)
With this considerations in Table 1.3 we can see the results in term of iterations for the different methods.

As expected, the symmetric FETI-1LM is the most suited method for this type of problem, however we wanted to show the non symmetric version of this method because the problem solved when working with contact models (one of our main objectives) is non symmetric, and in this case the performance is comparable to the one of Hybrid-FETI. As for the FETI-2LM, the very poor performance is explained by the type of divisions used, in this particular one, since the information is transmitted neighbour by neighbour, in order for it to cross the whole structure, we need at least the number of subdomains before achieving convergence.

Even if the FETI-2LM is at it worst, this case represents a good example of the problems we may find in applications, with two blocks separated by a small part with really bad conditioning (contact, non-conforming meshes, etc). Here the performance of the Hybrid method clearly outperforms that of FETI-2LM, meaning that we have ameliorated in terms of speed, the more robust method, which will be our next goal when we try to solve problems where the convergence for the FETI-1LM method is not assured.

Contact Problem

Our next test case will try to expose the advantages of this new method against both base ones. To achieve this, we will present and solve a contact problem.

Contact problems are often found in the structural engineering context, particularly when analyzing the assembly of different substructures. They are characterized by a non-penetration condition in an active area of contact which is not know a priori, this condition is modeled as constraint of a minimization problem. For these reasons, these problems may lead to stiff non-linear systems of equations.

Several simplification to the general contact models can be done to reduce the computational complexity of the solvers needed. In this line, we can name a still large class of contact problems that are characterized by the fact that one of the contacting bodies is much more stiffer than the other. Also, considering that the displacement of the contact body is constrained in only one direction we then have the so called unilateral contact problems. We also assume that the displacements are small and the deformation is linearly elastic (this is for presentation purposes, as in practice a large displacements model will be used). The final assumption to do is to not consider the effects of friction, which will simplify the contact constraints that we will establish later.

We start the formulation of this simplified contact model by recalling the 2D linear elasticity problem in a standard finite element framework. The equations that model this problem can be written as

-div(σ (u)) = f, in Ω u = 0 in ∂Ω D σ (u) • n = 0 in ∂Ω \ ∂Ω D (1.141) where σ ij (u) = 2µε ij (u) + λδ ij div(u) ε ij (u) = 1 2 ∂u i ∂x j + ∂u j ∂x i µ = E 2(1 + ν) λ = Eν (1 + ν)(1 -2ν) (1.142)
in this context µ and λ are the Lamé coefficients, E and ν are the Young's Modulus and Poisson ratio respectively. The different materials are defined by their values of E and ν, variation of this values in a generated structure will create the heterogeneities along or across the interfaces. This model correspond to the static linear elasticity, however for the next part, the displacements and several other components of the model will have a temporal dependency that we will not note, as we will just be testing the Hybrid method to solve this standard static model in a "fixed" final time T .

Lets present one of the classic mathematical formulations of a contact problem. We start as usual with a domain Ω that denotes the initial state of a linear elastic body, ∂Ω C will be the part of its boundary that is a potential area of Let u denotes the displacement field (in one direction), the discontinuity of u in the normal direction at point M ∈ ∂Ω C can be written as

t n =0 t n ≤0 δu = c 0 δu ≤ c 0 ∂Ω C ∂Ω C ' M M' n C ' n C a) c 0 (M) ∂Ω C ∂Ω C ' M M' n C b) ∂Ω C M n C c) t n =t n . t=σ(M)n
δu(M) = u(M) • n C + u(M ) • n C = u(M) • n C -u(M ) • n C = (u(M) -u(M )) • n C (1.143)
we write the non-penetration condition in the next form

δu(M) -c(M) 0 ≤ 0 (1.144)
where c(M) 0 denotes the initial gap, that we assume small, at the point M in Ω (see center image of Figure 1.14).

Let t = σ (M)n be the traction vector at point M, and t n = t • n the normal component of t. If ∂Ω C and ∂Ω C are in contact, the pressure at M is positive, which implies that t n ≤ 0 (see right image in Figure 1.

14).

As already said, we assume unilateral conditions for the contact problem, which leads to the Signorini -Fichera type of conditions, defined by:

• If δu(M) -c(M) 0 =
0 then t n ≤ 0 and the contact is said to be active.

• If δu(M)c(M) 0 < 0 then t n = 0 and the contact is said to be inactive. This relations can be rewritten equivalently as conditions on ∂Ω C , so for the unilateral frictionless contact problem, we have

t n ≤ 0 δu(M) -c(M) 0 ≤ 0 t n • (δu(M) -c(M) 0 ) = 0 (1.145)
In a finite element framework, we need to state this contact model in a different way, hence we present now one of the variational formulations of this problem. In this case we will formulate the previous model as a minimization problem with variational inequalities.

We define the following spaces

V = {v ∈ (H 1 (Ω)) 2 : v| ∂Ω D = 0} (1.146)
and

U = {v ∈ V : v • n -c(M) 0 ≤ 0 on ∂Ω C } (1.147)
We define then the bilinear and linear forms

a(u, v) := 2 Ω µε(u) : ε(v)dx + Ω λ(∇ • u)(∇ • v)dx b(v) := Ω f, v dx (1.148)
where

f, v = 2 i=1 f i v i (1.149)
and finally

F(v) = 1 2 a(v, v) -b(v) (1.150)
then for a given force f ∈ V the variational formulation of our contact problem is

Find u ∈ U such that a(u, v -u) ≥ f (v -u), ∀v ∈ U (1.151) or F(u) ≤ F(v) ∀v ∈ U (1.152)
With this presentation, a minimization method can be applied to solve this constrained problem.

For the solution of contact problems, we can not forget that two-types of nonlinearities exist, the geometrical ones and in the materials. In the geometric nonlinearities we encounter the ones derived from large displacements, deformation and contact constraints. An iterative incremental solution procedure is used then to obtain non-linear solution as a series of linear ones.

Discretization of this problem via finite elements, plus the Penalty method or the Lagrange multiplier method [START_REF] Zhong | Static Contact Problems -A review[END_REF] to impose the contact condition is one of the form of solving this contact problem, however different methods can be constructed starting from this. In this context, different FETI-like type of methods have been developed to solve this non-linear problem [4], [START_REF] Dureisseix | A numerically scalable domain decomposition method for the solution of frictionless contact problems[END_REF], [START_REF] Dostál | Scalable FETI Algorithms for Frictionless Contact Problems[END_REF]. They are all based in the basic linear FETI method.

To state our comparative test case, we will use some results of a particular contact solver thanks to Alexandros Markopoulos, [START_REF] Dostál | Scalable TFETI algorithm for the solution of multibody contact problems of elasticity[END_REF]. So from now on, we will only consider a FETI-like with the penalty method to solve the contact problem, since there is a better treatment of the contact constraint with no extra variables at the cost of ill-conditioned problems, that we will try to overcome.

For a contact problem, after the solver convergence we obtain a final active area of contact, with this we can "lock" the substructures in this final configuration, just by considering both contact bodies as a single conforming structure. This characteristic comes from the fact that a node-to-segment discretization is done in the contact area, that also treats with the large displacements nonlinearities of the model [START_REF] Wriggers | Discretization, Large Deformation Contact[END_REF]. This type of discretization, along with the penalty method, implies that, virtual contact elements are added on one of the contact bodies, this elements from a mechanical point of view represent the elements that connect this body to an imaginary string, that when the contact is active, should have a very high stiffness, in order to impose the non penetration condi-tion.

When iterating to solve the non linearities, at the beginning the penalty parameter is not necessarily very high and non-physical penetration may occur, but when approaching to the final computed loads we must increase the value of it, in order to nullify this penetration and have a correct imposition of the contact condition. Is this treatment, the one that produces the very ill-conditioning of the local stiffness matrices, when the contact area is active.

For this reason that, unlike FETI-1LM, the FETI-2LM is a more suited method to solve this kind of problem. The difference is that in FETI-2LM the operator is provided with augmentation or regularization matrices, added to this local stiffness problems, thanks to the generalized Robin condition imposed in the interfaces. This matrices are usually taken as approximations of the Schur complement of the neighbour subdomain, meaning that even if the stiffness of one neighbour subdomain is orders of magnitude larger than the other (this happens when the penalization parameter is acting in contact) it will be taken into account by its neighbour.

To see this regularization, we just need to recall the definition of the FETI-2LM operator for a two subdomain case

       I I -(A (1) b + A (2) b )(S (2) bb + A (2) b ) -1 I -(A (2) b + A (1) b )(S (1) bb + A (1) b ) -1 I        (1.153)
where we take A b ≈ S 1 bb in order to have an operator less affected by the ill-conditioning produced by the penalization method.

After the computation of the final contact area, we can perform a standard stress analysis of that configuration under the computed external loads produced by the interaction of both structures, including the penalization in the local stiffness matrices that share an interface with active contact. In this static problem we can compare the two standard FETI method against the new Hybrid solver.

In the first problem to solve, we have a block in top of a semi-circle, see Figure 1.15. The top block is the stiffer body and the bottom one is assumed to

E top = 10000, ν top = 0.4 E bottom = 6, ν bottom = 0.3 (1.154)
The total domain is divided into six subdomains, two for the bottom structure and four for the top one, then a vertical gravity-like force is applied, i.e f = (0, -10). The Dirichlet condition is applied to the base of the bottom structure, and in order to give an initial stability, the two top corners of the block are also fixed.

Then, in Figure 1.16 we can see the solution, after applying the contact solver. This allows to compute the final active contact zone, that is the bottom of the square block, in which the ill-conditioning is concentrated. The forces needed to obtain such a bending are also given by the solver.

One of the main issues shown in the previous test case was the definition of the 1LM or 2LM interfaces. When solving contact problems, we have the advantage of working with two separate and, a priori, known structures. This allows to recognize each subdomain as part of one or the other structure, which at the same time automatizes the process of marking the different interfaces in the following way

       Mark Γ (ij) as 1LM if Ω (i) , Ω (j) are in the same body Mark Γ (ij) as 2LM Otherwise (1.155)
With this process, plus the Robin or Dirichlet preconditioner, diminishing the stopping criterions to errors less than 10 -6 and a patch size of 3 and path depth of 5 for the augmentation matrices in the 2LM interfaces, we can launch our test for the final configuration of the contact problem, see results in Table 1.4.

The following final two cases of this chapter, model a physically more stable structure, where the same top block is now over an square bottom structure. This time we have augmented the number of subdomains to emphasize the advantages of the Hybrid method.

The top and first structure is divided into 16 subdomains and the bottom and second into 9 subdomains for the first case and 25 subdomains for the second. Dirichlet conditions are imposed in the base of the second structure and the top block is completely free, as the contact zone can offer enough initial stability, see Figure 1.17.

Again the contact solver gives us the final contact interfaces plus the forces In the results we can see that for both cases, the number of iterations for the Hybrid method is smaller than for the other two methods. Compared to the FETI-1LM we have gained in robustness since we only see convergence in the much more stable structure that is the second case, but also we improved the convergence ratio in this case, showing a better general performance for this type of problems. If we compare against the FETI-2LM method, we can see that in the biggest case the number of iterations used by the Hybrid method is less than half of the ones needed originally and at the same time we manage to keep the robustness given by this method.

Finally we point out that when the number of subdomains augments, we only see a small increase in the total iterations, in accord with the results shown in the first example of the two material bar.

Conclusion

In general the results given by the Hybrid-FETI method show us a positive first approach for the development of a general FETI method that tries to keep the good performance of the basic FETI-1LM method and the robustness of FETI-2LM. This advantages are clearer in the context of solving the ill-conditioned problems arising from the simplest contact models. In this cases the advantages in terms of easier implementation and improved convergence shown in the numerical results gives us a very good start in the search of general better algorithms.

The next steps for this method is to perform test in bigger and less "academic" test, but since our preliminary results show and improvement, we can only hope the same for the bigger and more complex cases. Also, a different type of problems can be tested, which are a class of non-conforming problems where the use of mortar elements can affect the conditioning of the local matrices forcing the use of FETI-2LM in the complete structure [START_REF] Roux | A FETI-2LM Method for Non-Matching Grids[END_REF], where in practice only a small part of the problem needs such a treatment, making or method also appropriate to this type of problems.

Chapter 2

Block FETI methods Different improvements to make more robust and efficient FETI methods have been developed in recent years, but for some harder engineering problems the performance is still very poor. To face this challenging problems, the Simultaneous-FETI or S-FETI was developed, first in [START_REF] Rixen | Substructuring and dual methods in structural analysis[END_REF] for a two subdomain case, then generalized for arbitrary configurations in [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF]. The main characteristic of the S-FETI method it's the generation from the preconditioning of more than one search direction for the Conjugate Gradient method, leading to a more efficient and in general more robust method.

With this new method, as usual in the practical use of FETI methods, a full reorthogonalization is mandatory, forcing the storage of all the new search directions. Depending on the problem, specially when the number of subdomains increases, the total number of directions to store can be a drawback due to memory limitations, making the S-FETI method impractical for this type of problems.

In this chapter we treat the problem of memory limitation of the S-FETI method using the sparse properties of the block of vectors from which we construct the search directions, in order to reduce the total memory allocated by the method. The strategy will be to rebuild at every iteration the search directions using this sparse blocks and some small coefficient matrices.

The introduction of this reconstruction steps will add computations that from a local point of view will be more expensive, but since this cost is small 81 compared to the time spend in the communication of the different processes, the total time used can be as good as the original method. The simplest of the implementations of this idea will not be enough to achieve the good performance needed, but it will show the pertinence of the storage used and the precision achieved in relation to the method with no special storage. Later we will show some optimizations based on the parallel properties of the code that improve the computation time, making a more efficient algorithm. This chapter start showing the basics of preconditioning in FETI, then we show how the ideas in the preconditioning steps lead to the formulation of the basic S-FETI method for two subdomains. We continue by presenting the formal general method, including its implementation and cost. Then we show in detail the way of sorting search directions in the presence of linear dependency between them, we also add a second strategy to achieve this. We continue with the new storage, reconstruction of the search directions and the optimizations in the code to finally show some numerical results to compare the different algorithms.

Introduction and preliminarities

The development of the S-FETI method begins in [START_REF] Rixen | Substructuring and dual methods in structural analysis[END_REF] when searching for improvements to the Dirichlet preconditioner in FETI. For this reason, in this section, we will follow the same path from the basics of the preconditioner that will lead to the first version of S-FETI.

Dirichlet preconditioner for two subdomains

We will start by revisiting the Dirichlet preconditioner at its basic, for a two subdomain case. The notation in this case will be analogous to the previous chapter.

The objective of a preconditioner for the PCPG Algorithm 5 used in FETI, is to build in every iteration p a correction to the Lagrange multipliers given an interface compatibility error, i.e. a jump denoted as g. (In the case of the elasticity problem this is a jump in the displacement, but depending on the problem this may change. Also we drop the notation as an iterative vector g p for simplicity).

Let us consider the two subdomain case as in Figure 1.2 from section 1.1.1.

Where the associated unknown are

x (s) =        x (s) i x (s) b        , s = 1, 2 (2.1)
the subscripts i, b represent interior and interface nodes respectively. The jump is then written as g = B (1) x

(1) b + B (2) x (2) b (2.2)
with B the signed boolean operator that defined this continuity condition as in Equation 1.33. We can then build a continuous interface approximate solution by averaging the already computed x (s)

x(1

) b = 1 2 x (1) b -B (1) T B (2) x (2) b x(2) b = 1 2 x (2)
b -B (2) T B (1) x

(1) b

(2.

3)

The operator B (s) T B (q) is the correspondence between the numbering of the interface nodes in Ω (s) and Ω (q) . The interface approximate solution corrections are

δx (1) b = x(1) b -x (1) b = - 1 2 B (1) T g δx (2) b = x(2) b -x (2) b = - 1 2 B (2) T g (2.4)
If we modify the internal nodes to satisfy the internal equilibrium (the first equation in the local problem 1.29) then we have the correction so the interface response is not uniquely defined. So again, an averaging is done as previously with the approximate solution, then the interface response correction will be chosen as

δx (s) i = -K (s) -1 ii K (s) ib δx
z = - 1 2 B (1) δf (1) b + B (2) δf (2) b (2.7)
This correspond to a new search direction in the CG algorithm. If we look at this definition, using 2.4 and 2.6, we can also write it as

z =        1 4 2 s=1 B (s) S (s) bb B (s) T        g (2.8)
which is equivalent to the definition of the Dirichlet preconditioner of subsection 1.1.3 with the scaling matrices in consideration. If the two matrices S (s) bb are the same, for example in a totally symmetric splitting of a simple square, then the preconditioner is the exact inverse of the operator FETI, if not the case, we only have locally the exact inverse of this operator.

The difference between this Dirichlet preconditioner and his basic non-scaled form lies in what its defined as a Consistent preconditioner. In the next part, we will give the definition and develop this subject in order to generalize the preconditioner to arbitrary partitions and heterogeneous problems.

Consistent preconditioners Preliminaries

We begin by recalling the signed boolean matrices needed to formulate the problem in general subdomain division. For an arbitrary partition into N s ≥ 2, we define

Γ (s) = ∂Ω (s) \∂Ω (2.9)
this is the interface boundary of Ω (s) for s = 1, . . . , N s , then we define the global interface as

Γ = ∪ s Γ (s) (2.10)
Next we consider the signed boolean matrix B (s) as the matrix that maps the nodes from the local interface Γ (s) into the global one Γ . The sign of the represented nodes is such that the opposite sign is in the position that represents the same global node on the neighbour subdomain.

The global interface will again be divided into interface edges Γ j defined as

Γ j := Γ (sq) = ∂Ω (s) ∩ ∂Ω (q) , ∀s, q = 1, . . . , N s (2.11)
The crosspoints are then nodes shared by more than two edges. It follows that B (s) : Γ (s) → Γ can be partitioned as

B (s) =                    B (s) Γ 1 B (s) Γ 2 . . . B (s) Γ n (s)                    (2.12)
with n (s) the number of neighbours of Ω (s) and where B (s) Γ j is the restriction of B (s) to Γ j . We define also the global assembly operator

B = B (1) • • • B (N s ) (2.13)
Next we define the multiplicity of a node in Γ j as m j , and for each node its values is

m j = |neighbours| + 1 (2.14)
the multiplicity varies in each node of the edge, but for simplicity we will denote as a single one for each Γ j . In general, m j ≤ 2 and for a crosspoint m j > 2. Since one lagrange multiplier is used to glue any pair of d.o.f in an edge Γ j , there are exactly (m j -1) lagrange multipliers applied to each d.o.f in any edge.

From equation 2.12 and knowing that there are nodes shared by more than two subdomains we can see that

B (s) T Γ j B (s) Γ j = (m j -1)I (2.15)
then for each subdomain we have that

B (s) T B (s) = B (s) T Γ 1 . . . B (s) T Γ n (s)              B (s) Γ 1 . . . B (s) Γ n (s)              = B (s) T Γ 1 B (s) Γ 1 + • • • + B (s) T Γ n (s) B (s) Γ n (s) = (m 1 -1)I + • • • + (m n (s) -1)I = diag(m j -1) (2.16)
Finally we can write

B (s) T B (s) + I = diag(m j ) (2.17)
where diag(m j ) is the diagonal matrix with the multiplicity of the nodes in the edges Γ j ⊆ Γ (s) .

Remark: This form of gluing connecting d.o.f introduces redundancies in the compatibility constraints at the crosspoints (continuity). However we will see that this redundancy is essential for an efficient preconditioner.

Consistent preconditioner

With previous definitions, we can analyse the preconditioner in terms of a general partition of the domain in N s > 2 subdomains.

Recalling the physical interpretation of the projected gradient in the FETI-1LM method, we know it is a jump of the displacement field across the subdomain interface boundaries

g = N s s=1 B (s) x (s) b (2.18)
So, from a mechanical point of view the objective of a preconditioner F-1 based on local problems is to generate a correction of the Lagrange multipliers z and its corresponding local interface forces B (s) T z = B (s) T [ F-1 g] in order to reduce the jump g as much as possible.

The basic preconditioning, either Dirichlet or Lumped, applied to the projected gradient g is defined in general terms as

z = F-1 g = N s s=1 B (s) S (s) bb or K (s) bb B (s) T g (2.19)
this operator is build in a three-step procedure, starting with g 1 We define the approximate solution corrections δx The first condition is fulfilled if we use an averaging process as the one shown for the two subdomain case. This idea of imposing a common value (e.g. the average) to every node independently on the interface, will also be valid for the general case, i.e. considering crosspoint.

Graphically, the three steps construction can be seen in Figure 2 b . Nevertheless, we require that z be constructed in such a way that as the δf

(s) b approach equilibrium, B (s) T z is exactly δf (s) b .
Lets consider the global interface response correction δf b as the vector defined in the interface, such that δf

(s) b = δf b | Γ (s) (2.24)
We can say now that a Consistent Preconditioner is the one that it has 1. Consistent approximate solution increments δx

(s) b , i.e. N s s=1 B (s) x(s) b = N s s=1 B (s) x (s) b + δx (s) b = 0 (2.25) 2. Consistent Lagrange multiplier corrections z, i.e. if δf b ∈ Im(B) (2.26) then, for each s = 1, . . . , N s B (s) T z = δf (s) b (2.27)
The Equation 2.26 is the interface equilibrium as it express that the sum of this boundary interactions δf At a subdomain level the condition 2.26 can also be written

δf (s) b -B (s) T N s r=1 r s B (r) δf (r) b = 0 (2.29)
With this definition, it can be expected that preconditioners that fulfills it are better than the ones who does not, this was shown in [START_REF] Daniel | A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems[END_REF]. In this sense we will explain the modifications done to the basic Dirichlet preconditioner in order to make it a consistent one, leading to a new structure that will be generalized in order to have a preconditioner also suited for heterogeneous problems (across the interface).

Remark:

1 The basic Dirichlet or lumped preconditioner are not consistent. If we look at the definition of δx (s)

b we have that

N s s=1 B (s) x (s) b + δx (s) b = g + N s s=1 B (s) B (s) T g 0 (2.30)
which violates the compatibility condition for the solution correction.

Furthermore, if we take the construction

z = N s s=1 B (s) δf (s) b (2.31)
and we assume that Equation 2.29 holds, then we have

B (s) T z = B (s) T N s r=1 B (r) δf (r) b = B (s) T B (s) δf (s) b + B (s) T N s r=1 r s B (r) δf (r) b = B (s) T B (s) + I δf (s) b (2.32)
and using Equation 2.17

B (s) T z = diag(m j )δf (s) b δf (s) b (2.33) 2
The natural self-equilibrium of the Lagrange multipliers can be expressed replacing δf

(s)
b by B (s) T z in Equation 2.29

B (s) T z -B (s) T N s r=1 r s B (r) B (r) T z = 0 (2.34) Ω 2 Ω 4 Ω 1 Ω 3 x (2) 2 g x (1) 2 x (1) 1 x (2) 1 x (4) 1 x (3) 1 x Figure 2.2 -4 Subdomain problem
from where we have

B (s) T = B (s) T N s r=1 r s B (r) B (r) T (2.35)
this is a direct consequence of the fact that the compatibility conditions are defined between any pair of connecting d.o.f.

Consistent preconditioner for homogeneous problems

For homogeneous problems, consistent preconditioners are constructed by extending the averaging scheme used in the two subdomain case, see subsection 2.1.1. For clarity we will start showing it for a 4 subdomain case as in Figure 2.2 then generalizing it for arbitrary meshes.

We assume that all subdomains have similar stiffness matrices, i.e. are build for subdomains of similar materials, geometrical and discretization properties. Therefore, the construction of the compatible interface solution is done by imposing an averaging as follows

x(1) 1 = x(2) 1 = x(3) 1 = x(4) 1 = x (1) 1 + x (2) 1 + x (3) 1 + x (4) 1 4 x(1) 2 = x(2) 2 = x (1) 2 + x (2) 2 2 (2.36)
which implies that the following corrections are consistent

δx (1) 1 = x(1) 1 -x (1) 1 = x (2) 1 -x (1) 1 + x (3) 1 -x (1) 1 + x (4) 1 -x (1) 1 4 δx (2) 1 = x (1) 1 -x (2) 1 + x (3) 1 -x (2) 1 + x (4) 1 -x (2) 1 4 . . . δx (1) 2 = x (2) 2 -x (1) 2 2 δx (2) 2 = x (1) 2 -x (2) 2 2 (2.37)
If the interface nodal responses are computed as in Equation 2.21, then, in general, they will not be in equilibrium. For instance if we see the node 2, we have that δf

(1) 2 δf (2)
2 . So the Lagrange multiplier correction is constructed by doing the same averaging between the results of the Dirichlet-to-Neumann mapping, as follows

z 7 = -δf (1) 2 + δf (2) 2 2
(2.38)

The minus sign in here comes from the fact that this vectors have opposite directions, because by convention, the entries of B (s)

Γ j (Γ j connecting subdomains s and q) are +1 if s > q and -1 otherwise. We chose to do same averaging to keep a symmetric preconditioner. Doing the analogous averaging for the crosspoint of multiplicity 4, we have the corrections

z 1 = -δf (1) 1 + δf (2) 1 4 , z 2 = -δf (2) 1 + δf (3) 1 4 z 3 = -δf (3) 1 + δf (4) 1 4 , z 4 = -δf (1)
1 + δf (4) 1

4

z 5 = -δf (1)
1 + δf

(3) 1 4 , z 6 = -δf (2)
1 + δf

(4) 1 4
(2.39)

and the evaluation of z = F-1 g is complete.

We can check that this Lagrange multiplier correction is consistent, in fact if then

B (1) T z =       -z 1 -z 4 -z 5 -z 7       =           1 4 δf (1) 1 -δf (2) 1 + 1 4 δf (1) 1 -δf (4) 1 + 1 4 δf (1) 1 -δf (3) 1 δf (1) 1 -δf (1) 1 2           =        δf (1) 1 δf (1) 2        (2.41)
And the same happens in analogous way with the other subdomains, so for the 4 subdomain case the averaging process shows the construction of a consistent preconditioner (Dirichlet or Lumped).

Lets show the extension of this averaging process to arbitrary meshes. Consider any edge Γ i , then the compatible solutions x(s) Γ i can be constructed by averaging all the nodes belonging to Γ i x(s)

Γ j = r:Γ i ⊆{Γ ∩Γ (r) } 1 m i x (r) Γ i (2.42)
where m i is the multiplicity of nodes in the edge Γ i . Then we build the following consistent corrections δx (s)

Γ j = x(s) Γ j -x (s) Γ j = r:Γ i ⊆{Γ ∩Γ (r) } r s 1 m i x (r) Γ i -x (s) Γ i (2.43)
and extended to the total subdomain interface we have

δx (s) b = -B (s) T diag 1 m j N s r=1 B (r) x (r) b = -B (s) T E (s) g (2.44)
where E (s) is a diagonal matrix with values that correspond to the multiplicity of nodes in the edges that intersect Ω (s) . This matrix can be seen as a scaling matrix.

Now we compute the interface responses δf (s)

b associated with this corrections δx (s) b , the same as before we simply use Equation 2.21, as for the consistent Lagrange multiplier corrections, we generalize the averaging showed in 2.38 and 2.39 to obtain

z = F-1 g = diag 1 m i N s r=1 B (r) δf (r) b (2.45)
Again, lets check the consistency of this correction. We assume the interface equilibrium condition 2.26 and using also 2.17 we have

B (s) T z = diag 1 m i N s r=1 B (s) T B (r) δf (r) b = diag 1 m i B (s) T B (s) δf (s) b + diag 1 m i N s r=1 r s B (s) T B (r) δf (r) b = diag 1 m i (diag (m i ) -I) δf (s) b + diag 1 m i δf (s) b = δf (s) b (2.46)
from which we have the consistency of the Lagrange multiplier correction z.

With this, we can write the consistent version of the Dirichlet and Lumped preconditioners

D -1 = s E (s) B (s) S (s) bb B (s) T E (s) (2.47) L -1 = s E (s) B (s) K (s) bb B (s) T E (s) (2.48)

Consistent preconditioner for heterogeneous problems

The previous scaling works in homogeneous problems when the stiffness matrices of neighbouring subdomains are similar, and we can think that the compatible field is in the "middle" of them. In heterogeneous problem this is not the case, as we can presume that the compatible solution will be closer to the subdomain with a higher stiffness (from a mechanical point of view), making the previous scaling to degrade in terms of convergence for this kind of problems.

To correct this behaviour, we will define a more general scaling, also consistent, that can acknowledge this differences.

We start again by explaining the idea in a 4 subdomain case, as in Figure 2.2. In this configuration we define the more general compatible field as

x(1) 1 = x(2) 1 = x(3) 1 = x(4) 1 = β (1) 1 x(1) 1 + β (2) 1 x(2) 1 + β (3) 1 x(3) 1 + β (4) 1 x(4) 1 x(1) 2 = x(2) 2 = β (1) 2 x(1) 2 + β (2) 2 x (2) 2 (2.49)
where the β (s) k terms are smoothing or weighting coefficients defined in every interface d.o.f. k of subdomain Ω (s) . This coefficients must also be constrained by

β (1) 1 + β (2) 1 + β (3) 1 + β (4) 1 = 1 β (1) 2 + β (2) 2 = 1 (2.50)
with this, the corrections vanish when the solution before weighting is already then

B (1) T z =       -z 1 -z 4 -z 5 -z 7       =         β (2) 1 + β (4) 1 + β (3) 1 δf (1) 1 -β (1)
1 δf

(2)

1 + δf (3) 1 + δf (4) 1 β (2) 2 δf (1) 2 -(1 -β (2) 2 )δf (2) 2         =        δf (1) 1 δf (1) 2        (2.54)
The same happens for the other subdomains.

This process defines a generalization of the weighting procedure for an arbitrary number of subdomains. The analysis is analogous as previous, but considering a more general scaling matrix β (s) , defined as the diagonal matrix of weighting coefficients of the interface d.o.f. belonging to the neighbours of Ω (s) .

The choice of this coefficients can be done using a physical criteria. In the case of an elasticity problem (works for any elliptic PDE problem) we decouple all the interface d.o.f. assuming that each one of them is connected to a stiffnessfree covering subdomain via a spring, in other words we see the stiffness of each subdomain lumped to its interface, so the lumped stiffness matrix is diagonal equal to diag K (s) bb , this weighting procedure is usually called Superlumped.

In general the coefficients are computed as the ratio between the stiffness of an interface d.o.f. in some subdomain and the sum of all the stiffness of the connected to this d.o.f. in all neighbouring subdomains. For example, for the 4 subdomain case, we have

Γ i = diag K (s) Γ i          r:Γ j ⊆{Γ ∩Γ (r) } diag K (s) Γ j          -1 (2.56)
So finally the preconditioners Dirichlet and Lumped respectively, added this Superlumped scaling can be written as • Construction of diagonal matrices β (s) .

D -1 = s β (s) B (s) S (s) bb B (s) T β (s) = s B(s) S (s) bb B(s) T (2.57) L -1 = s β (s) B (s) K (s) bb B (s) T β (s) = s B(s) K (s) bb B(
• At every multiplication by the preconditioner we add two extra by this local diagonal scaling, which are negligible in parallel computation.

Simultaneous FETI

As seen in previous subsection the Dirichlet (or Lumped) preconditioner are two correctors of the Lagrange multipliers imposed in the interface, they use an averaging or scaling on the interface to make this correction consistent. This averaging is necessary because the interface computed responses δf The idea is now to use this two different corrected responses independently, as search directions for the Conjugate Gradient method, as shown in [START_REF] Rixen | Substructuring and dual methods in structural analysis[END_REF].

We start again by showing this in a two-subdomain division, so lets consider the p iteration of the CG method applied to the FETI-1LM method (no projection is consider in this simple case). Then, we first solve the local Dirichlet problem consisting in enforcing on the interface the compatibility gap g p and we compute the boundary reaction which write

δf (s),p = S (s) bb B(s) T g p , s = 1, 2 (2.60)
In the classical Dirichlet preconditioner we evaluate and approximate correction to the Lagrange multipliers by computing a weighted inter-subdomain vector based on δf (s),p . Let us now consider this boundary conditions on each side of the interface as descent direction for the CG algorithm by writing the new update of the interface Lagrange multiplier λ λ p+1 = λ p + ρ (1),p B(1) δf (1),p + ρ (2),p B(2) δf (2),p (2.61)

Assuming now that, instead of using this boundary reactions directly, we use δλ (1),p = B(1) δf (1),p +

2 s=1 p-1 l=0 α (s),l δλ (s),l δλ (2),p = B(2) δf (2),p + α (2),p δf (1),p + 2 s=1 p-1 l=0 γ (s),l δλ (s),l (2.62)
such that δλ (2),p is orthogonal to δλ (1),p and that both directions are orthogonal to all previous directions δλ (1),p T Fδλ (2),p = 0 δλ (s),p T Fδλ (r),l = 0, r, s = 1, 2 l = 0, . . . , p -

Then, the new update is

λ p+1 = λ p + ρ (1)
,p δλ (1),p + ρ (2),p δλ (2),p g p+1 = g p + ρ (1),p Fδλ (1),p + ρ (2),p Fδλ (2),p (2.64) the direction coefficient ρ (s),p can be determined by this new orthogonality relations (g p+1 • δλ (1),p ) = (g p • δλ (1),p ) + ρ (1),p (Fδλ (1),p • δλ (1),p ) + ρ (2),p (Fδλ (2),p • δλ (1),p ) 0 = (g p • δλ (1),p ) + ρ (1),p (Fδλ (1),p • δλ (1),p )

(2.65) (g p+1 • δλ (2),p ) = (g p • δλ (2),p ) + ρ (1),p (Fδλ (1),p • δλ (2),p ) + ρ (2),p (Fδλ (2),p • δλ (2),p ) 0 = (g p • δλ (2),p ) + ρ (2),p (Fδλ (2),p • δλ (2),p ) (2.66) which implies

ρ (s),p = - (g p • δλ (s),p ) (δλ (s),p • Fδλ (s),p ) , s = 1, 2 (2.67)
that correspond to the same as the regular CG method. With this we can note that one of the fundamental features of this new algorithm is that even if at each iteration the minimization is done with respect to two descent directions, the cost of this is equivalent to the cost of a normal conjugate gradient direction. Actually, recalling that the directions come from a Dirichlet problem, we note that F B(1) δf (1),p = B(1) S (2.68)

F B(2) δf (2),p = B(1) S (1) -1 bb B(1) T + B(2) S (2) -1 bb B(2) T B(2) S (2) bb B(2) T g p = g p + B(1) S (1) -1 bb B(1) T B(2) δf (2),p (2.69)
And in this case we see that applying the FETI operator to both descent directions only requires one Neumann solution per subdomain.

In the generalized version, we don't expect to have this exact same property of cost equivalency between method, but we do expect a reduced number of multiplication by the FETI operator that will lead to a faster method in general.

We will now show the extension done in [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF] to this new method for arbitrary meshes.

The Algorithm

The idea is to exploit the additive structure of the preconditioner in FETI and generate several search directions, instead of one, in every step of the CG method. The basic S-FETI method generates one for each subdomain, but in a straightforward way this number can be two times the total number of local interfaces in the problem.

In this section the notation and definitions comes from the previous chapter, where the original FETI method is described. With this in mind, we can remember the construction of a search direction in the classical FETI method. We denote this direction w ∈ R n , n being the size of the interface, used by the Conjugate Gradient algorithm.

The consistent Dirichlet operator, defined in 2.57, is first applied to the global residual vector g

w = D -1 g =       s D (s) -1       g (2.70)
where

D (s) -1 = B(s) S (s) bb B(s) T (2.71)
then the vector w is orthogonalized with respect to previous search directions to generate the new one.

The idea of S-FETI is to improve the minimisation process in CG by spanning a search space, not from the addition of local contributions, but from each of this terms separately, letting the process to choose the best combination. Although more costly, this optimal combination minimizes the residual in the space generated by

Z = D (1) -1 g, D (2) -1 g, . . . , D (N s ) -1 g (2.72)
where N s is the number of subdomains.

Each one of this columns is then projected and orthogonalized to give the block W of N s columns, where each column correspond to a new search direction. The update of the solution λ of the CG algorithm at the p iteration is done simultaneously (analogous to Equation 2.64) for all this orthogonal vectors in a similar way that the classic CG, by adding the linear combination of these directions that minimizes the error in the F-norm λ p+1 = λ p + W p ρ p (2.73) with ρ p ∈ R N s such that

ρ p = (W T p FW p ) + Z T p g p (2.74)
Since the classical search direction is such that

w p = W p (1, . . . , 1) T ∈ range(W p ) (2.75)
this new approximation is better than the usual one given by the CG algorithm.

We can see the need of find the inverse of the matrix W T p FW p ∈ R N s ×N s . Because W p is formed from local contributions, we would expect that it is fullranked. This is not always the case (details will be given in later sections) and we only have a symmetric positive semidefinite matrix that can be pseudo-inverted.

To avoid this pseudo inversion, another equivalent option is to eliminate some directions to recover a smaller full-ranked family in W p . The approximated solution will be the same, but fewer vectors would need to be stored at every iteration.

To build this smaller full-ranked family, a rank revealing Cholesky factorization with complete pivoting is used, giving

N T (W T p FW p )N = LL T (2.76)
where N is a permutation matrix and

L =       L 0 × 0       (2.77)
with L a lower triangular matrix of full rank. The F-orthogonalization of the directions in the block W p can be done by

W p ← W p N       L-T 0       FW p ← FW p N       L-T 0       (2.78)
where we are suppressing the redundant directions. The computation of the optimization parameters ρ p is now

ρ p = W T p g p (2.79)
In Algorithm 8 we can find the description for this method. The definition of the projection and the matrices to compute λ 0 are defined as in previous FETI methods [START_REF] Axelsson | A Black Box Generalized Conjugate Gradient Solver with Inner Iterations and Variable-Step Preconditioning[END_REF]. for i = 0 to p do 15:

Algorithm 8 S-FETI algorithm 1: Initialization 2: λ 0 = AG[G T AG] -1 (-R T c
Φ i = -W T i FW p+1
16:

W p+1 = W p+1 + W i Φ i 17: FW p+1 = FW p+1 + FW i Φ i 18:
end for 19: end loop We note that when building the block of search directions in every iteration we are not making this block orthogonal to the ones computed in previous steps, meaning that in fact we are losing the short recurrence of the Conjugate Gradient method, adding the need of doing a full orthogonalization at every iteration. In any case this drawback is only theoretical since in practice this full orthogonalization is also needed in order to keep the numerical accuracy of FETI methods in general.

With this orthogonalization it can be proved the following minimization property

Theorem 2.1. The approximate solution computed by the p iteration of the S-FETI method minimizes the error λ p -P λ F over the space

λ p ∈ p-1 i=0 span W i (2.80)
where ⊕ is the direct sum and W i is defined by Algorithm 8

Proof. The proof is done in [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF] following the usual demonstration for CG, see [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF].

One particularity of this method is that we are no longer minimizing over a Krylov space. This is because at every iteration, the approximate solution is updated in the different directions given by the optimal combination of all local preconditioners, making the coefficients ρ to change from one iteration to another. This fact won't allow to find a bound (heuristic) for the number of iterations, but in any case we will expect a good robustness that can be explained by the similarities of this method and the FETI-Geneo method [START_REF] Spillane | Automatic spectral coarse spaces for robust finite element tearing and interconnecting and balanced domain decomposition algorithms[END_REF].

Remark The S-FETI method has been shown to be very efficient on hard problems where the classical FETI require many iterations. The convergence is comparable to the one in FETI-Geneo algorithm where a coarse space is constructed by solving in each subdomain and at every iteration the generalized eigenvalue problem

S (s) bb v (s) -µB (s) T D -1 B (s) v (s) = 0 s = 1, . . . , N s (2.81)
This problem allows to isolate the part of the solution on which the preconditioner is not sufficiently efficient for the iterative solver to perform well. The vectors detected are the ones where the restriction of the global preconditioner is not a good approximation of the non-assembled local component S (s) bb of the FETI operator. In S-FETI the solution space comes from the successive applications of the local non-assembled components B (s) S (s) bb B (s) T and the assembled FETI operator F, so the block of search directions spans a space where the local effects are taken in account. It is then similar to the deflated space where the Geneo iterations take place, and thus convergence is expected to be very fast.

Extension to local interface division

In the regular first version of the S-FETI method, the preconditioner is decomposed in the local contributions

D (s) -1 = B (s) S (s) bb B (s) T (2.82)
where S (s) bb is the Schur complement of the internal nodes into the interface ones. This implies that we use the complete local interface Γ (s) to build the different search directions, but different interactions can occur for the same subdomain, depending on the characteristics of each of his neighbours. If all of the neighbours have similar properties, then the use of the complete local interface is justified. If this is not the case, for example in Figure 2.3, then we can again divide the local interface into the different interface edges, one for each neighbour, to take into account this differences, increasing the number of search directions built and improving the approximation of the local behaviour, hence ameliorating the convergence of the method.

The formal construction of this search directions comes from the definition of the preconditioner, and the definition of the interface edge given in previous

E 1 =0,1 E 3 =1000 E 2 =10
Γ (1) 1 then in the S-FETI method we consider the local contributions separately

Γ (2) 1 Γ (2) 2 Γ (1) 2 Γ (3) 1 Γ (3
D (s) -1 = B(s) S (s) bb B(s) T (2.84)
from which we created the N s different search directions. This operator is applied in the CG method to the computed gradient g, defined in the total interface Γ .

To do this we first multiply by the assembling scaled matrix B(s) T and then we compute the forces needed to have this displacement

δf (s) bb = S (s) bb B(s) T g (2.85)
We know that this forces are defined in the local interface Γ (s) , but since their interactions may change from one neighbour to another, we create new vectors that reflect this interactions separately by considering the restrictions to the interface edges Γ (sq) = ∂Ω (s) ∪ ∂Ω (q) δf (s)

Γ (sq) = δf (s) bb | Γ (sq) s = 1, . . . , N s , q = 1, . . . , n (s) (2.86)
where n (s) is the number of neighbour subdomains in Ω (s) . We extend this vector by zero to match the size of the interface Γ (s) f (s)

Γ (sq) =                                    0 . . . 0 δf (s) Γ (sq) 0 . . . 0                                    (2.87)
each one of this vectors is rescaled to recover the consistent property, and with this form the columns of the block used to build the different search directions

D (s) -1 Γ (sq) g = B(s) f (s) Γ (sq)
s = 1, . . . , N s , q = 1, . . . , n (s) (2.88) so finally the search space will now be generated by the block

Z = . . . , D (s) -1
Γ (sq) g, . . . s = 1, . . . , N s , q = 1, . . . , n (s) (2.89) in this case we won't be having one column from each subdomain, but instead there will be one for each local edge (or neighbour). This will allow to span an even bigger search space at the cost of computing all the extra directions.

Graphically we can see the difference between this block and the regular one,

Ω 1 Ω 3 Ω 2 z 1 z 2 z 3 z 4
Subdomain division ⇒ Z =       z 1 z 2 0 0 z 3 z 4       Interface division ⇒ Z =       z 1 z 2 0 0 0 0 z 3 z 4       (2.90)
We note here that when the problem is divided in much more subdomains, the block Z will show an sparse pattern, since in every column, only the values associated to Γ (s) (or in particular Γ (sq) ) are non zero. This fact will be a very important in the practical implementation of the method, as we will see later in this section.

As done in the regular S-FETI to compute the final search directions, each of the columns of Z needs to be projected and then orthogonalized with the previous directions. Then we use the same rank revealing strategy 2.78 to compute the inverse of W T p FW p and eliminate the useless directions. The algorithm is analogous to the previous one and is described in 9.

Both algorithm for S-FETI, shows an implementation in general lines of the method, but several optimizations can be done by using different strategies, all of this in order to reach for the maximum performance of this method.

Algorithm 9 S-FETI algorithm interface

1:

Initialization 2: λ 0 = AG[G T AG] -1 (-R T c) 3: g 0 = P T (Fλ 0 -d) 4: Z 0 = [. . . , D (s) -1
Γ (sq) g 0 , . . . ], s = 1, N s q = 1, . . . , n (s) 5: W 0 = P Z 0 6: loop Iterate p = 0, 1, 2, ... until convergence Γ (sq) g p+1 , .

. . ], s = 1, N s q = 1, . . . , n (s) 13:

W p+1 = P Z p+1 14:

for i = 0 to p do 15:

Φ i = -W T i FW p+1
16:

W p+1 = W p+1 + W i Φ i 17: FW p+1 = FW p+1 + FW i Φ i 18:
end for 19: end loop

Cost and implementation of S-FETI

One of the characteristic of this method is the reduction in the number of total iterations versus the classical FETI, at the cost of the extra computations in each step. However, even if the cost of in each iteration does require a much larger computational effort, several strategies in the implementation of the code leads to efficient local and global optimizations, in order to reduce this cost differences. We start by summarizing the main changes in the cost of S-FETI vs FETI-1LM.

1. The number of exchanges phases is the same. The communications between neighbours are almost identical, with one less interface exchange when computing the Dirichlet preconditioner.

The main difference is that the exchanges in S-FETI involve more data, increasing to N T × N T for the matrices W T p FW p , Φ and to N T for the vector ρ, where N T is either the number of subdomains (in the basic S-FETI) or N T = O(N s ) (actually is two times the number of local interfaces) when using the local interface subdivision. Also the full reorthogonalization process is now done in a block way, so we change the vector operations by matrix operations, so again, same number of exchanges and computation, but with more data involved.

The advantages from keeping the number of interchanges, but increasing the size of the data are more clear if we consider the time consumption of the communications within an MPI implementation [START_REF] Chan | Collective Communication: Theory, Practice, and Experience: Research Articles[END_REF].

In general the cost of communication between processes (in our case, the total number of processes is the number of subdomains) can be described by

T (n) = α + nβ (2.91)
where n :number of data items.

α :startup time.

β :transmission time per data item.

(2.92)

The start-up cost is due to both hardware and software overhead on the sending and the receiving process. Typically, α is four to five orders of magnitude greater than β, where β is on the order of the cost of an instruction.

All this considerations, gives us some general lower bounds in time for the communications, the most important one in the start-up time, which depends on the number of processes (subdomains), and is independent of the type of communication, either send-receive or reduce operations. This is

T (n) ≥ log 2 (p)α, p : processes (2.93)
For this reason, avoiding communication is a main issue in parallel programming and therefore, in the S-FETI method, big part of the speedup comes from the reduction in the total number of exchange phases (as they are the same as in FETI per iteration, but the total iterations are reduced).

2. The addition of a Cholesky factorization of the small N T × N T dense symmetric positive semidefinite matrix W T p FW p is now needed. In the next section we will give more details about this and we will show an alternative to this procedure.

3. The number of stored directions is increased. Since the number of search directions increases in every iteration, we now store block of directions W p and FW p instead of single vectors.

A priori this difference may not be very significant, but as we will see later, it will turn in a mayor drawback for some cases, specially when using the modified version of S-FETI, where even more directions are created.

4. In general, the most costly part of a FETI algorithm is the computation of solutions of the local Neumann and Dirichlet problems in each subdomain.

In the case of the Dirichlet problem for the preconditioning, there is no difference between the methods, but the F operator must now be applied to the N T columns in W p , meaning that, a priori, N T local Neumann problems should be solved at every iteration. However, this computations can be performed efficiently, see Table 2.1.

As already said in the first point, communication cost is mainly driven by the start-up time, and when multiplying F by a block instead of a vector, the number of interchanges does not change, only the size of it, meaning that we will avoid the significant start-up time, for each iteration reduced with this method.

Another point to consider is that block operations are proportionally less expensive than single vector ones -e.g. N times a matrix-vector product versus a matrix-matrix product with N columns in the second matrix-, this is because the computation time is driven by the memory access. More details of this will be in the following implementation part.

Finally, but more important in our case, is the clever use of the locality of data. Lets note that the matrix Z p+1 is a sparse matrix, this is due to the fact that each column associated to the preconditioner coming from the subdomain Ω (s) is non zero only in the local interface Γ (s) (or Γ (sq) ) while in contrast W p+1 it is not, because of the projection and orthogonalization. Moreover, we can observe that

FW j+1 = FP Z j+1 + j i=0 FW i Φ ij = FZ j+1 -FAG[G T AG] -1 G T Z j+1 + j i=0 FW i Φ ij (2.94)
where we remember that A is a symmetric matrix, that can be the taken as the preconditioner, some scaling matrix or the identity matrix in the simplest case.

The theoretical optimal choice for A is to use the preconditioner, either Dirichlet or Lumped A = D -1 or A = L -1 respectively. Both choices are rather computationally expensive compared to some other alternatives.

From here on, we will choose

A = N s s=1 B(s) diag K (s) bb B(s) T (2.95)
This is the so called Sumper Lumped scaling, which has the property of being computationally inexpensive. This comes from the fact that we only need to multiply locally by a diagonal matrix. At the same time, the convergence ratio of the FETI method is improved, see [START_REF] Farhat | Implicit Parallel Processing in Structural Mechanics[END_REF].

Previous equation shows that in every iteration we only need to solve the localized problems FZ p+1 using the sparsity of Z p+1 and the computation of the projection is done by computing the also sparse matrix (FAG) during the initialization once and for all.

More details of this computations will be given in the details of the implementation to come next.

Remark: With all this differences, the final extra cost per iteration remains very small when N T is not too large. Also, the total of Neumann problems to solve in each subdomain can be done simultaneously as a block, greatly reducing the cost of the method. This extra costs needs to be compared with the final number of iterations given by the method, where a reduction of the order of N T is expected, making the total time of this method an improvement compared to the computation time in FETI-1LM.

To understand how some of the computations are done in order to improve the efficiency of S-FETI, we can now give some details of the practical implementation 1. Lets return to the last point in the cost of the S-FETI method, the simulta-neous forward-backward substitutions. By definition, the columns of Z p related to Ω (s) (the s-column in the regular S-FETI) are nonzero only on the interface Γ (sq) , so given a column in Z p , its product by F requires solving a Neumann problem only in the subdomain itself and its neighbour Ω (q) . On the other hand, given one subdomain Ω (s) , the workload associated with the computation of FZ p is N Neumann solves, where

Subdomain division ⇒ N = neighbours + 1 Interface division ⇒ N = 2 × neighbours (2.96)
In particular, this is much fewer than the rank of Z T p , which is in the order of the number of subdomains (the exact number of subdomain in the first version).

If we compute the multiple local solutions at once, the efficiency is greatly improved on a multi-core machine. In fact, the difference between doing N forward-backward substitution for the separate vectors and doing a single substitution with N -rhs (right hand side) lies in the number of memory accesses between them. In the second case we do almost the same number of accesses than a single rhs (which is the number of non zero entries of the factorized matrix), meaning that even if we multiply the arithmetic complexity by N we will not have the bottleneck that is the memory access.

To understand this we need to see how the memory access on a single or multi-core machine works. On a single core machine, we identify 5 levels of memory from which the data must move before doing a computation at the top. This leveles can be seen as a pyramid, where the top of the pyramid contains the memory that is the fastest, and also the smallest, and the bottom of the pyramid contains the memory that is the slowest but also the largest.

The best routines to take advantage of such a pyramid are the Basic Linear Algebra Subprograms of level 3 kernels (BLAS3). The reason they can take the greatest advantage of this hierarchy is that their memory interaction can be organized in a fashion such that it takes advantage of the various Number of cores Substitutions Time (s) 1 1 0.7 12 12 1

Table 2.1 -Time for forward-backward substitution on multi-core processor.

stages of the pyramid, this leads to have an order of magnitude more computation than memory interaction. In Level-2 BLAS we have O(n 2 ) data moves and O(n 2 ) operations, whereas in BLAS3 we have the same number of data moves O(n 2 ) but we improve the operations to O(n 3 ).

The pyramid concept does not completely map to modern multi-core optimization, however the difference is that on multi-core, there are often many cores at the top of the pyramid, and they may have a shared memories between them.

As an example given in [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF] we can name the case of a sparse matrix of dimension 2 × 10 5 on a 12-core Intel Nehalem processor, Santa Clara, California, US. In Table 2.1 we can see the time for a single forward-backward substitution versus 12 simultaneous substitutions on 12 cores using the PARDISO solver [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with {PARDISO}[END_REF]. With this results we can see that a good local optimization reduces the impact of the extra computations in S-FETI on multi-core machines.

2. Another important issue is the parallel implementation of the projector P . Lets recall its definition

P = I -AG(G T AG) -1 G T (2.97)
as previously said, we will consider A as the super lumped scaling to avoid expensive extra computational cost, the projector is then denoted by P A .

Even if this projection is not local, it only performs a low-rank correction. With this, the computation of P Z p implies a small extra cost.

In practice what it is done, is first consider

P A Z p = Z p -(AG)α p where α p solves (G T AG)α p = G T Z p (2.98)
This will guarantee that G T P I Z p = 0. We recall the definition of G G = B (1) t (1) R (1) . . .

B (N s ) t (N s ) R (N s ) (2.99)
where R (s) = ker(Ω (s) ). Comparing the definition of G and Z p we can see that their sparse pattern are very similar, in fact in both cases for each column associated to the subdomain Ω (s) the pattern is the same, meaning that, given a column

Z (sq) p of Z p , G T Z (sq) p
is computed by applying only dot products by the columns of G corresponding to Ω (s) and its neighbours.

The matrix (G T AG) is factorized during the initialization, so the computation of α p is done via a forward-backward substitution. Furthermore, it can be done in parallel, each subdomain Ω (s) solving for his own columns (the ones associated to Γ (sq) , q = 1, . . . , n (s) ) the following system

(G T AG)α (sq) p = G T Z (sq) p (2.100)
We know that α p is a dense matrix with the same number of rows that the rank(G T ) and with N T columns (the same as Z p ), but in any case, once it is computed we only require a low-rank correction of Z p in each subdomain, because only a few columns of G are nonzero in Ω (s) , with this we compute the total correction

P A Z p = Z p -AGα p (2.101)
To understand better this procedure, and the general implementation for the treatment of each sparse block, we will define what we call the coarse modes, this modes correspond to the non-null vectors in each one of the columns of the matrices Z, G, FZ, FAG. They are stored in each local interface of each subdomain, and have a correspondence with the named global matrices.

In the left part of Figure 2.5 we show from a subdomain point of view the different coarse modes stored. They are the local (but shared with neighbours) modes that describes Z or G, so for each interface we save one local mode and one from its neighbour. This double storage allows to perform computations of, for example G T Z, in each subdomain by only performing dot products between the coarse modes saved in each interface of each subdomain. This way of storage is generalized to allow local computations between this 4 matrices only by performing operations between the different modes stored this way.

In the case of the multiplication by the operator F, we perform in each subdomain a forward-backward substitution of each one of the modes previously described or in the basic case they are considered as an single mode (when doing the local subdivision it is 2 × neighours), and then also storing the resulting (more numerous) modes in each interface.

Also in the left drawing we see the subdomains sharing the interfaces involved in multiplications by the F operator.

In the right part of Figure 2.5 we show the coarse modes that describes FZ or FAG stored in one of the interfaces of certain subdomains. For each local interface we save the modes corresponding to the local modes after the forward-backward substitution plus the modes of the neighbour after the substitution.

Again, this implementation in form of coarse modes, allows the computations of matrix products such as G T (FAG), only by performing dot products between this modes. This implementation also explains how the single multi-rhs substitution is done, just by arranging (copy) the coarse modes in a block matrix. The same way we can improve multiple dot products computations by arranging the modes in order to perform matrix-vector products or matrix-matrix products that improve the performance of this type of computations (See previous point).

3. The previous process it is also important in the calculations to obtain FW p . Lets recall that the orthogonalization process that leads to FW p is applied to the block FP Z p , so we need first to compute this block. To do so, we use

SUBDOMAIN INTERFACE

Ω (s) Ω (q) Ω (s) Figure 2.

-Subdomain point of view (Left):

In red the coarse modes describing Z or G owned by subdomain Ω (s) , in black the interfaces where F times the red modes is also non null.

Interface point of view (Right):

In doted lines the coarse modes describing FZ or FAG owned by local interface Γ (sq) between subdomains Ω (s) and Ω (q) , in black the modes where, due to Γ (sq) , Z T (FZ) is non null. the matrix α p previously computed in the previous point, the we need to do

FP Z p = FZ p -(FAG)α p (2.102)
Again what we are doing here is a low-rank correction, but this time of FZ p . The difference is that,to do this corrections, we now use the sparse pattern of (FAG) (each column is nonzero in Ω (s) and its neighbours) instead of the one in G. The sparse matrix (FAG) is calculated at the initialization once and for all.

4. The last item is the one regarding the orthogonalization and full reorthogonalization procedures for the search directions. As said in the previous item, once the block P Z p is computed they must be F-orthogonalized to obtain the search directions W p and its corresponding block FW p . This can be achieved by using the classical modified Gram-Schmidt procedure, which will require the use of many communications between subdomains, typically Message Passing Interface (MPI) reductions of dimension 1.

To avoid this many communications we can compute (P Z p ) T (FP Z p ) by using the local contribution of each subdomain obtained with the BLAS3 and then doing only one MPI reduction. The Cholesky factorization is then applied to compute the block of directions W p as in line [START_REF] Barka | Scalability of FETI-2LM methods on HPC clusters for antenna RCS applications[END_REF] of Algorithm 9.

As for the full reorthogonalization process, we use a block modified Gram-Schmidt procedure with this same BLAS3 strategy in order to have a number of MPI reductions equal to the ones used in classical FETI, but changing from one vector to a matrix and also changing the size of the reductions from 1 to at most N T × N T .

The evaluation of all this considerations for a good implementation will be shown in the numerical results at the end of this chapter.

Sorting search directions in S-FETI

The block version of the FETI method is proved to be a great improvement in comparison with his basic version, all based in the idea of enlarging the search space of the solution. As a consequence of this extension, we have changes in the algorithm that will impact the execution time. We already name in the previous chapter the main differences between them, we studied the cost difference between FETI and S-FETI where we point out the need to inverse the matrix W T FW .

This matrix is in general positive semi definite because of the linear dependency that may occur in the directions built so we can only aim to compute a pseudo inverse, but for this same reason it can also be used to sort the different directions between the useful ones and the rest. Even if, a priori, the computation of this pseudo inverse is cheap in terms of time consumption his impact in the final results of the algorithm needs a detailed study.

Linear dependence in block FETI directions

We start this section by recalling the search block built from the preconditioner, that is computed in every iteration as

W p = P Z p + p-1 i=0 W i Φ ip (2.103)
This block is orthogonalized with respect to previous blocks W i already stored. Also the columns of each of the stored block, need to form an F-orthonormal family, which is the norm used in CG algorithm, i.e. for every block they need to fulfill

W T p FW p = I (2.104)
From this we know that we also have to find a decomposition of W T p FW p that allows us to modify the basic blocks to build this orthonormal family, thus the need to study the properties of this matrix.

Since F is symmetric it's clear that this matrix is also symmetric, but it is not so evident that will be positive definite. The lack of this property comes from the fact that a linear dependence between the columns of W p can appear at certain point of the main iteration. This dependence comes mainly from two issues, the first is the convergence in some of the local interfaces and second is the working precision.

To give a better understanding of the first issue, we have to consider that the global interface is built from the smaller local interfaces that exist between two neighbour subdomain, and the characteristics of each of this local interfaces varies from one to another, all being of variable size (as in number of d.o.f.) and also with different mechanical properties along the same. Both characteristic have an impact in the convergence of the method.

For the simpler or smaller local interfaces we expect a faster convergence, for example in a 2D elasticity problem, if the interface is made of 1 node (e.g. in a corner) the convergence is achieved in 2 iterations, so the columns created from this interface (in the subdivided version of S-FETI) will be no longer useful and the linear dependence appears.

For the problem of representation in limited precision arithmetic, depending on the configuration, but also closely related to the phenomenon of achieved convergence, we can have directions whose differences are so small that it's difficult to distinguish it from numeric noise, meaning that in practice we have the same direction, so the use of them is redundant.

With both previous issues in consideration it is clear that the matrix W T p FW p is only positive semi-definite so no full inverse exists. We have two choices from here, either compute a pseudo inverse and use the complete block with redundancies, or we can sort the directions in order to build a smaller but full block with the exact same information but less computational effort in the operations to follow.

The second choice is in practice the best one. To make this sorting, in [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF] a Cholesky factorization is proposed, but we can also add a second way that consist in apply a diagonalization process to the matrix W T p FW p . Both procedures will be detailed in the next sections.

Cholesky factorization with complete pivoting

For the process of F-orthonormalization done to the block of search directions, as already said, we start with the decomposition of the W T p FW p matrix. Since this matrix is symmetric we can use the Cholesky decomposition, instead of a regular QR decomposition, meaning that at each step of the algorithm a lower triangular matrix L p is built. This matrix is such that

W T p FW p = L p L T p (2.105)
This decomposition only works if the previous matrix is of full ranking, since this is not our case, we need to use the alternative Cholesky decomposition with complete pivoting (symmetric pivoting) described in [START_REF] Hammarling | LAPACK-Style Codes for Pivoted Cholesky and QR Updating[END_REF]. The algorithm for this factorization is proved to work as a rank revealing procedure. This is achieved by computing an square permutation matrix N , in this case the decomposition can be written as

N T W T p FW p N = L p L T p (2.106)
and if we look in detail the matrix L p we know that it can be described as

L p =       Lp 0 × 0       (2.107)
where Lp ∈ R r×r (2.108) this smaller matrix is lower triangular with positive elements in the diagonal and it has a full rank r ≤ N T , which is at the same time the rank of the W T p FW p matrix.

The existence of the permutation matrix (proved in [START_REF] Nicholas | Accuracy and Stability of Numerical Algorithms[END_REF], Thm. 10.9) allows to sort the directions that are useful from the rest. This is simply done by taking the first r columns of W p N .

In summary, at each iteration of the method, the directions to be used and stored are defined as

W p ← (W p N ) L-1 p 0 T (2.109)
this block of sorted directions fulfills the 2.104 condition, in fact

L-1 p 0 (W p N ) T F(W p N ) L-1 p 0 T = L-1 p 0 L p L T p L-1 p 0 T = L-1 p Lp 0 L-1 p Lp 0 T =       I 0       (2.110) with I ∈ R r×r .
The algorithm for the computation of the Cholesky factorization, introduces the use of a new parameter ε > 0 close to zero, to determine at which point we are in the presence of linear dependence. If this parameter is too small we can be working within the numerical noise and losing precision in the total convergence, but on the contrary using one that is too big can have the impact of losing directions that may be important for the convergence.

Finding the optimal ε is not easy task and is usually chosen based in numerical experiences, but as a basic consideration, we'll try to be closer to the "losing" directions approach because, even if we can have a slower convergence, eventually we will recover this directions in later iterations. On the other hand, working within the range of the natural noise of the method, can only introduce sources of error harder to control, and may even be ending in blowups.

Usually a regular double precision implementation should be enough, but introducing new sources of numerical instabilities needs to be done carefully, so later in numerical tests we'll compare the use of quadruple precision in the computation of the Cholesky decomposition for this matrix, as an alternative for accurately capture the directions needed.

Diagonalization of search directions block

The second way to treat the sorting of independent search directions can be done using a diagonalization process. At first we now that this procedure will lead to the decomposition

W T p FW p = E p D p E -1 p (2.111)
Where E p is the matrix of the eigenvectors in the columns, and D p is a diagonal matrix with the eigenvalues in the diagonal. In our case, the matrix is real and symmetric so we have that E -1 p = E T p . Hence, if we are looking for a block that fulfills the condition 2.104 we can use the previous decomposition to build it. First, we know that

W T p FW p = E p D p E T p (2.112) this implies that E T p W T p FW p E p = D p (2.113)
and since the matrix D p is diagonal, we can easily find its inverse and even more we can divide it into the multiplication of two diagonal invertible matrices, meaning that we build

E T p W T p FW p E p = D 1/2 p D 1/2 p (W p E p D -1/2 p ) T F(W p E p D -1/2 p ) = I (2.114)
And so our F-orthogonalized search directions can be constructed by doing

W p ← W p E p D -1/2 p (2.115)
The previous diagonalization process works if the matrix W T p FW p is invertible, but we only have that is positive semi-definite, so the eigenvectors are greater or equal than zero.

The common algorithms for the computation of eigenvalues can be used to obtain them in descending (or ascending) order. This fact will allow us to build a smaller F-orthonormal full block of search directions by taking into account the biggest of this eigenvalues and their associated eigenvectors, so again, the construction is the same as previous, but we will only consider the r biggest eigenvalues to build then the block

W p ← W p Ẽp D-1/2 p (2.116)
where

rank(W T p FW p ) = rank( ẼT p ) = rank( Dp ) = r (2.117)
meaning that no information will be missing from this sorting.

The treatment of the linear dependencies, is still limited by some ε > 0 parameter, but now this is done in a clearer way, because the algorithm for the computation of the eigenvectors and eigenvalues will give them in decreasing order, so the biggest values, and thus the most important ones, are the first columns of E p so we can eliminate the smallest eigenvalues and vectors below the threshold ε. In this case the dependence between this parameter and the final result is less important than in the cholesky factorization because we do not eliminate certain directions, but instead we build a smaller full block containing all the relevant information.

The impact in time of the diagonalization process is minimal, because of the size of the matrix we are solving depends on the number of subdomains. Also the use of LAPACK algorithm with BLAS optimization, allows a fast computation in every subdomain. Remark: Later in the numerical results we will also show the impact in the use of the Lumped preconditioner, as an economical alternative to add into this new decomposition.

Memory usage in S-FETI

In previous sections we stipulate that one of the difference between the S-FETI method and the classical FETI lies in the total memory use of both of them. For the classical FETI, only two vectors are stored at each iteration, they correspond to the search direction and its multiplication by F. On the other hand, in S-FETI the number of search directions built in every iteration is increased, making mandatory the storage of blocks W p and FW p . This difference is in fact one main drawback in the implementation of S-FETI.

The total number of directions stored depends and grows linearly as the number of iterations, this value is in every step usually the number of columns in Z p (in practice, due to linearity of the directions, we expect a smaller number that correspond to Rank(Z T p )). In 2D problems this is not much of an issue, but in 3D even in small cases this needs to be considered, for example if our domain is a small cube divided into 5 × 5 × 5 smaller cubes, we see in 2.2 that the number of stored directions at each step is of 1600 directions approximately, a fact that comes only from the geometry of the problem. Now for the total memory needed this will depend on the size of the triangulation of the mesh, in this case of about 10 5 elements in each subdomain.

To be able to use this method in bigger configurations we will change the usual, direct storage and use a different strategy. We know that the time of a parallel application is limited by the exchange of information between processes, so if we add some, not too expensive, extra computation we hope to have an implementation as fast as the original one, but without the memory constraint.

New sparse storage

The idea of this new storage is to, reconstruct the search directions at every iteration, that is, to compute the projection and the full reorthogonalization of the columns in Z p , by storing matrices that are either sparse or small ones. This is done by saving at each iteration the sparse matrices, Z p and FZ p , instead of the full W p = P Z p and FW p = P FZ p . With this sparse matrices we also compute and store some smaller coefficient matrices that will allow this reconstruction. The size of this new matrices will depend on the number of used search directions at each iteration, and will be independent of the size of the mesh, so in general, it will be a small number compared to the size of the problem.

The memory needed for the storage will now be limited in a more important way on the problem configurations, that determines the maximum total search directions in every iteration, rather than the size of the discretization. Compared with the total unknowns of the problem, this number (N T ) is usually a small one, because the number of neighbours of each subdomain is also limited, so the problems that we can solve can increase in its total size and have a smaller impact in the memory used by it, allowing us to use the S-FETI method in bigger problems.

The fact that we no longer store the complete search directions will also have an impact in the precision achieved with this implementation, but we hope to keep the good accuracy of the S-FETI method or at least be close enough to make this a useful algorithm, we will later discuss this fact in the numerical results.

Reconstruction of search directions

In this section we will show the details used by the algorithm of the sparse storage for the S-FETI method.

Idea of the reconstruction

Use the definition of the projector and full reorthogonalization process, to try to find a recurrence for the coefficients that define each search direction. Lets consider the sparse matrices Z p from 2.89 and lets recall the definition of the projector

P = I -AG(G T AG)G T (2.118)
with A a symmetric matrix, that can be the identity, the preconditioner or some scaling matrix. We also know that to apply this projector we use low-rank corrections as in 2.101. Starting from Z p we have

Z 1 =⇒ P Z 1 = Z 1 + (AG)D 1 (2.119)
its clear that

G T P Z 1 = 0 (2.120)
and that the computation of D 1 is done by solving

(G T AG)D 1 = -G T Z 1 (2.121)
meaning that the value of D 1 is

D 1 = -(G T AG) -1 G T Z 1 (2.122)
Since no previous directions were computed, we have our first search direction block, and the projection coefficients that define it

W 1 := Z 1 ∆ 1 := D 1 P W 1 := P Z 1 = Z 1 + (AG)∆ 1 (2.123)
In the next iteration and from the definition of the S-FETI method, we build the sparse block Z 2 , then we apply the same process for this block

Z 2 =⇒ P Z 2 = Z 2 + (AG)D 2 (2.124)
where

G T P Z 2 = 0 (2.125)
and D 2 is computed by solving

(G T AG)D 2 = -G T Z 2 (2.126)
so again, its value is

D 2 = -(G T AG) -1 G T Z 2 (2.127)
Now, we need to consider the full reorthogonalization against the directions built in the previous step. Hence we have

W 2 = Z 2 + W 1 Φ 12 =⇒ P W 2 = P Z 2 + P W 1 Φ 12 (2.128)
and to compute Φ 12 we use the orthogonality properties of the new directions and the previous ones in the CG algorithm

(FP W 1 ) T P W 2 = 0 (2.129)
replacing P W 2 we have

(FP W 1 ) T (P W 1 )Φ 12 = -(FP W 1 ) T P Z 2 (2.130)
this problem is solved in order to compute Φ 12 whose value is

Φ 12 = -(P W 1 ) T (FP W 1 ) -1 (FP W 1 ) T P Z 2 (2.131)
then, in 2.128 we use the projection definition again, to obtain

P W 2 = Z 2 + (AG)D 2 + (Z 1 + (AG)D 1 )Φ 12 = Z 2 + Z 1 Φ 12 + (AG)(D 2 + D 1 Φ 12 ) (2.132) so P W 2 = Z 2 + Z 1 B 12 + (AG)∆ 2 (2.133)
where

B 12 := Φ 12 and ∆ 2 := D 2 + D 1 Φ 12 (2.134)
but also, since

W 2 = Z 2 + Z 1 B 12 (2.135)
the we have that

P W 2 = W 2 + (AG)∆ 2 (2.
P W p = Z p + p-1 j=1 Z j B jp + (AG)∆ p (2.139)
With the previous recurrences in mind, we will formalize and explicit the computation of the coefficient matrices B jp and ∆ p that needs to be stored. This matrices along with the sparse ones Z p define the construction of the search directions.

Remark: From previous construction we can extend the sparse multiplication against the operator F, namely from Equation 2.133 we have

FP W 2 = FZ 2 + FZ 1 B 12 + (FAG)∆ 2 (2.140)
and then in the general form, using Equation 2.139 we have

FP W p = FZ p + p-1 j=1 FZ j B jp + (FAG)∆ p (2.141)

Formal construction

At the p iteration, we apply to Z p the low-rank correction definition of

P P Z p = Z p + (AG)D p (2.142)
with

D p = -(G T AG) -1 G T Z p (2.143)
also from the definition of P , we apply it this time to the reorthogonalized block

W p P W p = W p + (AG)∆ p (2.144)
later we will compute ∆ p , not from the classical definition, but recursively.

On the other hand and to avoid notation problems, lets consider some iteration j. We replace in the previous equation 2.144 the computation of W j (W p in there) by the Gram-Schmidt orthogonalization process applied to

Z j P W j = Z j + j-1 i=1 Z i B ij + (AG)∆ j (2.145)
with B ij to be computed later, recursively. If we replace this in 2.138 that actually correspond to the formal construction of the orthogonalized and projected search directions, then we have

P W p = P Z p + p-1 j=1         Z j + j-1 i=1 Z i B ij + (AG)∆ j         Φ jp = Z p + (AG)D p + p-1 j=1 Z j Φ jp + p-1 j=1 j-1 i=1 Z i B ij Φ jp + p-1 j=1 (AG)∆ j Φ jp =          Z p + p-1 j=1 Z j Φ jp + p-1 j=1 j-1 i=1 Z i B ij Φ jp          + (AG)          p-1 j=1 ∆ j Φ jp + D p          (2.146)
from this equation and 2.144 it's clear that we can define the computation of W p and ∆ p as

W p = Z p + p-1 j=1 Z j Φ jp + p-1 j=1 j-1 i=1 Z i B ij Φ jp (2.147) ∆ p = p-1 j=1 ∆ j Φ jp + D p (2.148) Developing the computation of W p W p = Z p + p-1 j=1 Z j Φ jp + p-1 j=1 j-1 i=1 Z i B ij Φ jp = Z p + p-1 j=1 Z j Φ jp + p-2 i=1 Z i         p-1 j=i+1 B ij Φ jp         = Z p + p-1 i=1 Z i Φ ip + p-2 i=1 Z i         p-1 j=i+1 B ij Φ jp         (2.149)
and thus our orthogonalized search direction is

W p = Z p + p-1 i=1 Z i B ip (2.150)
with

B ip := Φ ip + p-1 j=i+1 B ij Φ jp (2.151)
Finally, using the computation of W p from 2.150 and ∆ p defined recursively in 2.148, the projected search directions at the p iteration can now be computed as

P W p = W p + (AG)∆ p (2.152)
So we can now reconstruct the directions only from the sparse matrices Z j , j = 1, . . . p and the small matrices B ij , ∆ j and Φ ij , i = 1, . . . , j, j = 1, . . . p that are stored instead of the full larger matrices W p or FW p .

At each iteration with start with the Z p matrix that is built from previous gradient and with it we can easily compute the coefficients D p that do not need to be stored (we can, for example use the same memory space that uses ∆ p and later overwrite it).

To define the computation of the coefficients in Φ ip we need to add some extra calculus, that cannot be based in the storage of any mesh dependant full vector. To achieve this, lets consider at each iteration the generalization of 2.129, thats based in the general orthogonal properties of the CG algorithm (FP W j ) T P W p = 0 for j < p

(FP W j ) T (P Z p + p-1 j=1 P W j Φ jp ) = 0 (2.153) which implies (FP W j ) T (P W j )Φ jp = -(FP W j ) T P Z p for j < p (2.154)
The left side matrix correspond to the matrix symmetric positive semi-definite that needs to be pseudo inverted at each iteration. To simplify the computations we apply the cholesky decomposition to it, giving

L j L T j Φ jp = -FP W j T P Z p (2.

155)

The L j matrix need to be stored at every iteration, but then again this are small triangular matrices of size rank(Z T j ).

The right hand side, is the one needing a special treatment as it depends on full terms that we do not want to store. For the P Z p term, the definition of P is enough to avoid the full storage, as for the term FP W j , instead of using it directly, we also need to use the reconstructions that comes from the sparse matrices.

From the construction of P W j in 2.145 and the fact that F is a linear operator, we have

FP W j = FZ j + j-1 i=1 FZ i B ij + (FAG)∆ j (2.156)
From this, as expected, we also need to store the sparse term FZ j . Finally the right term in 2.155 is now computed

FP W j T (P Z p ) =         FZ j + j-1 i=1 FZ i B ij + (FAG)∆ j         T         Z p + (AG)D p         = Z T j (FZ p ) + j-1 i=1 B T ij Z T i (FZ p ) + ∆ T j (FAG) T Z p + + Z T j (FAG)D p + j-1 i=1 B T ij Z T i (FAG)D p + ∆ T j (AG) T (FAG)D p (2.157)
We have to compute the six previous terms at every iteration for j = 1, p -1. With this we can know obtain from 2.155 the values of Φ jp and the new implementation is completed.

Next we present a summary of the computations needed and the matrices stored in this new algorithm for the S-FETI method. Z p = [. . . , F(s) -1 Γ (sq) P T g p , . . . ], s = 1, N s q = 1, . . . , n (s) 6:

Q p = FZ p 7: D p = -(G T AG) -1 G T Z p 8:
for j = 0 to p -1 do 9:

S jp = Z T j Q p 10:
end for 11:

T p = Z T p (FAG)

12:

for j = 0 to p -1 do 13:

Φ = S T jp + j-1 i=0 S T ip B ij + T p ∆ j 14: Φ = Φ T + T j D j + j-1 i=0 B T ij T i D p + U T j D p 15: Φ jp = -(L j L T j ) -1 Φ 16:
end for 17:

∆ p = D p + p-1 j=0 Φ jp 18:
U p = (AG) T (FAG)∆ p 19:

for j = 0 to p -1 do 20:

B jp = Φ jp + p-1 i=j+1 B ji Φ ip 21:
end for 22:

W = Z p + p-1 j=0 Z j B jp + (AG)∆ p Compute projected
and reorthogonalized blocks 23:

FW = Q p + p-1 j=0 Q j B jp + (FAG)∆ p 24: N p L p L T p N T p = W T (FW ) 25: ρ = -( Lp LT p ) -1 (W N p ) T g p Lp ∈ R r×r is full ranked 26: λ p+1 = λ p + (W N p )ρ 27:
g p+1 = g p + (FW N p )ρ

28: end loop

Compared with the S-FETI method, we just replaced the projection and full reorthogonalization processes to build W and FW , the rest is analogous, including the use of different preconditioners as well as the rank revealing strategy.

Implementation details and exploitable parallelism

In previous section, in order to save memory usage, we introduced several extra computations to the basic S-FETI algorithm. If we want this method to be as fast as the full version, all this new computations needs to be implemented in the most efficient way possible. The algorithm shown in 10 describes a simple presentation of the construction of the search directions, that can be implemented in a straitforward way, but here we will give some considerations that can give an important time boost to the implementation of this method.

Before explaining this, we will explain the sparsity of the new extra terms 2.158 that we need to compute and store, lets recall them

(G) T (FG), Z T j FZ p , Z T p (FG) (2.160)
with j = 1, . . . , p -1 and for simplicity we consider A = I. If we look at this three terms we note that each one shares very similar structures, this comes from the definitions of Z and G G = B (1) t (1) R (1) , . . . , B (N s ) t (N s ) R (N s ) Z = B (1) ∆f (1) , . . . , B (N s ) ∆f (N s )

(2.161) again for simplicity, no scaling is considered in Z, we consider the complete local interface in Z (without any subdivision) and we think in R (s) = Ker(Ω (s) ) as a single vector in Ω (s) . In both cases the extension to local interface division and kernels with more than one vector is straightforward, because each column associated to a subdomain s can be consider to have the exact same structure.

With this definitions we see that G and Z have the same sparse structure given by the definition of boolean matrices B (s) , to see it more clearly, lets take as example the Figure 2.6, where the domain is divided in 6 subdomains in a b , and divide a vector w along the total interface into local interface blocks

w 1 w 2 w 3 w 4 w 5 Ω 1 Ω 2 Ω 3 Ω 4 Ω 5 Ω 6
w =                       w 1 w 2 w 3 w 4 w 5                       (2.162)
we see that the size and structure of w is the same as

B (s) v (s)
b , for any s = 1, . . . 6. With this, each column of G and Z, defined in the whole interface, have the following structure

B (1) v (1) b =                       × 0 0 0 0                       , B (2) v (2) b =                       × × 0 0 0                       , B (3) v (3) b =                       0 × × 0 0                       , B (4) v (4) b =                       0 0 × × 0                       , B (5) v (5) b =                       0 0 0 × ×                       , B (6) v (6) b =                       0 0 0 0 ×                       , (2.163)
Here we see that in each column only a small part of the interface is non zero, and it correspond to the d.o.f. in the column s shared by the subdomain s (Ω (s) ∩ Γ 0).

The next structure to consider is the one of the terms FG and FZ, that are the multiplications of previous blocks by the FETI operator. From the definition of F we can see the structure of each previous blocks multiplied by F depends on the multiplication of the boolean matrices B (s) T B (q) , ∀s, q = 1, . . . , N s (2.164)

When subdomains s and q are neighbours, this multiplication is non zero, so when multiplying Z or G by F only local interfaces corresponding to the neighbours are non zero. In our example the structures are

FB (1) v (1) b =                       × × 0 0 0                       , FB (2) v (2) b =                       × × × 0 0                       , FB (3) v (3) b =                       × × × × 0                       , FB (4) v (4) b
=                       0 × × × ×                       , FB (5) v (5) b
=                       0 0 × × ×                      
, FB (6) v

(6) b =                       0 0 0 × ×                       , ( 2.165)
With previous blocks, the final structure for the terms (G) T (FG), Z T j FZ p , Z T p (FG) is given, in this case, by

                            × 0 0 0 0 × × 0 0 0 0 × × 0 0 0 0 × × 0 0 0 0 × × 0 0 0 0 ×                                                   × × × 0 0 0 × × × × 0 0 0 × × × × 0 0 0 × × × × 0 0 0 × × ×                       =                             × × × 0 0 0 × × × × 0 0 × × × × × 0 0 × × × × × 0 0 × × × × 0 0 0 × × ×                             (2.166)
here we see that the terms that are non zero, correspond in each subdomain to the subdomain itself plus its neighbours plus the neighbours of neighbours. This is a general fact that comes from the definitions of Z, G and F

(G T FG) ik = v (i) T b B (i) T       s B (s) S (s) -1 bb B (s) T       B (k) v (k) b (2.167)
in here the term (G T FG) ik is non zero if B (i) T B (s) and B (s) T B (k) are non zero, which occurs when the subdomains i-s and s-j are neighbours.

Even though the last term is less sparse that both terms Z and FZ (G and FG) its total size is much smaller than previous blocks since they are square matrices of size N T .

Optimization of the code

With this structures in mind we can now give the details of several possible optimization to any basic implementation of Algorithm 10. Lets remember that as part of the new orthogonalization process we have to compute the six terms in 2.157, so the improvements come, in one hand, from the fact that some matrix vector computations can be done as block, using the advantages of BLAS3, and on the other hand some of the coefficient matrices computations can be done in parallel.

This changes to optimize the code can be summarized as:

• Concatenated allocation in memory of the terms Z j with j = 1, ..., p -1 so that the computation of Z T j (FZ p ) can now be done as a single block for all j.

To understand this, lets see how the computation of each term Z T j (FZ p ) is done locally. Lets consider a column associated to the subdomain s of Z p , defined in the complete interface, but taking only local values

z (s) p = B (s) f (s) Γ (sq)
(2.168) where f (s) Γ (sq) is defined in 2.87. Then, each value of Z T j (FZ p ) is computed as

(Z T j FZ p ) ik = z (i) T j Fz (k) p = z (i) T j       s F (s)       z (k) p = s z (i) T j F (s) z (k) p (2.169)
where, as usual,

F (s) = B (s) S (s) -1
bb B (s) T and as seen previously, the values that are non zero are the ones where i-s and s-k are neighbours. This implies that in each subdomain s, we need to do an exchange with every neighbour and store in s the vector z (k) p , then we compute, exchange and store F (s) z (k) p at every iteration.

In general, each subdomain s will compute and store the sparse vectors

z (s) p , F (s) z (s) p , z (k) p , F (s) z (k) p , where k ∈ neighbour(s).
Finally each subdomain s will compute the local dot products

z (i) j F (s) z (k) p , ∀i, k ∈ neighbour(s) ∪ {s} (2.170)
then an assembly into (Z T j FZ p ) ik plus a global reduction are done to compute the total matrix Z T j (FZ p ). This are small sparse matrices stored in every subdomain. This computation, is done for each j = 1, . . . , p-1 at the p iteration, meaning, that we repeat this process p -1 times at each iteration, with the corresponding load of doing the p -1 global reductions and (p -1)-times the dot products.

In order to improve the time of this calculations what we can do is to take advantage of block implementation, as already done in this chapter. To do so, in every subdomain Ω (s) we do a pseudoallocation that will concatenate the newly computed and exchanged local vectors z

(i) p into a block of vectors z (i) 1 , . . . , z (i) p-1 , ∀i ∈ neighbour(s) ∪ {s} (2.171)
with this we can compute the terms in Equation 2.170 as a block

z (i) 1 , . . . , z (i) p-1 T F (s) z (k) p , ∀i, k ∈ neighbour(s) ∪ {s} (2.172)
we can use the optimized libraries BLAS3 to obtain a better performance than usual dot products. Since every Z j shares the same structure, a multiple assemblage can be done into the different Z T j FZ p , and, for each j at the same time, only a single reduction is now done (slightly bigger in size). We then reduce the number of global exchanges from p -1 to 1 in each iteration, improving them, and thus the time cost.

• If we look at the following term in the loop to compute

Φ jp j-1 i=1 B T ij T i D p (2.173)
we can see that the matrix D p is independent of the addition, so it can be treated as

        j-1 i=1 B T ij T i         D p (2.174)
so the term between parentheses is independent of the iteration p and can be computed in previous iterations.

Therefore, instead of doing computations of the type (FAG)D p , we can compute and save

V p ←         p-1 i=1 B T ip T i         (2.175)
This computation is added after the matrices B ip are calculated. We store this sum as a single matrix to use it directly in the next iterations. This term is another small matrix, so the memory cost is still controlled.

• Thanks to the first point in here, the matrices needed for the computations of the terms Φ jp for all j < p, are global and shared in each subdomain, including the matrices coming from the Cholesky factorization. Hence, a prori, each process compute the same terms.

This involves the computation of several small matrix products, that increase in numbers at each iteration. To avoid this increasing computations, the calculation of Φ jp for each j = 1, p -1, can be done in parallel.

Each process j can now be in charge of each matrix Φ jp , up to a limit of N s matrices (the total number of processes). This means that after the iteration p > N s the first process will need to do the computation of two matrices Φ jp and then for each extra iteration, another process will be in charge of the new computations. Since the number of iterations of the S-FETI method is expected to be low, there will be some rare cases where this will happen.

Even if, a priori, this parallelization will speed up the total time, we cannot forget that a global exchange is added after, so that every process can have access to the new computed Φ jp matrices. The impact of this optimization will be tested in the numerical examples.

Numerical results

In this section we want to focus in the S-FETI method, by testing the different implementations and changes done to the method that are presented in this chapter. The objective is to see the comparative advantages or disadvantages that can show each one of them.

To begin with, we will solve the Poisson problem in 3D with Dirichlet condition in a part of the boundary. Let Ω ⊆ R 3 be a bounded domain, let also ∂Ω D ⊆ ∂Ω be a part of the boundary where Dirichlet conditions will be imposed. Given f ∈ L 2 (Ω), the problem is

Find u ∈ H 1 (Ω) such that              -ν∆u = f in Ω u = 0 on ∂Ω D ∂u ∂n = 0 on ∂Ω \ Ω D (2.176)
with ν ∈ R + , a parameter that characterizes the medium or material and therefore will produce the differences in the local stiffness matrices of the subdomains when a partition of Ω is created.

Let V := {v ∈ H 1 (Ω) : v| ∂Ω D = 0}.
Then, the weak or variational formulation of previous problem is

Find u ∈ V such that ν Ω ∇u • ∇vdx = Ω f v, ∀v ∈ V (2.177)
This problem is discretized using bi or tri-linear Q 1 finite elements functions, leading to the known global system of equations

Kx = f (2.178)

FETI and S-FETI

Using the system arising from the finite element method for the Poisson problem we will start by revisiting some of the numerical results shown in [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF], but adding some small differences that we will consider to be our basic starting S-FETI algorithm when doing the comparative performance of the different implementations.

In terms of the method itself, the local interface subdivision will be the default, mainly because we want to reduce to a minimum the number of iterations by augmenting the number of search directions created each step. Also a small variation will be used in the detection of the kernels and the solutions of local problems, usually the PARDISO solver ( [START_REF] Cosmin | An augmented incomplete factorization approach for computing the Schur complement in stochastic optimization[END_REF], [START_REF] Cosmin | Real-time stochastic optimization of complex energy systems on high-performance computers[END_REF]) for the sparse local solutions is used in this type of implementation, but in our case we will change this to the use of a different solver, the one called DISSECTION Solver. This different solver allows a more robust computation of the kernels and recently it has shown an improved general performance for double and quadruple computations, compared to the usual Intel PARDISO implementations, see [START_REF] Suzuki | A dissection solver with kernel detection for symmetric finite element matrices on shared memory computers[END_REF].

In this first case, for the 3D Poisson problem, we model the cube Ω = [0, 1] 3 , with null Dirichlet boundary in ∂Ω D = {(x, y, z) ∈ ∂Ω : x = 0 ∨ x = 1}, also f = 1.

We divide the cube into 125 smaller cubes, and define ν = 10 0 for the whole domain, see Figure 2.7, then each subdomain is divided into 75 thousand elements approximately, We consider as a stopping criterion the same as previous chapter, that is the global relative error and the relative solution jump norm across the interfaces. In some cases precision measures will be performed, if that is not the case the criteria will be to stop when both errors are less than 10 -6 .

We consider the "Super lumped" scaling for the projection, i.e P = P A = I -AG(G T AG) -1 G T with A as in Equation 2.95, and the consistent weighted Dirichlet preconditioner as the global preconditioner.

The results are noted in Table 2.3. We can see here that the number of iterations is greatly reduced, but at the cost of adding also the extra search directions. This change tries to reduces the impact of communications between processes which is usually the bottleneck in parallel algorithms. In the next case, if we use the local interface subdivision that generates more search directions per iteration, the difference between solvers can be exposed, as the number of local solutions is augmented and therefore the total time cost will depend more on the solver than in the regular S-FETI or simply FETI.

S-FETI

The following test, as well as all the time performance test of this chapter, were done with a fortran-mpi implementation, on a machine SGI UV 2000 with 32 CPUs Intel Xeon 64 bits EvyBridge E4650 of 10 cores each one, with a frequency of up to 2.4 GHz. In the graphic of Figure 2.8 we can see that independent of the size problem, the total time is reduced when using the Dissection solver. This speedup however is at an extra memory cost. From now on, both local interface subdivision and Dissection will be the default for the following S-FETI test.

Decomposition of W t FW

Following the propositions of section 2.2 in this part we will compare the different strategies to sort the search directions built in S-FETI. We can name both Cholesky and Eigenvalues Decompositions of the matrix W T FW .

For this problem, we use the same cube as before, with 125 smaller cubes, but we make an heterogeneous case, where we change the parameter ν to be 10 0 and 10 3 in a checkerboard configuration, as seen in Figure 2.9 and we use the local interface subdivision for the S-FETI method.

In this configuration, what we want is to produce a "high enough" number of search directions that can allow us to see the differences between sorting methods and the sensibility to the zero value parameters. Is this process which introduces new sources of numerical instability, so as a first approach, although very time costly, the use of quadruple precision in the decomposition process may improve the quality of the same and possibly reduce the numerical noise. We test for different values of ε and we check the number of directions kept as In Table 2.4 we see that as the parameter ε is closer to zero, clearly the number of directions saved augments, and thanks to the full reorthogonalization we also gain in iteration numbers, but the extra directions used are not entirely the ones that we hope to capture with this method, this is show in the minimal error column where bigger values are obtained if we use more of them. Also a small gain in precision is obtained when we use the quadruple precision, but since the behaviour of both implementations is the same and the cost of the quadruple computations is much more elevated, we will perform only double implementations with the already named value for ε.

A significant more precise method (smaller minimal error) is obtained when we do not consider values of ε too close to zero, meaning values that may be within the limits of just noise, so a value of, for example ε = 10 -4 is a much better alternative than ε = 10 -6 or any value smaller than that.

Remark: The solver Dissection also comes with a quadruple implementation for solving the local problems and doing the computation of the kernel. However there were no significant advantages in this context that could suggest the use of this still expensive new implementation. Simple Jacobi one step improvement used to augment the precision of the local solver are enough to obtain the same results as a complete quadruple computation performed by Dissection.

Next we want to compare the performance of the alternative eigenvalue decomposition. From here on, only double precision implementation will be considered.

One of the advantages searched in changing the decomposition is to make the method less sensitive to parameter variations, in this case, the value of ε > 0 used to identify linear dependency and also to diminish the number of saved directions without loosing in performance.

When using the Cholesky decomposition we can add an extra parameter that can recognizes when a gap between the decreasing pivots is too big. This gap is an indicator that we may enter into the limits of the noise directions, making this decomposition somehow more robust, as we avoid the noise at the cost of eliminating some other maybe useful directions. In any case this "lost" directions are usually recovered in later iterations.

For the Eigenvalue decomposition, we have also added an extra parameter to limit the number of directions saved. It works after the computation with max tolerance (given by the machine precision) of the eigenvalue algorithm. Then it limits the stored directions to the ones associated to the relative eigenvalues that are bigger than this new tolerance ε. This way of sorting directions, is a much more stable algorithm for sorting the directions than the previous Cholesky decomposition.

In the Table 2.5 we show the number of directions kept with each zero value until max precision for each method is achieved.

Here we can see that the bigger values of ε reduces the number of directions kept, but the convergence is still assured, as the important information is kept in the directions associated to the bigger eigenvalues. The eigenvalue decomposition in this case is less sensitive to variations in the zero parameters and in all the cases shows a reduction in the number of directions kept and total iteration number, meaning that at a small cost (in parallel computation it is negligible) we can have a method with less memory charge, less total iterations and improved precision.

We also insist in the fact that a much more precise method is obtained when using values of ε not so close to zero, just ε = 10 the minimal error and at the same time reduce the stored directions.

S-FETI with sparse storage

In this section we want to test the new sparse storage implementation for the S-FETI method, in terms of both time and memory usage.

The model problem is again the 3D Poisson equations, solved in the cube Ω = [0, 1] 3 , with the same boundary conditions and source as the previous examples.

We also use 125 subdomains with local interface subdivisions, in order to have at least the same amount of search directions as previous case. The parameter ν will represent the homogeneous case, meaning that we will have one material as shown in Figure 2.7. The time, memory and precision measures will be more important, rather than the number of iterations since this term should not change between the two sparse implementations and we hope to have similar results between the sparse and the full versions.

We start by testing the speed and precision achieved by the two versions of the SPARSE-S-FETI method, the first with a straitghtforward implementation when operating with the sparse blocks and the second one with optimizations in memory allocation to use the BLAS3 routines and reduce the bottleneck of doing several matrix-vector operations instead of just one matrix-matrix product. We The results in terms of memory and time are shown in Table 2.6, where we measure the total time used and also the maximum between the memory used by each process. We can see an speed up of roughly two times for both examples with 100 and 150 thousand elements approximately in each subdomain. We did not added the expected results in terms of number of iterations and precision, where there is no difference between both implementations. Since no precision is lost, plus the gain shown in the memory column and also in the time obtained, makes the second implementation to be the one used from here on.

Method

The objective of the sparse storage implementation is to keep a memory usage controlled in order to be able to use this method for the biggest applications. In this part we will compare the memory usage of the SPARSE-OPT implementation versus the regular FULL one.

In the Figure 2.10 we can see the maximum memory usage between process of both methods for the same problem, but with two different number of total elements in each one. Since we are interested in memory performance we will force to save all the search directions, with no regarding whatsoever of convergence, this is because we want an upper estimate in terms of memory versus the number of iterations. In this example we are using a number of 2072 directions in each iteration.

For the FULL-S-FETI method, we have a linear relation between the number of iterations (starting from 2) and the memory usage, or in this case the number of search directions stored, since every vector uses the same size. For the sparse storage, this iteration/memory relation is in general quadratic, however this relation is between the number of iterations and the size of the coefficient matrices stored, which are, as already explained dense but much smaller than a regular full block of search directions and is also independent of the size of the mesh (the size of the full vectors). In any case there is still an increasing linear relation between the sparse stored directions (Z p to recall) and the iterations, but as explained previously the size of them is less than the full block and it depends on the neighbour connectivity of each subdomain.

This results shows that the difference between the SPARSE storage and the FULL implementation is more clear as we augment the number of elements in each subdomain. In the first case, the FULL case is a better option, but as we get close to sizes of around 300 thousand elements in each subdomain the need of the sparse storage is more evident if we want to be able to solve this and other larger problems.

The following test to do, once the memory advantages of the sparse storage were established, is a time consumption test between the SPARSE-OPT and regular method to see the cost of all the extra computations.

The problem to solve is the same checkerboard configuration as previous with two materials, ν 1 = 10 0 and ν 2 = 10 4 . The results are shown in Table 2 where we see that the differences in precision and number of iterations are negligible. In the same way, the time difference is small considering the size of the problems, and it goes down whenever the problem gets bigger. We can say then that the new implementation is suitable to test big cases where the memory is a limitation.

General comparison

Finally we want to test the difference between the two existent preconditioners for S-FETI, the Dirichlet and Lumped. We test for memory, number of iterations and total time for the computations.

Added to this, we want to compare the sparse and full storage implementations with this preconditioners, as well as the two decomposition for sorting the search direction. The idea is to have global view that can gives some notion of which method should we use in different situations.

-Type of implementation: FULL-S-FETI or SPARSE-OPT-S-FETI -Preconditioner: Dirichlet or Lumped.

-Decomposition method for W T FW : Cholesky or Eigenvalue.

-Elimination of corner crosspoint interfaces: Yes or no.

We added the last point since in practical engineering applications, the corner interfaces are really rare and almost nonexistent. This is due to the fact that the domain division are done automatically by a mesh partitioning algorithm, such as METIS ( [START_REF] Karypis | METIS, a Software Package for Partitioning Unstructured Graphs and Computing Fill-Reduced Orderings of Sparse Matrices[END_REF]). This type of algorithm, with a regular cartesian mesh, does not create the corner interfaces that we usually see in square or cubes division of a domain. To see the difference in term of number of local solves and total search directions, we can "eliminate" this corners and test the method in a more realistic way, with non of the search directions created by this corner interfaces. In our particular case, we lower the number of total search direction per iteration from 2072 to 600.

One if the crucial difference that we want to test is the one produced by the different preconditioners available, in this case the consistent weighted Dirichlet and the, also consistent Lumped preconditioner. If we look at both definitions in Equation 2.57, Equation 2.58, and knowing that a multiplication by the Schur complement implies the resolution of a local problem, we can name two of the big differences between them, this is, the computational and memory cost, where for both characteristics the Dirichlet is much more expensive, however we expect a reduced number of iteration for this preconditioner, since it consider the variation of the complete subdomain and not only of the nodes in the interface.

The test case will be the exact same as previous, with 125 subdomains in a checkerboard configuration for the ν parameter and with 150 thousand elements per subdomain. In this case we will perform several computations and, as already said, we will measure the time and memory performance of them. Each result will be a different combination of four variables In the Table 2.8 we show the results obtained, we note that due to the characteristic of the problem, meaning a number of subdomain and size of the mesh not high enough, we obtain results that does not go with the lines of what to expect, as we see a reduction in time for the SPARSE-OPT, but with a higher memory usage. In Table 2.9 a remarkable speedup is obtained for this case when using the SPARSE implementation.

The elimination of corner interfaces also produces a reduction in the total time, making this results closer to a real implementation with automatic subdomain partitioning. In terms of memory usage, there is no big difference between methods, since the size of the problem is not big enough to make them clearer, as was the case presented before in this section. In any case we still recall that depending on memory availability the use of one or the other is recommended.

Conclusion

Several implementation changes have been introduced in this chapter, most of them showing an improved performance compared to the basic S-FETI. We added a new way of sorting search directions, the Eigenvalue decomposition can show to be a more robust option that reduces the number of directions stored in each iteration. In terms of memory usage we developed a new sparse storage that allows to use this method in bigger applications without losing significant performance. The use of the Dissection solver with his kernel detection algorithm are a good alternative for the largely used PARDISO solver. The Lumped preconditioner for this method, even if the results shown here are not definitive to say in which cases should be used, is still a good alternative to reduce memory and computation efforts, further testing will be done to this preconditioner as part of the future work. Finally, other changes and comparatives can still be performed, leaving space for improvements to this method and its implementations. This is all part of the future work for this method.

Chapter 3

Block strategies as a preconditioner

As a continuity on the development of the domain decomposition methods, in particular the FETI method, always searching for a fast and accurate algorithm that can be used in as much problems as possible, we will try now a different approach to the improvements already done to the S-FETI, using some of the ideas exposed in the previous chapter.

The main advantage shown by the S-FETI method is the construction of an enlarged search space in the Conjugate Gradient algorithm where the capture of the largest eigenvalues of the FETI operator, even for some hard problems, can be achieved in very few iterations. The price to paid for this is the use of large amounts of memory that can render impractical its use for some big problems, usually found in real applications, in this context we already developed a way to circumvent this limitation by using a sparse storage of the directions.

In this chapter, we will try to mimic the precision and memory use of the original FETI method, improving the convergence ratio by building a new preconditioner based on the enlarged space generated by the S-FETI method. This new preconditioner is a priori more expensive, but hopefully the faster convergence expected will compensate this and we may have a useful method for practical use. This chapter starts with some preliminaries to understand the idea of the new preconditioned FETI method, later we will show a natural extension of this idea that mixes it with the S-FETI method to end up with numerical results of 157 both new methods.

Introduction and preliminaries

The FETI method is based in finding the solution of an interface problem using the Preconditioned Conjugate Gradient Method. This iterative method (CG) has been largely studied over the years, leading to different results in its research from which several generalizations were made to it.

Within this generalizations we can name the Flexible Conjugate Gradient, the Multipreconditioned CG and the Method of Conjugate Directions. Although the last one was developed before the CG method it can still be considered as one of the generalization of it.

We will now show in detail, two of this methods and their relation of some of them with previous FETI methods. This relations will allow to understand the roots of the new FETI methods proposed in this chapter and justify its convergence.

Method of Conjugate Directions

We consider the solution of the general system

Fλ = d (3.1)
with F symmetric positive definite.

This method, shown in [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF] can be considered as a more general method than the CG method.

In this method the search directions Choose a direction w p+1 such that (w p+1 , Fw j ) = 0, ∀j = 0, . . . , p 10: end loop is not precise, in the sense that no formula is given for choosing the search directions. Depending on how to build this directions, different methods can be obtained (including the CG method). Later we will explicit this definition in order to formulate the new FETI method.

w 1 , w 2 , . . . ( 3
The basic properties of this method are given by the following theorems Theorem 3.1. The search directions w 1 , w 2 , . . . are mutually conjugate. The residual vector g p is orthogonal to w 1 , w 2 , . . . , w p-1 . The inner product of w p with the residuals g 0 , . . . , g p is the same. All this can be written as

(w i , Fw j ) = 0, i j (3.4) (w i , g p ) = 0, i = 0, 1, . . . p -1 (3.5) (w p , g 0 ) = (w p , g 1 ) = • • • = (w p , g p ) (3.6)
The scalar ρ p can be given by the formula

ρ p = - (w p , g 0 ) (w p , Fw p ) (3.7)
Computing the inner product

(w i , g 0 ) = -α 0 (w i , Fw 0 ) -• • • -α p-1 (w i , Fw p-1 ) = -α i (w i , Fw i ) (3.15)
an using 3.7 we imply that α i = ρ i and hence λ = λ m .

Finally, we can see that the approximated solution improves with the number of iterations and the error is minimized at each step. Theorem 3.3. Let λ be the solution of 3.1 and f be the error function

f (v) = λ -v F = (λ -v, F(λ -v)) = (v, Fv) -2(v, d) + (λ, d) (3.16)
then, the point λ p minimizes f (v) on the line v = λ p-1 + αw p-1 . In fact, the point λ p minimizes f (v) on the space λ 0 + span{w 0 , . . . , w p-1 } (3.17)

and this space contains the vectors, λ 0 , . . . , λ p

Proof. If w is a search direction, then we have

f (v + αw) = f (v) -2α(w, g) + α 2 (v, Fv) (3.18)
where

g = d -Fv = F(λ -v) (3.19)
Considered as a function of α the function f (v + αw) has a minimum value at

α = ρ, where ρ = (w, g) (w, Fw) (3.20)
Comparing this, with the definition of ρ p in the algorithm, we have the minimization result.

For the last part, lets consider the point v ∈ λ + span{w O , . . . , w p-1 }, then we have

v = λ 0 + α 0 w 0 + • • • + α p-1 w p-1 (3.21) applying f f (v) = f (λ 0 ) - p-1 i=0 [2α i (w i , g 0 ) -α 2 i (w i , Fw i )] (3.22)
At its minimum we have

α i = (w i , g 0 ) (w i , Fw i ) (3.23)
and so α i = -ρ i , for all i = 0, . . . , p -1. Then is clear that the minimum point is λ p .

Considering the previous theorems, as long as we do the full reconjugation, we have the liberty to build the search directions in any form we want.

Flexible Conjugate Gradient

As a continuation of the previous Conjugate Directions method, and within the Conjugate Gradient framework, we can name the General CG or Flexible CG (as denoted in [START_REF] Notay | Flexible Conjugate Gradients[END_REF]). We consider again the problem 3.1, solved with the Preconditioned CG method, the idea is to add at each iteration of the CG method a preconditioner step with the aim of accelerate the convergence by using suitable search directions, usually this is achieve by finding the solution of a problem of the type Bx = y (3.24) so that the conditioning of the matrix B -1 F is smaller than the conditioning of F.

In PCG the preconditioner B -1 is symmetric positive-definite and fixed, in order to keep the convergence properties of CG. In the Flexible-CG we relax this conditions and change the preconditioning step to a general mapping .25) this mapping also tries to approximate the inverse of F in order to improve the converge. This more general mapping does not need to have and explicit form as it can be a recursive non-linear mapping. We will simply consider it as a general procedure to approximate the inverse of F.

g → B[g] ( 3
With this preconditioning procedure and using a full reorthogonalization process we can define the algorithm for the Flexible CG method in 12. This for i = 0 to p do

11: γ = - (w i , Fw p+1 ) (w i , Fw i ) 12:
w p+1 = w p+1 + γw i 13:

end for 14: end loop algorithm explicits the construction of the F-conjugate search directions by introducing a preconditioning procedure and applying the Modified Gram-Schmidt procedure in order to build the new directions.

It is also essentially the same CG algorithm used by FETI, the main difference lies in the preconditioner. If we consider B as the Dirichlet or Lumped preconditioner, we can recover the usual FETI method.

Remark: The detailed convergence properties of this algorithm, including bounds for the error, depend on the definition of the mapping B, were conditions on coercivity and boundness are assumed [START_REF] Axelsson | Generalized Conjugate Gradient Methods[END_REF]. In our case, and in order to keep the mapping as general as possible, we will only consider the previously shown convergence of the CD-method to validate our new FETI method with recursive preconditioner.

FETI with recursive preconditioner

The Simultaneous-FETI method can be considered as a Multipreconditioned CG method [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF], this is because the fact that several directions are created in the preconditioning step, thanks to the additive nature of the Dirichlet preconditioner in FETI. Closely related to this method we have the Flexible CG in which recursive or non-linear preconditioners can be used at the cost of doing a full reorthogonalization of the search directions (truncated orthogonalization are also possibles). In the FETI framework there are no practical problems in the use of this extra computations, since we are forced to do them in order to keep a good convergence (See [START_REF] Roux | Acceleration of the Outer Conjugate Gradient by Reorthogonalization for a Domain Decomposition Method for Structural Analysis Problems[END_REF]).

Following this idea, and since we already do a full reorthogonalization, we can consider the general method of Conjugate Directions which give us total liberty in the creation of the search directions, as long as they are conjugate between them, more precisely F-conjugate (F being the FETI operator, which is symmetric positive definite). This last part made us think in consider a sort of mixing of this methods, where we consider some special preconditioner that changes from iteration to iteration, but it is a much better preconditioner than the usual Dirichlet in the FETI method, since it is built straightforward from the one in S-FETI where clearly a better approximation is obtained at every iteration. Even if we are not exactly in the conditions to consider this new method as a Flexible CG method, we can see it as a Conjugate Directions method, where the convergence is assured. Now we will show the details of this new method, and later its extension to an "in between" method, if we consider the S-FETI method and this new one.

One direction from S-FETI block

We derive this new method starting with the FETI basic algorithm, and modifying the preconditioning step w p = F-1 g p (3.26) with F-1 being any suited preconditioner (mainly the Dirichlet or Lumped one). Even though the multiplication for this precondition matrix is not computed in a direct way, we can still consider this step as simple and not so expensive (details about this have been discussed in previous chapters), the main difference will be in this step, where a more elaborated preconditioner will be constructed. We change this multiplication by the Dirichlet or Lumped operator for the definition of a new procedure (or mapping) B based on the block construction done in S-FETI.

The idea behind this is to use of the minimization process of the CG method applied in the S-FETI method, where we construct a bigger search space, based in the combination of local terms. Lets recall this process in Algorithm 13 where after the computation of the search directions block, we obtain the optimization parameters and update the solution In this new method we will consider Algorithm 13 Minimization process in S-FETI 1:

. . . . . . the update vector W p ρ p as our unique search direction, for it contains all the information from the non local preconditioner, and it can help speed the process of capture of the eigenvalues. Meaning that at every iteration of the Conjugate Gradient method, we will define the search direction as

w p = W p ρ p = -W p (W T p g p ) (3.27)
The construction of W p was already explained when we define the S-FETI method in chapter 2 and it involves the solution of several Neumann and Dirichlet problem (if Dirichlet preconditioner is used), so it is more costly than the regular preconditioner. This extra cost compared to the gain in the convergence will be tested in the numerical examples.

Adding the construction of the block W p to our new search direction allows to explicit the definition of the procedure B[g] used as a preconditioner. This mapping is described in Algorithm 14. The notation in this procedure is Basically what we are doing is similar to what we do in any preconditioned CG algorithm, that is, a process from which we take the gradient and we apply a certain operator to build a search direction that approximates to the actual error F -1 g p , at the p iteration. The difference in this case is that the search direction does not come from a direct multiplication for a fixed matrix, but from a more complicated process that in the end makes the new direction to depend on the previous gradient.

The final method is shown in Algorithm 15. What we describe is a modification of the FETI method (using the Flexible CG instead of classical CG), in which we change the computation of the preconditioner with the procedure coming from S-FETI that allows to define B[•] in Algorithm 14. The definition of the projection operator P , the matrix A, the matrices G of the FETI coarse space and the block of vectors of the kernel R are defined in the same way as previous methods (FETI-1LM or S-FETI).

We hope that with this new method we are going to recover the accuracy and memory usage of the basic FETI method, by only storing one direction per iteration and at the same time diminishing the roundoff effects of the extra sources of numerical noise that can be found in the S-FETI method. As for the advantages, the number of iterations will be reduced, due to the improved way to for i = 0 to p do 11:

γ i = - (w i , Fw p+1 ) (w i , Fw i ) 12:
w p+1 = w p+1 + γ i w i

13:

Fw p+1 = Fw p+1 + γ i Fw i 14:

end for 15: end loop capture the largest eigenvalues, and we expect this will be enough to compensate for the extra cost on the mapping of the preconditioner.

Cost of the method

Considering the cost of the original FETI method and the one of the S-FETI method, we can give a comparative estimate of the performance of this method against the original FETI at each iteration.

We know by its definition that all we need to compare is in the preconditioning step in FETI versus the preconditioning procedure defined in Algorithm 14. We note that this is also because the procedure B[•] will also give the vector Fw p at the end of the procedure since its value correspond to FW ρ, hence no extra multiplication by F are needed after this part.

Based on the analysis of the cost in S-FETI we can say that the most expensive part here will also be the computation of the block FW . This includes a projection and several forward-backward computations, but we can use the same strategies to minimizes the computations as done in S-FETI, meaning that for each iteration we only add to a process the workload of N Neumann solves, with N being either neighbours + 1 or 2 × neighbours depending on the construction of the blocks. The projection is also done using low-rank corrections. For details in the cost of S-FETI and thus this preconditioning procedure, see the section 2.1.5 in the previous chapter.

Finally we can say that the cost of each iteration is approximately an iteration of S-FETI (which is at the same time have a few extra cost compared to FETI) plus an iteration of FETI without the Neumann solution, but since we can not estimate the number of iteration needed for convergence (hopefully much less than FETI), we will analyse in the numerical examples the real performance of this method and test different ways to compute the block W used in the mapping B[•], namely the diagonalization process in the treatment of the term W T FW and the use of the Lumped preconditioner.

Linear combination of directions from block

In this section we will extend the previous idea of using the added combination W ρ built in the S-FETI method as a single direction, to an "in between" method, where we will save in addition to the unique previous directions a reduced number of vectors containing the information of the block of search directions generated by S-FETI. This number will be a limited number of vectors between 1 and the total number of linear independent search directions in the iteration N S .

In order to achieve this, we will proceed in a similar way as before, only this time, instead of using the FETI method as a base we will use the S-FETI method and we will change the directions stored after the update of the gradient at the end of each iteration. First we will recall the construction of the block of directions W in S-FETI, but this time computed with the eigenvalues decomposition, since this will allow to sort the total independent search directions in decreasing order of importance. In this case, starting from the gradient we have that Algorithm 16 Block construction in S-FETI 1:

. . . . . . From here we know that the matrix W p is formed by columns of normalized vectors that are sorted in decreasing order from the ones that capture the bigger eigenvalues of the FETI operator to the ones that represent the smaller ones.

Is from this characteristic that we can use and save the more relevant directions. In the previous case we added each column to build a unique search direction

w p = W p ρ =            W 1 p . . . W N S p                        ρ 1 . . . ρ N S             = i ρ i W i p (3.28)
with ρ the vector of minimization parameters ρ = -W p g p . This time we split this sum and we will recover a new block of search directions, smaller than the previous one, and to try and not to lose information we add an extra direction either the vector W p ρ, this implies the construction of yet another block Wp defined by

Wp =            W p ρ ρ 1 W 1 p . . . ρ M W M p            (3.29)
Where M is a number that depends on the number of columns, fixed by the user, that we will study in detail in the numerical tests. What we try to accomplish here is the reduction of the relevant information to a limited number of search directions, improving convergence but with controlled memory usage, at expenses of the extra computations needed to build the original W p . As previous method, the definition of the FETI projector P is as usual. In Algorithm 17 we described the final algorithm, build from the basic S-FETI, but changing the construction of the stored block directions. The definition of the projection and all the matrices to compute λ 0 are defined as in previous FETI methods ( see 1. for i = 0 to p do 8:

Γ i = -W T i FW 9: W = W + Γ i W i 10: FW = FW + Γ i FW i 11:
end for 12:

EDE T = W T FW 13: W = W ED -1/2 14: ρ = -W T g p 15:
λ p+1 = λ p + W ρ 16:

g p+1 = g p + FW ρ 17: W p = W ρ, ρ 1 W 1 , . . . , ρ N W M 18:
FW p = FW ρ, ρ 1 FW 1 , . . . , ρ N FW M 19: end loop Since the base of this algorithm is S-FETI, its convergence is assured. This fact comes from the Theorem 2.1 where the result is independent on the preconditioner or the construction of the block of search directions, as long as they are F-conjugate, which clearly is our case.

Cost of the method

The cost of this algorithm, at every step, is derived directly from the cost of S-FETI in subsection 2.1.5, since we do the same construction to obtain the block of search directions. The only difference is the size of the stored vector that may be reduced as much as we want, however in the next section we will try to find the "optimal" value of needed stored directions.

Numerical results

Since the methods shown in this chapter can be considered as variations of the S-FETI algorithm, the characteristics of them in terms of problems that can or cannot solve are very similar to the original S-FETI. Is for this reason that a simple test case can and will be considered next, in order to make this first comparative between this methods.

The Poisson problem, or any other simple case, solved in 3D will allow us to do this first comparative of both method presented in this chapter.

We recall the variational formulation of this problem for some parameter ν ∈ R + and some bounded Ω ∈ R 3 .

Find u ∈ V such that ν Ω ∇u • ∇vdx = Ω f v, ∀v ∈ V (3.30)
where

V := v ∈ H 1 (Ω) : v| ∂Ω D = 0 (3.31)
and ∂Ω D is the section of the boundary where we impose u = 0.

A finite element discretization, using tri-linear Q 1 functions is used to approximate the solution of this problem.

We finish the problem by setting the characteristics of our model, which in this case are a cube domain Ω = [0, 1] 3 with boundary conditions in ∂Ω D = {x = 0} ∪ {x = 1}, the value of the source is f = 1 and the parameter ν will have an homogeneous value along the domain. The cube will be divided in 125 smaller cubes of equal size, with variable number of elements in each one, and the stopping criterion is the same as previous chapters.

Storage of single direction

The first method to test will be the one where all the search directions are reduced to a single one. In simpler terms, we reduce the S-FETI approach to a FETI method with a more expensive preconditioner. In terms of memory, we recover the same use as the FETI-1LM original, so no test related to memory will be performed.

Method

Our first and main approach will be to asset the number of iterations needed to converge for this method and compare it to the original FETI-1LM. Using the already described configuration, with ν = 10 -3 . we can see the results in Table 3.1.

In terms of iteration numbers we can see a reduction in the total number of them, compared to FETI-1LM, but no versus S-FETI. This goes with the lines of what we expected, meaning that an expensive preconditioner can improve the FETI-1LM method, but the information lost when storing the constructed directions as a single vector will reduce the convergence compared to S-FETI.

In any case, the excessive cost of this approach does not compensate for the good results shown when reducing iterations of FETI-1LM.

Storage of reduced directions

We continue with the extension of previous method, so in this case, we change the reduction of the search directions from a single vector to a small number of vectors. Since this vectors are computed in descending order of importance, we want to keep a minimal number of search directions without loss of convergence.

We continue to use the same test case, that is, the Poisson problem in the cube with 125 subdomains and fixed 75 thousand elements per subdomain. This is complex enough, mainly due to the number if directions produced in each iterations, so we can expose the differences between this method and the original S-FETI and FETI.

In the definition of this method, we have to choose a fixed number of search directions to keep at every iterations. In that sense a fixed number M was introduced, but what we actually do in practice is to save a percentage of the independent search directions. This value has a great impact in the method, so we need to do variation of it to try in find an "optimal" value if possible.

In Figure 3.1 we show the behaviour of the method in terms of total iterations versus the percentage of directions stored. Lets recall that the total number of search directions generated, without consider linear independence is 2072, from which we do a first reduction in number, to later keep a part of this reduced number.

In this image we see the behaviour of the method, with an initial fast reduction in the number of iterations, even if we store only a small percentage of directions, such as 10%. Values close to this can be considered as "optimal" for this method if we are looking for memory reduction and higher precision for the method. Depending on the problem to solve and the machine used, we choose the best suited value for every case.

Conclusion

The variations of the S-FETI method shown in this chapter, are a priori ideas worth testing, since the storage of the information produced by S-FETI can be an issue. Numerical results shows that in the actual conditions of the first variation, even if we obtain and improvement in terms of iterations to convergence, precision and memory, there is no gain total time due to the too much expensive cost of the implementations used. The second method arises as an alternative to certain particular cases, where at the price of more iterations we reduce the memory usage. In any case, we find important to state this results, to avoid the same work in them and at the same time, it can still be source of new development for the S-FETI and other FETI like domain decomposition methods.

Chapter 4

FETI-2LM with enlarged search space

In this final chapter we will now show the extension of the main idea in the S-FETI method, that is, enriching the search space of the FETI algorithm, but this time to one of the other FETI method shown in this thesis, in particular the FETI-2LM method.

The FETI-1LM and FETI-2LM methods share similarities in the iterative algorithm used, both of them are Krylov based, the Conjugate Gradient and the ORTHODIR method. This fact makes the idea of enriching the search space in FETI-2LM, a priori, a good one, but big part of the convergence properties of S-FETI comes from splitting the preconditioner to isolate the modes that make the convergence slower, meaning that the fast convergence can be expected, also because the relation of this method with the FETI-Geneo method which has a proven fast convergence speed.

In the FETI-2LM method we no longer have an optimal preconditioner, nor a method to compare, so we do not expect to have the same speed up as in S-FETI, but we will have in any case an large improvement in the convergence, due to the constructions of directions using local properties that tend to get lost in the FETI-2LM method, but also we will gain from the size of the new enriched space compared to the original one. Finally in terms of time spend, we will be applying the same implementation strategies as in S-FETI, hence a time speedup is most assured.

We start by recalling the FETI-2LM method, to then pass directly to the

Ω (2) Γ (3)
Ω (1) λ (2) λ (1) Robin Condition

Introduction

Most of the elements needed in this chapter, were already described throughout the previous chapters, but in any case we want to recapitulate some of the more important ones, in order to keep clarity in the new method to be constructed.

The FETI-2LM method

This method started as a solver for acoustic problems, first shown in [START_REF] Farhat | Two-level domain decomposition methods with Lagrange multipliers for the fast iterative solution of acoustic scattering problems[END_REF], then it was generalized in [START_REF] Roux | Optimization of Interface Operator Based on Algebraic Approach[END_REF] and [START_REF] Roux | Approximation of Optimal Interface Boundary Conditions for Two-Lagrange Multiplier FETI Method[END_REF] as a robust solver for problems coming from any Finite Element discretization of an elliptical partial differential equation.

It is based on the imposition of a Robin condition on the interface, see Figure 4.1, in order to "glue" the solutions of each subdomain. We will recall the basic formulation for a two-subdomain division, as the general formulation comes from using the two-subdomain case in each pair of neighbour subdomains.

Consider the linear problem, arising from a Finite Element discretization of a PDE

Kx = f (4.1)
We divide the global problem in two subdomains Ω (1) , Ω (2) and their interface Γ (3) . Then the global problem has the structure

            K (1) ii 0 K (1) ib 0 K (2) ii K (2) ib K (1) bi K (2) bi K (1) bb + K (2) bb                         x (1) i x (2) i x b             =             f (1) i f (2) i f (1) b + f (2) b             (4.2)
and the subdomain stiffness matrices and right hand sides are:

K (1) =        K (1) ii K (1) ib K (1) bi K (1) bb        , f (1) =        f (1) i f (1) b        K (2) =        K (2) ii K (2) ib K (2) bi K (2) bb        , f (2) =        f (2) i f (2) b        (4.
3) The FETI-2LM method comes from the imposition of independent generalized Robin condition on the interface Γ (3) .

The discrete local problem are then

       K (s) ii K (s) ib K (s) bi K (s) bb + A (s) bb               x (s) i x (s) b        =        f (s) i f (s) b + λ (s) b        , s = 1, 2 (4.4)
The augmentation matrix A (s) bb will be considered as an sparse approximation of the neighbour Schur complement, as seen in subsection 1.2.3.

In order to be the restrictions of the global problem, each local solution must satisfy two conditions, first the continuity condition: The definition and implementation of this method for general problems, comes from the computation of this operator between every two neighbour subdomains in a general configuration.

The previous condensed non symmetric problem is solved using the OR-THODIR iterative method. This is because in the general case, the operators is dense so its assembling is computationally inefficient. On the contrary, in order to multiply some vector by this operator (denoted F as usual FETI methods) we only need local data and exchanges between neighbours, hence an iterative method such as ORTHODIR is more suited.

We present in 18 the ORTHODIR algorithm used in FETI-2LM, as it will be the base to construct its generalized block version. for i = 0 to p do 12:

γ i = - (Fw i ) T (Fw p+1 ) (Fw i ) T (Fw i ) 13:
w p+1 = w p+1 + γ i w i

14:

Fw p+1 = Fw p+1 + γ i Fw i 15:

end for 16:

end loop 17: end loop

In the next section we will extend this algorithm to its block version using the decomposition of the gradient to build the block of search directions used to update the solution in every iteration.

The Block-2LM Algorithm

Following the work done in the development of the S-FETI method, we can apply the same strategy to the FETI-2LM method. In this case, the method called Block-2LM, will build a search space based on the separation of the gradient, with the difference that the gradient this time will no longer be preconditioned, but it can also be considered as a sum of local gradients that when separated will improve the mimic of the local subdomains behaviour. In any case, since in FETI-2LM we have two independent lagrange multipliers instead of one in each the interface, we can create as much directions as the S-FETI method in both of its version (subdomain or interface divisions).

Following the notation of chapter 2, the gradient g ∈ R n in the FETI-2LM method is

g =            
g (1) . . .

g (N s )             (4.12)
with N s being the number of subdomains. Also g (i) ∈ R m i , m i ≤ n are the values of the global gradient in the subdomain Ω (s) , this local gradients are independent from each other since there are two lagrange multipliers on the interface. This subdivision is enough to form a new block of directions, all by considering the fact that g can be written as

g = N s s=1 ḡ(s) (4.13)
where ḡ(s) ∈ R n is the extension by zero of each local vector, i.e.

ḡ(s) =                                    0 . . . 0 g (s) 0 . . . 0                                    (4.14)
with this division, our block of search directions will be defined by Z = ḡ(1) , ḡ(2) . . . , ḡ(N s ) (4.15)

Which is similar to the subdomain division in S-FETI.

In an analogous way, to form the interface subdivision, we can consider at the local level that g (s) can also be decomposed in smaller vectors, one for each of the subdomain neighbours with n (s) the number of neighbours of subdomain Ω (s) . Hence, the general gradient can also be written as g = N s s=1 n (s) i=1 g (s),i (4.17)

g (s) =             
where g (s),i ∈ R n is again an extension by zero, defined as g (s),i := From the equation 4.17, we can use the additive form of the gradient in order to define the new search space for the Block-2LM method. This space will be spanned at each iteration by the column vectors in Z = g (1),1 , . . . , g (1),n (1) , . . . , g (N s ),1 , . . . , g (N s ),n (N s ) (4.19)

                                   0 . . .
The iterative algorithm used for the FETI-2LM method is the ORTHODIR method, so the Krylov space built is formed from the successive matrix gradient product, orthonormalized using the F T F-norm. We are now spanning a different search space, but the use of this bigger space still implies the convergence of the method since it contains the original Krylov space, this result comes from the definition of Z and the fact that Span{g (1) + • • • + g (m) } ⊆ Span{g (1) , . . . , g (m) } (4.20) for any vectors g (i) ∈ R m , m ∈ N, i = 1, . . . , m.

Since the norm used in this method changes, from F to F T F-norm, the matrix that needs to be inverted also changes. This matrix is used to obtain the block of orthonormal vectors and also to obtain the optimal descent coefficients. In this case this matrix correspond to (FW ) T FW , where W is the block after applying the full reorthogonalization, i.e. at the p-iteration we have

W p = Z p + p i=0 W i Φ i (4.21)
The same linear dependency issue occurs for this matrix as in the S-FETI 2.2. Lets recall that this is due to the phenomenon of work within the computer precision or because we achieve convergence in some of the interfaces before the other. This dependency between columns in Z leads to a positive semi-definite (FW ) T (FW ) matrix, so as previous treatment, the Cholesky decomposition with partial pivoting will be used to select the directions to be used (The Eigenvalues decomposition is also valid). The algorithm is described in 19.

Implementation and cost of the method

The difference in cost between the FETI-2LM and Block-2LM are similar as the one from FETI and S-FETI, we can name the augmentation on the size of the information shared between neighbours, but with no incremental number of exchanges, and we also have the pseudo-inversion of the matrix (FW ) T (FW ), both of this changes are fairly cheap ones.

The most expensive part is again the multiplication by the FETI operator, this can be done with a simple straightforward implementation that solves a local for i = 0 to p do 16:

Φ i = -(FW i ) T (FZ p+1 ) 17:

W p+1 = Z p+1 + W i Φ i 18:
FW p+1 = FZ p+1 + FW i Φ i 19:

end for 20:

end loop 21: end loop problem for each column in Z, this can only serve for academical purposes, to know the precision and convergence behaviour of it, but if we want a practical method we need again to exploit the sparse structure of the search directions matrix Z.

We know from the definition that the columns g (s),i of Z are non null only on the interface of subdomain Ω (s) and the multiplication by the FETI operator only requires, at the subdomain level, the solution of a number of local problems equal to the number of neighbours, which can be done in a block implementation. With this consideration we can have the speed up to make this method useful in real problems, since the time spend during the resolutions we hope to recover it by reducing the iterations needed for convergence.

To understand the implementation of this method, we will again use the so called Coarse Modes, that correspond to the non-null vectors in each column of Z or FZ. This modes are stored locally in each subdomain, in this case we have one mode for each neighbour, this correspond to each of the vectors g (s),i loc defined in Equation 4.16, but at the same time we shared them with every neighbour, for further computations, e.g. multiplication by the FETI operator.

In Figure 4.2 we can see in the left, the described modes owned by the subdomains, where in dotted lines are the ones own by a subdomain, also we see the subdomains involved in the computing of FZ, and in the right, the dotted modes are the ones owned by the red interface, at the same time this modes are non null modes of FZ in the red interface, this comes from the definition of the operator.

To better understand this, lets look at the definition of the operator F for a single interface

F        λ (s) 1 λ (q) 1        =        λ (s) 1 + λ (q)
1 -(A (s) + A (q) )(S (q) + A (q) ) -1 λ (q) 1 λ (s) 1 + λ (q) 1 -(A (s) + A (q) )(S (s) + A (s) ) -1 λ 1 -(A (s) + A (q) )(S (s) + A (s) ) -1 λ (s) 1

       F       0 λ (q) 1       =        λ (q)
1 -(A (s) + A (q) )(S (q) + A (q) ) -1 λ

(q) 1 λ (q) 1        (4.23)
We note that each subdomain solves one local problem for each local interface (number of neighbours). All the local solutions plus the multiplication by local augmentation matrix are then send together to each of the neighbours, meaning we send the vectors (S (s) + A (s) ) -1 λ (s) 1 A (s) (S (s) + A (s) ) -1 λ (s) 1 (4.24) this for every local computation performed, allowing to multiply this by the local augmentation matrix associated to the interface and in general compute and store a local contribution to FZ. With this, each local interface stores the contributions to FZ done by the neighbour subdomain.

Finally, in terms of memory use, the same issue described in chapter 2 is expected with this implementation, but we can again use the same strategy to store the sparse directions and reconstruct them later, although we will leave that as a future work.

Remark:

We want to point out some of the differences between this method and S-FETI, mainly the fact that no projection is needed in this case, so no extra cost is added in this part and also the fact that the connectivity is more limited than in the S-FETI method, thanks to the FETI-2LM operator, hence we need less modes stored in each interface.

Numerical results

In this section, we will be performing another basic comparative testing, that will include, mainly the FETI-2LM and Block-2LM methods. Classical comparison between iteration number, as well as a full fortran-mpi implementation to measure for time performance will be analyzed.

To perform this benchmark, we do no need for a special complex case, as one of the main properties of both method is the robustness and is expected to be similar. We will leave the tests for complex cases as future work. From now on, we will focus on the study of the basic Poisson problem in 3D, with tri-linear Q 1 functions for the finite element discretization. For more details, see the numerical results in chapter 3 and chapter 2. We just recall the existence of the ν parameter, that also in this case will change between neighbour subdomains, meaning that an heterogeneous checkerboard type of configuration will be tested, and blocks of two different materials will be used.

Block-2LM vs 2LM

The global domain is again the cube Ω = [0, 1] 3 , several subdomain divisions will be done, all characterized by being smaller cubes of equal size. We want to This effects is less present in the S-FETI methods, because the computation in this part can be reduced by 2 since the values of the Lagrange multiplier are shared by neighbour subdomains, but in this case we are forced to work with full size blocks of search directions since two independent Lagrange multipliers are used in this case. Also the number of iterations goes up as high as 62 and the quadratic effects of full reorthogonalization produces the notorious slowdown.

We consider this implementation as a first step before doing the sparse storage implementation, the same way that was done in chapter 2. In that case the performance of this method, specially when doing the full reconjugation, should be several times faster, and hence a method to be used in practical applications.

Remark: As part of the future work, we leave the test of a variation of this method, with even more directions per iterations, the idea is to use the neighbour local gradient as a new direction doubling the number of directions built. The normal block of directions, using the notation of this chapter, is defined as Z = . . . , g (s),1 , . . . , g (s),n (s) , . . . , s = 1, . . . , N s (4.25) and when doubling the direction we will have Z = . . . , g (s),1 , g vec{(s),1} , . . . , g (s),n (s) , g vec{(s),n (s) } , . . . , s = 1, . . . , N s (4.26) where g vec{(s),i} is the extension by zero (using the same position as the non zero values of g (s),i ) of the local gradient coming from the neighbour subdomain of Ω (s) through his local interface i.

Conclusion

A new Block method as been developed for the FETI-2LM following the ideas used in S-FETI. Although the properties of FETI-2LM, such as the lack of a preconditioner or the fact that two lagrange multipliers are already imposed in the interfaces, make this method a less appropriate candidate for a block version than the FETI-1LM was. However the fact of enlarging the search space, still predicts a better convergence in terms of iterations for this method. The numerical results validate the reduction of iterations needed by this method to converge. Nevertheless this reduction is not compensated as in the S-FETI case, mainly because the total number of iterations achieved is not sufficiently small to keep the time spend in the full reconjugation as a factor controlled. Making this factor the major slowdown for the Block-2LM method.

Alternatives to the reconjugation must be found, mainly the sparse storage proposed in chapter 2, but we leave that as a future work.

Conclusion and perspectives

Throughout the work done in the four chapters exposed in this thesis, we have detailed several improvement to existing classic FETI method found in the literature that are currently used in real life applications that comes from the industry, namely structure, electromagnetism models among others.

The first part consisted in the development of the Hybrid-FETI method, born from the mixing of the FETI-1LM and FETI-2LM methods. This methods consist in the imposition of either Neumann or Robin conditions, as done in FETI-1LM and FETI-2LM respectively, in the local interfaces as way of "gluing" them.

Due to the particularities of both base methods, the Hybrid-FETI is well suited for solving problems where a small number of interfaces present a bad conditions beyond the heterogeneities. This is the case of contact problems, where the contact area presents an ill-conditioning which the FETI-1LM even with the consistent version of the Dirichlet method does not assure the convergence. This problems are usually solved using the FETI-2LM due to the formulation with two independent lagrange multipliers provided by this method, that handles the issue of the contact area. Since the use of the FETI-2LM is only required for the contact area, we can improve the global convergence by treating the non contact interfaces with the regular and faster FETI method. This global liberty of choice for the interfaces is our new Hybrid-FETI method.

When imposing one or the other type of interfaces, we form subdomains where both methods are active, forcing us to create a preconditioner optimal for this subdomains. In the subdomains where only the Neumann boundary condition is imposed, we use the regular Dirichlet preconditioner to boost the global convergence. In this context, as an extension of the regular Dirichlet preconditioner, the Robin preconditioner was developed. This preconditioner differs from the Dirichlet only in the subdomains with mixed interfaces, in this subdomains the augmentation matrix present in the local problem is added to the Schur complement usually computed in the Dirichlet preconditioner.

Although we see and improvement with the Robin preconditioner, we leave as future work for this part the development of an additional preconditioner for the 2LM interfaces where no preconditioner is actually applied or in general a new preconditioner for the global operator of the Hybrid-FETI.

Looking at numerical results, the best one obtained by this new method are seen in the contact case where we reduced by almost two times the number of iterations of the only working base method, which is the FETI-2LM. This reductions should be directly applied to a reduction in terms of total time, but since we did not have the chance to test bigger cases, this will rest as a developing work. In this same line, we want to extend the application of the Hybrid-FETI to other problems solved only by the FETI-2LM method, for example the ones arising from the use of mortar elements to treat non-conforming meshes, where also an ill-conditioning is present in a localized zone, see [START_REF] Roux | A FETI-2LM Method for Non-Matching Grids[END_REF].

In the second part, we present the development of the S-FETI starting from his first version presented by [START_REF] Rixen | Substructuring and dual methods in structural analysis[END_REF] to his more precise general formulation, shown in [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF]. The development of the consistent Dirichlet preconditioner is also exposed, as this is a key part of the good behaviour of S-FETI. Details of the implementation were also explained, also needed to understand the new sparse storage developed in this chapter.

Following the presentation of the method, we introduced a variation to the original method. This change comes from the problem of treating the linear dependency usually present when constructing a block of search directions. This dependency comes from different elements, that are present in almost all problems. A good recognition of them can have a large impact in the final results. The Eigenvalue decomposition is proposed as an alternative to the basic Cholesky factorization that uses the basic S-FETI, the idea of this variation is to keep the good behaviour but using less directions in each iteration.

Confirmed by the numerical results we can see that the same convergence can be achieved, but using less directions in each iteration. This reduction comes from the fact that the Eigenvalue decomposition can order the orthogonal search directions in order of relevance an reduce the complete information contained in a block of directions, on the contrary the Cholesky decomposition only eliminates the ones considered problematic, a process that leads to lost of information at each iteration.

One of the goals of the Eigenvalue decomposition is try to reduce the number of directions used in each iteration, and therefore make the usage of memory S-FETI method less expensive. Even though there is an improvement present in here, the memory cost of this method makes it limited if we want to solve bigger cases, to solve this, a new sparse storage was developed. This implementation is based on the sparsity present in the block of search directions, where at a local level each block of directions contains non null columns in the ones representing the same subdomain and its neighbours. Using this information we can reconstruct the directions from a series of smaller coefficient matrices and at the same time we can reduce the multiplication of the sparse blocks by the FETI operator to faster matrix-matrix computations.

In the numerical results we show that the impact in total time done by the extra computations is negligible if we do a correct implementation. Several improvements to the basic implementation were introduced, leading to a new S-FETI with sparse storage as fast as the original one, however there is still room for new improvements as we have augmented the complexity of the algorithm by introducing several new sources of computations, that can be also improved. We leave that as future work.

The third chapter is dedicated to test different ideas regarding variations of the S-FETI, we can name the two main ones, which consist on one side of using the block of search directions build by S-FETI as a preconditioner for the original FETI method, and the second idea is to add to the search direction created in the first method a small part of the directions associated to the biggest eigenvalues when doing the decomposition of W t FW .

Both of the method show an improvement in term of iterations needed for convergence, but due to the expensive cost of the first one of them, it can only be considered for academical purposes an future development of methods derived from S-FETI. The second one, needs further tests to asset his value.

Finally, the Block-2LM method was developed as an extension of FETI-2LM, using the same ideas that lead to the S-FETI method. The construction of a bigger search space containing the original Krylov space assures convergence. Sparsity of the search directions helps to speed up this method, reducing drastically the number of forward-backward reductions done in every iteration.

The results confirm a reduction in the number of iterations, although not as important as the one seen in S-FETI and also with no reduction in the total time spend, however this method is considered as a first step, since an sparse storage version is also needed in this case, knowing that the reduction achieved is not small enough and that we will want to solve some big cases presented in problems arising from electromagnetism, as the ones seen in [START_REF] Barka | Scalability of FETI-2LM methods on HPC clusters for antenna RCS applications[END_REF], [START_REF] Roux | Block Krylov Recycling Algorithms for FETI-2LM Applied to 3-D Electromagnetic Wave Scattering and Radiation[END_REF].

Introduction 1 15 1

 115 Contributions of this thesis . . . . . . . . . . . . . . . . . . . . . . . . 2 Iterative methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Krylov methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Conjugate gradient . . . . . . . . . . . . . . . . . . . . . . . . . 10 ORTHODIR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 Parallelization of Krylov methods . . . . . . . . . . . . . . . . . Hybrid FETI method 17 1.1 Basic FETI method . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.1.1 Model problem and discretization . . . . . . . . . . . . 18 1.1.2 FETI one lagrange multiplier . . . . . . . . . . . . . . . 25 1.1.3 Local preconditioner . . . . . . . . . . . . . . . . . . . . 30 1.1.4 FETI Algorithms . . . . . . . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Three subdomain divisions. . . . . . . . . . . . . . . . . . . . . 1.10 Boundary conditions for both preconditioners. . . . . . . . . . 1.11 Two material bar. . . . . . . . . . . . . . . . . . . . . . . . . . . 1.12 Hybrid-FETI Iterations versus Elements number. . . . . . . . . 1.13 a) Regular interface marking. b) Extra covering marking. . . . 1.14 a) Contact problem. b) Initial gap. c) Contact pressure. . . . . 1.15 First example, initial configuration and subdomain division. . 1.16 First example, final configuration (solution). . . . . . . . . . . . 1.17 Second example, initial configuration. . . . . . . . . . . . . . . 2.1 Preconditioner construction . . . . . . . . . . . . . . . . . . . . 2.2 4 Subdomain problem . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Local interfaces with completely different subdomains (Young's modulus 0, 1 ≤ E ≤ 1000) . . . . . . . . . . . . . . . . . . . . . . 2.4 Three subdomain subdivision and computed corrections . . . . 2.5 Subdomain point of view (Left): In red the coarse modes describing

Algorithm 1 3 : v 1 = g 0 g 0 4 :

 134 Arnoldi iteration algorithm 1: Initialization 2: g 0 = Ax 0b loop Construction of the j + 1 vector of the base 5: w = Av j 6:

Algorithm 2 3 : v 1 = g 0 g 0 4 :

 234 Algorithm of Lanczos 1: Initialization 2: g 0 = Ax 0b loop Construction of the j + 1 vector of the base of Lanczos 5: w = Av j 6:

( 17 )

 17 Ax • x) -2(Ax • y) + (Ay • y) = (A(xy) • (xy)) > 0 (16) whenever x y. Now, if α ∈]0, 1[, then α(1α) < 0, hence 1 2 α(α -1) [(Ax • x) -2(Ax • y) + (Ay • y)] < 0 (Theorem 0.8. The functional J (x) = 1 2 (Ax • x) -(b • x)admits and absolute minimum x who also verifies Ax = b Proof. J is strictly convex and lower bounded , because J (x) → +∞ when x → +∞. It is obviously differentiable with a value of DJ (x) • y = (Ax • y) -(b • y) = ((Axb) • y)

Algorithm 3 8 :

 38 Conjugate gradient method 1: Initialization 2: g 0 = Ax 0b 3: w 0 = g 0 4: loop Iteration of the CG method 5:ρ p-1 = -(g p-1 • w p-1 )/(Aw p-1 • w p-1 ) 6:x p = x p-1 + ρ p-1 w p-1 7:g p = g p-1 + ρ p-1 Aw p-1 if (g p • g p )/(b • b) < 2 then 1 = (g p • Aw p-1 )/(Aw p-1 • w p-1 )12:

  the continuity relation 1.27 and use the fact that x

) 3 : 5 :

 35 g 0 = P D -1 P (Fλ 0d) 4: w 0 = g 0 Fw 0 = P D -1 P Fw 0 6: loop ORTHODIR Iteration from p = 0, 1, . . . until convergence 7: ρ p = -(Fw p ) T g p (Fw p ) T (Fw p ) 8: λ p+1 = λ p + ρ p w p 9: g p+1 = g p + ρ p Fw p 10: loop Construction of the p + 1 vector of the base F T F-orthonormal 11: w p+1 = g p+1 12:

  global block equation system comes from assembling both contributions giving

  are equivalent if the sum of the augmentation matrices A non singular. This is of course satisfied when the matrix A (s) b is the Schur complement of the rest of the domain and easily satisfied when A (s) b is any consistent approximation of the Dirichlet to Neumann operator of the rest of the domain.

Finally

Figure 1 . 4 -

 14 Figure 1.4 -Γ j division example and crosspoint detail.

Figure 1 . 5 -

 15 Figure 1.5 -One way splitting.

Theorem 1 . 1 .

 11 In a case of a two-domain splitting, the simple (Jacobi) iterative algorithm for 2-Lagrange multiplier with augmented term equal to the complete outer Schur complement defined as in Equation1.78 converges in one iteration at most.Proof. From the definition of A

  the system is the exact condensation of local augmented problem

  with these coefficients.

Figure 1 . 6 -

 16 Figure 1.6 -Nodes numbering in subdomain Ω (s) .

Figure 1 . 7 -

 17 Figure 1.7 -Interface patch of size p = 1 (red) with one and two layers d = 1, 2 (blue).

FormallyFigure 1 . 10 -

 110 Figure 1.10 -Boundary conditions for both preconditioners.

Figure 1 .

 1 Figure 1.12 -Hybrid-FETI Iterations versus Elements number.

Figure 1 .

 1 Figure 1.13 -a) Regular interface marking. b) Extra covering marking.

Figure 1 .

 1 Figure 1.14 -a) Contact problem. b) Initial gap. c) Contact pressure.

Figure 1 . 15 -

 115 Figure 1.15 -First example, initial configuration and subdomain division.

Figure 1 .

 1 Figure 1.16 -First example, final configuration (solution).

Figure 1 . 17 -

 117 Figure 1.17 -Second example, initial configuration.

5 )

 5 Multiplying by the local matrix, the interface response corrections on the inter-face Γ are computed by δf

b 21 )Figure 2 . 1 -

 2121 Figure 2.1 -Preconditioner construction

. 1 From

 1 Figure 2.1 (a) and (b) shows how the increments δx (s) b are chosen so that the corrected interface approximate solution x(s) b satisfy the interface compatibility. Following Equation 2.21 we compute the nodal interface responses required to maintain the increment δx (s) b . This mapping δf (s) b do not satisfy the interface equilibrium unless the corrected solutions x(s) b are the exact final solution. In Figure 2.1 (d) we compute the Lagrange multiplier corrections z.Interface Lagrange multipliers are naturally self-equilibrated in the sense that they result in interface values B (s) T z that are in equilibrium. Hence, it is in general impossible to define z such that B (s) T z restores the interface response corrections δf (s)

b

  acting on an interface d.o.f. are zero, i.e ∀i ∈ Γ ,

  And in general, this coefficients are computed for every edge Γ i by β (s)

  a compatible solution are, in general, different in two connected d.o.f.

B( 2 )

 2 T B(1) δf(1),p

) 3 :

 3 g 0 = P T (Fλ 0d) 4: Z 0 = [. . . , D (s) -1 g 0 , . . . ], s = 1, N s 5: W 0 = P Z 0 6: loop Iterate p = 0, 1, 2, ... λ p + W p ρ p 11: g p+1 = g p + P T FW p ρ p 12: Z p+1 = [. . . , D (s) -1 g p+1 , . . . ], s = 1,

) 2 Figure 2 . 3 -

 223 Figure 2.3 -Local interfaces with completely different subdomains (Young's modulus 0, 1 ≤ E ≤ 1000)

Figure 2 . 4 -

 24 Figure 2.4 -Three subdomain subdivision and computed corrections

= λ p + W p ρ p 11 :

 11 g p+1 = g p + P T FW p ρ p 12: Z p+1 = [. . . , D (s) -1

( 1 ( 2

 12 AG) T (FAG) : Computed once for all and sparse stored Z T j FZ p : Computed at every iteration for j = 1, p -1 and sparse stored Z T p (FAG) : Computed at every iteration and sparse stored (2.158) D p : Computed at every iteration and overwritten by ∆ p ∆ p : Small matrix stored every iteration (AG) T (FAG)∆ p : Small matrix stored every iteration Φ jp : Small matrices stored every iteration, for j = 1, p -1 B jp : Small matrices stored every iteration, for j = 1, p -.159) With all this considerations, we can describe the S-FETI method with this new search directions storage in Algorithm 10 Algorithm 10 S-FETI with Sparse Storage 1: Initialization 2: λ 0 = AG[G T AG] -1 (-R T c) 3: g 0 = Fλ 0d 4: loop Iterate p = 0, 1, 2, ... until convergence 5:

Figure 2 . 6 -

 26 Figure 2.6 -One-way split of six subdomains

Figure 2 . 7 -

 27 Figure 2.7 -Cube with 125 subdomains.

Figure 2 . 8 -

 28 Figure 2.8 -Time of Pardiso vs Dissection for different element number in each subdomain.

Figure 2 . 9 -

 29 Figure 2.9 -Checkerboard cube with 125 subdomains.

Figure 2 . 10 -

 210 Figure 2.10 -Max local memory usage in cube problem for 75 (left) and 300 (right) thousand elements per subdomain.

Algorithm 12 Flexible

 12 Conjugate Gradient Method 1: Initialization 2: Choose an approximation vector λ 0 3: g 0 = Fλ 0d 4: w 0 = B[g 0 ] 5: loop p = 0, 1, 2, ... until convergence 6: ρ p = -(w p , g p ) (w p , Fw p ) 7: λ p+1 = λ p + ρ p w p 8: g p+1 = g p + ρ p Fw p 9: w p+1 = B[g p ] 10:

2 :

 2 ρ p = -W T p g p 3: λ p+1 = λ p + W p ρ p 4:

Algorithm 14 Preconditioning procedure B 1 : 1 Γ

 11 Input Residual vector g 2: Z = [. . . , D (s) -(sq) g, . . . ] Local Dirichlet preconditioner 3: W = P Z FETI Projection 4: N LL T N T = W T FW Cholesky decomposition 5: W = W N L -T 6: ρ = -W T g 7: w = W ρ 8: Return B[g] ← w the same as in chapter 2 so D (s) -1Γ (sq) is the local preconditioner partitioned into neighbour interfaces and the matrix N is the permutation matrix coming from the Cholesky factorization applied to W T FW .

Algorithm 15 FETI

 15 block to one direction Algorithm 1: Initialization2: λ 0 = AG[G T AG] -1 (-R T c) 3: g 0 = (Fλ 0d) 4: w 0 = B[P T g 0 ] 5: loop Iterate p = 0,1, 2, ... until convergence 6: α p = -(g p , w p ) (w p , Fw p ) 7: λ p+1 = λ p + α p w p 8: g p+1 = g p + α p Fw p 9: w p+1 = B[P T g p+1 ] 10:

2 : 1 Γ

 21 W p = [. . . , D (s) -sq g p , . . . ] 3: EDE T = W T p FW p 4: W p = W p ED -1/2 T 5:

1 Γ

 1 1.2). Algorithm 17 FETI block with reduced saved directions algorithm 1: Initialization 2: λ 0 = AG[G T AG] -1 (-R T c) 3: g 0 = (Fλ 0d) 4: loop Iterate p = 0, 1, 2, ... until convergence 5: Z = . . . , D (s) -(sq) g p , . . .

Figure 3 . 1 -

 31 Figure 3.1 -Iterations versus percentage of computed directions. Example 1

Figure 4 . 1 -

 41 Figure 4.1 -Two subdomains with Robin Interface Condition

5 ) 7 )Ai

 57 The second condition is the interface equilibrium, which is the last line in the global block formulation K Now using the last line of the local problems, we can write the equilibrium condition as eliminating the inner unknowns on the Robin local problem, we can obtain the explicit relation between the trace of the solution on the interface x (s) b and the discrete augmented flux λ the condensed right-hand side. Replacing x (s) b by their values as functions of λ (s) b coming from previous explicit relation into the mixed equations 4.9, we can write the condensed interface problem of the form Fλ = d associated to the FETI-2LM method

Algorithm 18 FETI-2LM algorithm 1 :

 181 Initialization 2: λ 0 = 0 3: g 0 = (Fλ 0d) 4: w 0 = g 0 5: loop ORTHODIR Iteration from p = 0, 1, . . . until convergence 6: ρ p = -(Fw p ) T g p (Fw p ) T (Fw p ) 7: λ p+1 = λ p + ρ p w p 8: g p+1 = g p + ρ p Fw p 9: loop Construction of the p + 1 vector of the base F T F-orthonormal 10: w p+1 = g p+1 w p+1 = Fw p is replaced by the gradient 11:

s = 1 ,

 1 . . . , N s , i = 1, . . . , n (s)(4.18)

Algorithm 19 8 :Eliminates useless directions and 9 :F T F-orthogonalizes blocks 10 :

 198910 Block FETI-2LM method 1: Initialization 2: λ 0 = 0 3:g 0 = Fλ 0b 4: Z 0 = [. . . , g (s),i 0 , . . . ], s = 1, N s , i = 1, n (s) 5: W 0 = Z 0 6: loop Block ORTHODIR Iteration from p = 0, 1 . . . convergence 7: N p L p L T p N T p = (FW p ) T (FW p ) Cholesky factorization W p = W p N p L -T p FW p = FW p N p L -T p ρ p = -(FW p ) T g p 11:λ p+1 = λ p + W p ρ p 12:g p+1 = g p + FW p ρ p 13:loop Construction of the p + 1 vector of the base F T F-orthonormal 14: Z p+1 = [. . . , g (s),i p+1 , . . . ], s = 1, N s , i = 1, n (s)15:

  for the case of the Block-2LM, at the interface level, we have toSUBDOMAIN INTERFACEΩ (s) Ω (s)Ω(q)

Figure 4 . 2 -

 42 Figure 4.2 -Left:In red the coarse modes describing Z owned by subdomain Ω (s) . In dotted lines, the modes shared between Ω (s) and the subdomains involved in multiplying by the FETI-2LM operator. Right: In dotted lines the coarse modes describing FZ owned by local interface of Ω (s) and conversely they are the non null modes in the interface in red.

  . 33 1.2 FETI with two lagrange multipliers . . . . . . . . . . . . . . . . 35 1.2.1 FETI-2LM . . . . . . . . . . . . . . . . . . . . . . . . . . 35 1.2.2 Arbitrary mesh partition . . . . . . . . . . . . . . . . . . 38 1.2.3 Optimal Interface Boundary Conditions . . . . . . . . . 41 1.3 New FETI as an hybrid between one and two lagrange methods 48 1.3.1 Development . . . . . . . . . . . . . . . . . . . . . . . . 48 1.3.2 Extension to a general problem . . . . . . . . . . . . . . Preconditioner . . . . . . . . . . . . . . . . . . . . . . . . 58 1.3.4 Implementation . . . . . . . . . . . . . . . . . . . . . . . 61 1.4 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . 65 1.4.1 Two material bar . . . . . . . . . . . . . . . . . . . . . . 66 1.4.2 Contact Problem . . . . . . . . . . . . . . . . . . . . . . 71 1.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 Introduction and preliminarities . . . . . . . . . . . . . . . . . . 82 2.1.1 Dirichlet preconditioner for two subdomains . . . . . . 82 2.1.2 Consistent preconditioners . . . . . . . . . . . . . . . . . 84 2.1.3 Simultaneous FETI . . . . . . . . . . . . . . . . . . . . . 98 2.1.4 The Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 101 2.1.5 Cost and implementation of S-FETI . . . . . . . . . . . . 110 2.2 Sorting search directions in S-FETI . . . . . . . . . . . . . . . . 119 2.2.1 Linear dependence in block FETI directions . . . . . . . 120 2.2.2 Cholesky factorization with complete pivoting . . . . . 121 2.2.3 Diagonalization of search directions block . . . . . . . . 123 2.3 Memory usage in S-FETI . . . . . . . . . . . . . . . . . . . . . . 125 2.3.1 New sparse storage . . . . . . . . . . . . . . . . . . . . . 126 2.3.2 Reconstruction of search directions . . . . . . . . . . . . 126 2.3.3 Implementation details and exploitable parallelism . . 136 2.4 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . 142 2.4.1 FETI and S-FETI . . . . . . . . . . . . . . . . . . . . . . . 143 2.4.2 Decomposition of W t FW . . . . . . . . . . . . . . . . . . 146 2.4.3 S-FETI with sparse storage . . . . . . . . . . . . . . . . . 149 2.4.4 General comparison . . . . . . . . . . . . . . . . . . . . . 152 2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 Introduction and preliminaries . . . . . . . . . . . . . . . . . . 158 3.1.1 Method of Conjugate Directions . . . . . . . . . . . . . . 158 3.1.2 Flexible Conjugate Gradient . . . . . . . . . . . . . . . . 162 3.2 FETI with recursive preconditioner . . . . . . . . . . . . . . . . 164 3.2.1 One direction from S-FETI block . . . . . . . . . . . . . 164 3.2.2 Linear combination of directions from block . . . . . . . 168 3.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . 171 3.3.1 Storage of single direction . . . . . . . . . . . . . . . . . 171 3.3.2 Storage of reduced directions . . . . . . . . . . . . . . . 172 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Convergence of Hybrid-FETI preconditioner (Number of iterations) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Convergence of different marking Hybrid-FETI (Number of iterations) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Convergence of the different methods (Number of iterations) . 1.4 Convergence of the three contact examples for the different FETI methods (Number of iterations) . . . . . . . . . . . . . . . . . . Two subdomain splitting . . . . . . . . . . . . . . . . . . . . . . 1.2 Two subdomain divisions with duplicated nodes. . . . . . . . . 1.3 Multiple interface node. . . . . . . . . . . . . . . . . . . . . . .

	xx	Contents
	1.3.3 2 Block FETI methods 2.1 3 Block strategies as a preconditioner 3.1 List of Tables 1.1	81 157
		56
	xix

2.1 Time for forward-backward substitution on multi-core processor. 2.2 Direction stored in a 5x5x5 cube configuration . . . . . . . . . 2.3 Difference in subdomain division versus local interface division 2.4 Quadruple vs Double comparison . . . . . . . . . . . . . . . . . 2.5 Iterations/Search Directions results for different values of ε. . . 2.6 SPARSE vs SPARSE-OPT . . . . . . . . . . . . . . . . . . . . . . 2.7 FULL vs SPARSE-OPT . . . . . . . . . . . . . . . . . . . . . . . 2.8 Comparative between variations of S-FETI with corner interfaces. 125 subdomains and 150 thousand elements per subdomain. . 2.9 Comparative between variations of S-FETI without corner interfaces. 125 subdomains and 150 thousand elements per subdomain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Convergence comparative. . . . . . . . . . . . . . . . . . . . . . 4.1 64 subdomains . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 125 subdomains . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiii List of Figures 1.1 1.4 Γ j division example and crosspoint detail. . . . . . . . . . . . . 1.5 One way splitting. . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Nodes numbering in subdomain Ω (s) . . . . . . . . . . . . . . . . 1.7 Interface patch of size p = 1 (red) with one and two layers d = 1, 2 (blue). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 A contact problem.

  [START_REF] Farhat | The two-level FETI method Part II: Extension to shell problems, parallel implementation and performance results[END_REF] where x(s) is solution of the local discrete Neumann problem 1.29. Continuity relation 1.33 defines the condensed interface problem for FETI.

  indexes of the nodes inV Γ such that the minimum connectivity distance between each of these nodes and the node labelled j

	is lower or equal than p, p ∈ N}
	V p,d = { indexes of the nodes in V Ω (s) such that the minimum connectivity j
	j p distance between each of these nodes and the nodes in V
	is lower or equal than d, d ∈ N}
		(1.87)
	In other words, V	j p is a patch of radius p around some node j, and V p,d is the j
	neighbouring area of depth d of this patch.

11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

	7	8	9	10
	1	2	3	4	5	6
	= {1, 2, 8, 14, 20, 26, 25}		(1.88)
	V 13 1,2 = {1, 2, 8, 14, 20, 26, 25, 3, 9, 15, 21, 27, 33, 32, 31}
	These subsets correspond to the overlapping layers represented Figure 1.7.

Table 1 .

 1 1 -Convergence of Hybrid-FETI preconditioner (Number of iterations)namely the Dirichlet or Robin one. We will test this problem for three different number of subdomains. The number of global elements in each computation is constant of 75 thousand elements per subdomain. Also, following the definition of Algorithm 7, we calculate the augmentation matrices for the 2LM interfaces by setting the patch size as 1 and the depth size as 3.

	Number of Hybrid-FETI Hybrid-FETI
	subdomains	Dirichlet	Robin
	64	20	18
	125	20	20
	250	19	21

Table 1 .

 1

		Number of Regular Extra interface
		subdomains marking	marking
		64	18	79
		125	20	76
		250	21	64
	Table 1.2 -Convergence of different marking Hybrid-FETI (Number of iterations)
	Number of FETI-1LM Unsymmetric FETI-2LM Hybrid-FETI
	subdomains		FETI-1LM
	64	13	15	151	14
	125	14	14	283	15
	250	13	17	715	15

3 -Convergence of the different methods (Number of iterations)

Table 1 . 4

 14

	Subdomains	Global	FETI-1LM FETI-2LM Hybrid-FETI Hybrid-FETI
		equations			(Robin)	(Dirichlet)
	6	956	X	29	34	34
	25	31954	61	119	50	51
	41	5540	X	111	69	74

-Convergence of the three contact examples for the different FETI methods (Number of iterations)

needed. The different methods are then tested to compare the convergence rate. In the Table

1

.4 we can see the number of iterations in each method for different problem sizes.

Table 2 .

 2 3 -Difference in subdomain division versus local interface division

		Total	Iterations
	decomposition	Search Directions
	Local subdomain	1125	9
	Local interface	6356	4

Table 2 .

 2 4 -Quadruple vs Double comparison well as the total iterations.

	ε	Implementation Search directions	Minimal	Iterations
			per iteration	error	to error
	10 -4	QUAD	1247	6.0537x10 -5	8
	10 -4	DOUBLE	1245	9.0064x10 -5	8
	10 -6	QUAD	1356	1.7085x10 -3	6
	10 -6	DOUBLE	1356	1.7131x10 -3	6
	10 -10	QUAD	1498	6.4209x10 -3	5
	10 -10	DOUBLE	1494	1.6876x10 -3	5

Table 2 .

 2 -4 is more than enough to reduce 5 -Iterations/Search Directions results for different values of ε.

	ε	Implementation Search directions	Minimal	Iterations
			per iteration	error	to error
	10 -4	CHOLESKY	1245	9.0064x10 -5	8
	10 -4	EIGEN	1183	1.7590x10 -8	11
	10 -6	CHOLESKY	1356	1.7131x10 -3	6
	10 -6	EIGEN	1299	1.7082x10 -5	8
	10 -10	CHOLESKY	1494	1.6876x10 -3	5
	10 -10	EIGEN	1474	2.2169x10 -3	5

Table 2 .

 2 [START_REF] Axelsson | Variable-step multilevel preconditioning methods, I: Self-adjoint and positive definite elliptic problems[END_REF] -SPARSE vs SPARSE-OPT also change the full reorthogonalization process by making only one process to orthogonalize each saved block and then sending the information to the rest of the processes, see the last point in section 2.3.3.

		Elements	Time Max local memory
		Subdomain	(s)	(Gb)
	SPARSE-S-FETI	103823	571.1	4.91
	SPARSE-OPT-S-FETI	103823	289.4	4.63
	SPARSE-S-FETI	148877	648.1	6.57
	SPARSE-OPT-S-FETI	148877	339.8	6.13

Table 2 .

 2 8 -Comparative between variations of S-FETI with corner interfaces. 125 subdomains and 150 thousand elements per subdomain.

	Implementation Precon-	Decom-Max local Iterations Time
		ditioner position	memory		(s)
				(Gb)
	FULL	Dirichlet Cholesky	5.98	7	429.0
	FULL	Dirichlet	Eigen	5.93	7	432.4
	FULL	Lumped Cholesky	4.67	8	507.6
	FULL	Lumped	Eigen	5.02	11	825.0
	SPARSE-OPT	Dirichlet Cholesky	6.27	7	340.8
	SPARSE-OPT	Dirichlet	Eigen	6.32	7	331.4
	SPARSE-OPT	Lumped Cholesky	5.43	8	410.4
	SPARSE-OPT	Lumped	Eigen	7.77	11	812.5
	Implementation Precon-	Decom-Max local Iterations Time
		ditioner position	memory		(s)
				(Gb)
	FULL	Dirichlet Cholesky	5.21	12	119.6
	FULL	Dirichlet	Eigen	5.21	12	123.4
	FULL	Lumped Cholesky	4.46	24	366.9
	FULL	Lumped	Eigen	4.58	25	390.0
	SPARSE-OPT	Dirichlet Cholesky	4.43	12	59.7
	SPARSE-OPT	Dirichlet	Eigen	4.41	12	60.6
	SPARSE-OPT	Lumped Cholesky	4.31	24	158.3
	SPARSE-OPT	Lumped	Eigen	4.09	25	171.4

Table 2 .

 2 9 -Comparative between variations of S-FETI without corner interfaces. 125 subdomains and 150 thousand elements per subdomain.

  restriction is imposed in them. The norm considered here is the one induced by F, because, as already said, is a symmetric positive definitive operator.The resume of the algorithm is done in 11. We can see that the CD-method

	and no further Algorithm 11 Method of Conjugate Directions
	1: Initialization
	2: Choose an approximation vector λ 0
	3: g 0 = Fλ 0 -d
	4: Choose any non-zero vector w 0
	5: loop p = 0, 1, 2, ... until convergence
	6:	ρ p = -	(w p , g p ) (w p , Fw p )
	7:	λ p+1 = λ p + ρ p w p
	8:	g p+1 = g p + ρ p Fw p
	9:
				.2)
		are selected to be mutually conjugate, that is
			(w i , Fw j ) = 0,	∀i j	(3.3)

Table 3 .

 3 1 -Convergence comparative.

		Elements	Iterations
		Subdomains
	FETI-1LM	74088	17
	S-FETI	74088	7
	SINGLE_DIR-S-FETI	74088	14
	FETI-1LM	195112	18
	S-FETI	195112	7
	SINGLE_DIR-S-FETI	195112	15

We can use this same preconditioner as a first approximation for our method, in this case the definition will be similar to the one in FETI-1LM, but we apply it only on the interface edges Γ j marked as 1LM.

Acknowledgements

Thanks to my family, friends and people from the University. xvii Acknowledgements Contents

compatible. In this case, the consistent corrections are δx (1) 1 =

x(1) 1x

(1) 1

= β

(2) 1 x

(2)

1x

(1) 1

+ β

(3)

(1)

2x

(2) 2 (2.51) the interface responses are again computed using Equation 2.21, so the corrections for the Lagrange multipliers are build using the same weighting

1 δf

(1)

1 δf

(2)

1 δf

(3) 1

1 δf

(3)

1 δf

(1)

1 δf

(2) As for the convergence of the method we have the following theorem Theorem 3.2. The CD-method is an p-step method (p ≤ n) in the sense that at the p-iteration the estimate λ p is the solution λ

Proof. Let p be the first integer, such that

Clearly m ≤ n since the vectors are w 0 , w 1 , . . . are linearly independent. On the other hand, we may choose scalar α 0 , . . . , α p-1 such that

Hence, As usual the global error and solution jump across the interfaces will be our stopping criterion, when both of them are less to 10 -4 . We have diminished this value to be able to use the Block version in cases where the memory is an important constraint.

In the tables 4.1 and 4.2 we can see a general comparative of both methods. We have set the values for the parameter patch_size = 2 and patch_depth = 3, see Algorithm 7, also the number of elements per subdomain is change, as well as the total number of subdomains (cubes). The values of the parameter ν are 10 0 and 10 5 in a checkerboard configuration Figure 2.9.

In the results we can see the reduction in the number of iteration when comparing the two methods, however its values are not small enough to compensate for the time spend in each iteration. This slowdown in the total time is mainly due to the full reconjugation, where after a certain number of total iterations the number of matrix-matrix products needed is too costly.

Index

Hybridization of FETI Methods

Abstract

In this work new domain decomposition methods and new implementations for existing methods are developed. A new method based on previous domain decomposition methods is formulated. The classic FETI [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF] plus FETI-2LM [START_REF] Farhat | Two-level domain decomposition methods with Lagrange multipliers for the fast iterative solution of acoustic scattering problems[END_REF] methods are used to build the new Hybrid-FETI. The basic idea is to develop a new algorithm that can use both methods at the same time by choosing in each interface the most suited condition depending on the characteristics of the problem. By doing this we search to have a faster and more robust code that can work with configurations that the base methods will not handle it optimally by himself. The performance is tested on a contact problem. The following part involves the development of a new implementation for the S-FETI method [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF], the idea is to reduce the memory usage of this method, to make it able to work in larger problem. Different variation for this method are also proposed, all searching the reduction of directions stored each iteration of the iterative method. Finally, an extension of the FETI-2LM method to his block version as in S-FETI, is developed. Numerical results for the different algorithms are presented.

Keywords: numerical analysis, domain decomposition methods, algebra, scientific computation LJLL Laboratoire Jacques-Louis Lions Laboratoire Jacques-Louis Lions -4 place Jussieu -Université Pierre et Marie Curie -Boîte courrier 187 -75252 Paris Cedex 05 -France