Hussam Al Daas
email: aldaas.hussam@gmail.com

Laura Grigori
email: laura.grigori@inria.fr

Pierre Jolivet
email: pierre.jolivet@enseeiht.fr

Pierre-Henri Tournier
email: tournier@ljll.math.upmc.fr

A MULTILEVEL SCHWARZ PRECONDITIONER BASED ON A HIERARCHY OF ROBUST COARSE SPACES *

Keywords: domain decomposition, multilevel, elliptic problems, subspace correction AMS subject classifications. 65F08, 65F10, 65N55

In this paper we present a multilevel preconditioner based on overlapping Schwarz methods for symmetric positive definite (SPD) matrices. Robust two-level Schwarz preconditioners exist in the literature to guarantee fast convergence of Krylov methods. As long as the dimension of the coarse space is reasonable, that is, exact solvers can be used efficiently, two-level methods scale well on parallel architectures. However, the factorization of the coarse space matrix may become costly at scale. An alternative is then to use an iterative method on the second level, combined with an algebraic preconditioner, such as a one-level additive Schwarz preconditioner. Nevertheless, the condition number of the resulting preconditioned coarse space matrix may still be large. One of the difficulties of using more advanced methods, like algebraic multigrid or even two-level overlapping Schwarz methods, to solve the coarse problem is that the matrix does not arise from a partial differential equation (PDE) anymore. We introduce in this paper a robust multilevel additive Schwarz preconditioner where at each level the condition number is bounded, ensuring a fast convergence for each nested solver. Furthermore, our construction does not require any additional information than for building a two-level method, and may thus be seen as an algebraic extension.

1. Introduction. We consider the solution of a linear system of equations (1.1) Ax = b, where A ∈ R n×n is a symmetric positive definite (SPD) matrix, b ∈ R n is the righthand side, and x ∈ R n is the vector of unknowns. To enhance convergence, it is common to solve the preconditioned system

M -1 Ax = M -1 b.
Standard domain decomposition preconditioners such as block Jacobi, additive Schwarz, and restricted additive Schwarz methods are widely used [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF][START_REF] Chan | Domain decomposition algorithms[END_REF][START_REF] Cai | A restricted additive Schwarz preconditioner for general sparse linear systems[END_REF]. In a parallel framework, such preconditioners have the advantage of relatively low communication costs. However, their role in lowering the condition number of the system typically deteriorates when the number of subdomains increases. Multilevel approaches have shown a large impact on enhancing the convergence of Krylov methods [START_REF] Spillane | Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps[END_REF][START_REF] Dolean | An introduction to domain decomposition methods[END_REF][START_REF] Brezina | Algebraic Multigrid Based on Element Interpolation (AMGe)[END_REF][START_REF] Mandel | Multispace and multilevel BDDC[END_REF][START_REF] Jones | AMGe Based on Element Agglomeration[END_REF][START_REF] Chartier | Spectral AMGe (ρAMGe)[END_REF][START_REF] Kalchev | Upscaling of mixed finite element discretization problems by the spectral AMGe method[END_REF][START_REF] Adams | Ultrascalable Implicit Finite Element Analyses in Solid Mechanics with over a Half a Billion Degrees of Freedom[END_REF][START_REF] Heinlein | A Three-Level Extension of the GDSW Overlapping Schwarz Preconditioner in Three Dimensions[END_REF][START_REF] Kong | A highly scalable multilevel Schwarz method with boundary geometry preserving coarse spaces for 3D elasticity problems on domains with complex geometry[END_REF][START_REF] Toivanen | A multilevel feti-dp method and its performance for problems with billions of degrees of freedom[END_REF][START_REF] Pechstein | Finite and boundary element tearing and interconnecting solvers for multiscale problems[END_REF]. In multigrid and domain decomposition communities, multilevel methods have proven their capacity of scaling up to large numbers of processors and tackling ill-conditioned systems [START_REF] Zampini | PCBDDC: A class of robust dual-primal methods in PETSc[END_REF][START_REF] Badia | Multilevel balancing domain decomposition at extreme scales[END_REF][START_REF] Jolivet | Scalable domain decomposition preconditioners for heterogeneous elliptic problems[END_REF]. While some preconditioners are purely algebraic [START_REF] Brezina | Algebraic Multigrid Based on Element Interpolation (AMGe)[END_REF][START_REF] Jones | AMGe Based on Element Agglomeration[END_REF][START_REF] Chartier | Spectral AMGe (ρAMGe)[END_REF][START_REF] Marques | Tuning the coarse space construction in a spectral AMG solver[END_REF][START_REF] Notay | An aggregation-based algebraic multigrid method[END_REF][START_REF] Henson | BoomerAMG: A parallel algebraic multigrid solver and preconditioner[END_REF][START_REF] Adams | Ultrascalable Implicit Finite Element Analyses in Solid Mechanics with over a Half a Billion Degrees of Freedom[END_REF], several multilevel methods are based on hierarchical meshing in both multigrid and domain decomposition communities [START_REF] Toselli | Domain Decomposition Methods -Algorithms and Theory[END_REF][START_REF] Chan | Domain decomposition algorithms[END_REF][START_REF] Mandel | Multispace and multilevel BDDC[END_REF][START_REF] Heinlein | A Three-Level Extension of the GDSW Overlapping Schwarz Preconditioner in Three Dimensions[END_REF][START_REF] Kong | A highly scalable multilevel Schwarz method with boundary geometry preserving coarse spaces for 3D elasticity problems on domains with complex geometry[END_REF]. Mesh coarsening depends on the geometry of the problem. One has to be careful when choosing a hierarchical structure since it can have a significant impact on the iteration count [START_REF] Kong | A highly scalable multilevel Schwarz method with boundary geometry preserving coarse spaces for 3D elasticity problems on domains with complex geometry[END_REF][START_REF] Mandel | Multispace and multilevel BDDC[END_REF]. In [START_REF] Kong | A highly scalable multilevel Schwarz method with boundary geometry preserving coarse spaces for 3D elasticity problems on domains with complex geometry[END_REF], the authors propose a multilevel Schwarz domain decomposition solver for the elasticity problem. Based on a heuristic approach and following the maximum independent set method [START_REF] Adams | Parallel Multigrid Solver for 3D Unstructured Finite Element Problems[END_REF], they coarsen the fine mesh while preserving the boundary in order to obtain a two-level method. This strategy is repeated recursively to build several levels. However, they do not provide a bound on the condition number of the preconditioned matrix of the multilevel method. Multilevel domain decomposition methods are mostly based on non-overlapping approaches [START_REF] Toselli | Domain Decomposition Methods -Algorithms and Theory[END_REF][START_REF] Chan | Domain decomposition algorithms[END_REF][START_REF] Mandel | Multispace and multilevel BDDC[END_REF][START_REF] Kong | A highly scalable multilevel Schwarz method with boundary geometry preserving coarse spaces for 3D elasticity problems on domains with complex geometry[END_REF][START_REF] Zampini | PCBDDC: A class of robust dual-primal methods in PETSc[END_REF][START_REF] Badia | Multilevel balancing domain decomposition at extreme scales[END_REF][START_REF] Pechstein | Finite and boundary element tearing and interconnecting solvers for multiscale problems[END_REF][START_REF] Toivanen | A multilevel feti-dp method and its performance for problems with billions of degrees of freedom[END_REF]. Two-level overlapping domain decomposition methods are well studied and provide robust convergence estimates [START_REF] Spillane | Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps[END_REF][START_REF] Dolean | An introduction to domain decomposition methods[END_REF][START_REF] Bjørstad | Does SHEM for Additive Schwarz Work Better than Predicted by Its Condition Number Estimate?[END_REF]. However, extending such a construction to more than two levels while preserving robustness is not straightforward. In [START_REF] Borzì | Parallel algebraic multilevel Schwarz preconditioners for a class of elliptic PDE systems[END_REF], the authors propose an algebraic multilevel additive Schwarz method. Their approach is inspired by algebraic multigrid strategies. One drawback of it is that it is sensitive to the number of subdomains. In [START_REF] Heinlein | A Three-Level Extension of the GDSW Overlapping Schwarz Preconditioner in Three Dimensions[END_REF], the authors suggest applying the two-level Generalized Dryja-Smith-Widlund preconditioner recursively to build a multilevel method. In this case, the condition number bound of the two-level approach depends on the width of the overlap, the diameter of discretization elements, and the diameter of the subdomains. They focus on the preconditioner for the three-level case. One drawback of their approach is that the three-level preconditioner requires more iterations than the two-level variant. In this paper, the only information from the PDE needed for the construction of the preconditioner consists of the local Neumann matrices at the fine level. These matrices correspond to the integration of the bilinear form in the weak formulation of the studied PDE on the subdomain-decomposed input mesh. No further information is necessary: except on the fine level, our method is algebraic and does not depend on any coarsened mesh or auxiliary discretized operator. For problems not arising from PDE discretization, one needs to supply the local SPSD matrices on the finest level. In [START_REF] Daas | A class of efficient locally constructed preconditioners based on coarse spaces[END_REF], a subset of the authors propose a fully algebraic approximation for such matrices. However, their approximation strategy is heuristic and may not be effective in some cases.

Our preconditioner is based on a hierarchy of coarse spaces and is defined as following. At the first level, the set of unknowns is partitioned into N 1 subdomains and each subdomain has an associated matrix A 1,j = R 1,j AR 1,j obtained by using appropriate restriction and prolongation operators R 1,j and R 1,j respectively, defined in the following section. The preconditioner is formed as an additive Schwarz preconditioner coupled with an additive coarse space correction, defined as,

M -1 = M -1 1 = V 1 A -1 2 V 1 + N1 j=1 R 1,j A -1 1,j R 1,j ,
where V 1 is a tall-and-skinny matrix spanning a coarse space obtained by solving for each subdomain j = 1 to N 1 a generalized eigenvalue problem involving the matrix A 1,j and the Neumann matrix associated with subdomain j. The coarse space matrix

is A 2 = V 1 AV 1 .
This is equivalent to the GenEO preconditioner, and is described in detail in [START_REF] Spillane | Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps[END_REF] and recalled briefly in section 2. The dimension of the coarse space is proportional to the number of subdomains N 1 . When it increases, factorizing A 2 by using a direct method becomes prohibitive, and hence the application of A -1 2 to a vector should also be performed through an iterative method.

Our multilevel approach defines a hierarchy of coarse spaces V i and coarse space matrices A i for i = 2 to any depth L + 1, and defines a preconditioner

M -1 i such that the condition number of M -1 i A i is bounded. The depth L + 1 is chosen such that the
This manuscript is for review purposes only.

coarse space matrix A L+1 can be factorized efficiently by using a direct method. At each level i, the graph of the coarse space matrix A i is partitioned into N i subdomains, and each subdomain j is associated with a local matrix A i,j = R i,j A i R i,j obtained by using appropriate restriction and prolongation operators R i,j and R i,j , respectively.

The preconditioner at level i is defined as,

M -1 i = V i A -1 i+1 V i + Ni j=1 R i,j A -1 i,j R i,j ,
where the coarse space matrix is

A i+1 = V i A i V i .
One of the main contributions of the paper concerns the construction of the hierarchy of coarse spaces V i for levels i going from 2 to L, that are built algebraically from the coarse space of the previous level V i-1 . This construction is based on the definition of local symmetric positive semi-definite (SPSD) matrices associated with each subdomain j at each level i that we introduce in this paper. These matrices are obtained by using the local SPSD matrices of the previous level i -1 and the previous coarse space V i-1 . They are then involved, with the local matrices A i,j , in concurrent generalized eigenvalue problems solved for each subdomain j that allows to compute the local eigenvectors contributing to the coarse space V i .

We show in Theorem 5.3, section 5, that the condition number of M -1 i A i is bounded and depends on the maximum number of subdomains at the first level that share an unknown, the number of distinct colors required to color the graph of A i so that span{R i,j } 1 j Ni of the same color are mutually A i -orthogonal, and a user defined tolerance τ . It is thus independent of the number of subdomains N i .

The main contribution of this paper is based on the combination of two previous works on two-level additive Schwarz methods [START_REF] Daas | A class of efficient locally constructed preconditioners based on coarse spaces[END_REF][START_REF] Spillane | Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps[END_REF]. The coarse space proposed in [START_REF] Spillane | Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps[END_REF] guarantees an upper bound on the condition number that can be prescribed by the user. The SPSD splitting in the context of domain decomposition presented in [START_REF] Daas | A class of efficient locally constructed preconditioners based on coarse spaces[END_REF] provides an algebraic view for the construction of coarse spaces. The combination of these two works leads to a robust multilevel additive Schwarz method. Here, robustness refers to the fact that at each level, an upper bound on the condition number of the associated matrix can be prescribed by the user a priori. The rest of the paper is organized as follows. In the next section, we present the notations used throughout the paper. In section 2, we present a brief review of the theory of one-and two-level additive Schwarz methods. We extend in section 3 the class of SPSD splitting matrices presented in [START_REF] Daas | A class of efficient locally constructed preconditioners based on coarse spaces[END_REF] in order to make it suitable for multilevel methods. Afterwards, we define the coarse space at level i based on the extended class of local SPSD splitting matrices associated with this level. Section 4 describes the partitioning of the domain at level i + 1 from the partitioning at level i. In Section 5, we explain the computation of the local SPSD matrices associated with each subdomain at level i + 1. We compute them using those associated with subdomains at level i. Section 6 presents numerical experiments on highly challenging diffusion and linear elasticity problems in two-and three-dimensional problems. We illustrate the theoretical robustness and practical usage of our proposed method by performing strong scalability tests up to 8,192 processes.

Context and notation. By convention, the finest level, on which (1.1) is defined, is the first level. A subscript index is used in order to specify which level an entity is defined on. In the case where additional subscripts are used, the first subscript always denotes the level. For the sake of clarity, we omit the subscript cor-This manuscript is for review purposes only.

responding to level 1 when it is clear from context, e.g., matrix A. Furthermore, the subscripts i and j always refer to a specific level i and its subdomain j, respectively.

The number of levels is L + 1. Let A i ∈ R ni×ni denote symmetric positive definite matrices, each corresponding to level i = 1, . . . , L + 1. We suppose that a direct solver can be used at level L + 1 to compute an exact factorization of A L+1 .

Let B ∈ R p×q be a matrix. Let P ⊂ 1; p and Q ⊂ 1; q be two sets of indices. The concatenation of P and Q is represented by [P, Q]. We note that the order of the concatenation is important. B(P, :) is the submatrix of B formed by the rows whose indices belong to P . B(:, Q) is the submatrix of B formed by the columns whose indices belong to Q. B(P, Q) = (B(P, :)) (:, Q). The identity matrix of size p is denoted I p . We suppose that the graph of A i is partitioned into N i nonoverlapping subdomains, where N i n i and N i+1 N i for i = 1, . . . , L. We note that partitioning at level 1 can be performed by using a graph partitioning library such as ParMETIS [START_REF] Karypis | Multilevel k-way partitioning scheme for irregular graphs[END_REF] or PT-SCOTCH [START_REF] Chevalier | PT-SCOTCH: A tool for efficient parallel graph ordering[END_REF]. Partitioning at greater levels will be described later in section 4. In the following, we define for each level i = 1, . . . , L notations for subsets and restriction operators that are associated with the partitioning. Let Ω i = 1; n i be the set of unknowns at level i and let Ω i,j,I for j = 1, . . . , N i be the subset of Ω i that represents the unknowns in subdomain j. We refer to Ω i,j,I as the interior unknowns of subdomain j. Let Γ i,j for j = 1, . . . , N i be the subset of Ω i that represents the neighbor unknowns of subdomain j, i.e., the unknowns at distance 1 from subdomain j through the graph of A i . We refer to Γ i,j as the overlapping unknowns of subdomain j. We denote Ω i,j = [Ω i,j,I , Γ i,j ], for j = 1, . . . , N i , the concatenation of interior and overlapping unknowns of subdomain j. We denote ∆ i,j , for j = 1, . . . , N i , the complementary of Ω i,j in Ω i , i.e., ∆ i,j = Ω i \ Ω i,j . In Figure 1.1, a triangular mesh is used to discretize a square domain. The set of nodes of the mesh is partitioned into 16 disjoint subsets Ω 1,j,I , which represent a non-overlapping decomposition, for j = 1, . . . , 16 (left). On the left, a matrix A 1 whose connectivity graph corresponds to the mesh is illustrated. The submatrix A 1 (Ω 1,j,I , Ω 1,j,I ) is associated with the non-overlapping subdomain j. Each submatrix A 1 (Ω 1,j,I , Ω 1,j,I ) is colored with a distinct color. The same color is used to color the region that contains the nodes in the non-overlapping subdomain Ω 1,j,I . Note that if two subdomains j 1 , j 2 are neighbors, the submatrix A 1 (Ω 1,j1,I , Ω 1,j2,I ) has nonzero elements. For j = 1, . . . , N i , we denote by n i,j,I , γ i,j and n i,j the cardinality of Ω i,j,I , Γ i,j and Ω i,j respectively.

Let R i,j,I ∈ R n i,j,I ×ni be defined as R i,j,I = I ni (Ω i,j,I , :).

Let R i,j,Γ ∈ R γi,j ×ni be defined as R i,j,Γ = I ni (Γ i,j , :).

Let R i,j ∈ R ni,j ×ni be defined as R i,j = I ni (Ω i,j , :).

Let R i,j,∆ ∈ R (ni-ni,j )×ni be defined as R i,j,∆ = I ni (∆ i,j , :).

Let P i,j = I ni ([Ω i,j,I , Γ i,j , ∆ i,j ], :) ∈ R ni×ni , be a permutation matrix associated with the subdomain j, for j = 1, . . . , N i . The matrix of the overlapping subdomain j, R i,j A i R i,j , is denoted A i,j . We denote D i,j ∈ R ni,j ,×ni,j , j = 1, . . . , N i , any set of non-negative diagonal matrices such that

I ni = Ni j=1 R i,j D i,j R i,j .
We refer to {D i,j } 1 j Ni as the algebraic partition of unity. Let V i ∈ R ni×ni+1 be a tall-and-skinny matrix of full rank. We denote S i the subspace spanned by the columns of V i . This subspace will stand for the coarse space associated with level i.

By convention, we refer to S i as subdomain 0 at level i. Thus, we have n i,0 = n i+1 .

This manuscript is for review purposes only. Left: a triangular mesh is used to discretize the unit square. The set of nodes of the mesh is partitioned into 16 disjoint subsets, non-overlapping subdomains, Ω 1,j,I for j = 1, . . . , 16. Right: Illustration of the matrix A 1 whose connectivity graph corresponds to the mesh on the left. The diagonal block j of A 1 corresponds to the non-overlapping subdomain Ω 1,j,I . Each submatrix A 1 (Ω 1,j,I , Ω 1,j,I ) is colored with a distinct color. The same color is used to color the region of the square that contains nodes in Ω 1,j,I .

The interpolation operator at level i is defined as:

R i,2 : Ni j=0 R ni,j → R ni (u j ) 0 j Ni → Ni j=0 R i,j u j .
(1.2) Finally, we denote V i,j the set of neighboring subdomains of each subdomain j at

level i for (i, j) ∈ 1; L × 1; N i . V i,j = {k ∈ 1; N i : Ω i,j ∩ Ω i,k = ∅}.
As previously mentioned, partitioning at level 1 can be performed by graph partitioning libraries such as ParMETIS [START_REF] Karypis | Multilevel k-way partitioning scheme for irregular graphs[END_REF] or PT-SCOTCH [START_REF] Chevalier | PT-SCOTCH: A tool for efficient parallel graph ordering[END_REF]. Partitioning at further levels will be defined later: the sets Ω i,j,I , Ω i,j,Γ , Ω i,j , and ∆ i,j for i > 1 are defined in subsection 4.2. The coarse spaces S i as well as the projection and prolongation operators V i and V i are defined in subsection 3.2. We suppose that the connectivity graph between the subdomains on each level is sparse. This assumption is not true in general, however, it is valid in structures based on locally constructed coarse spaces in domain decomposition as we show in this paper, see [18, Section 4.1 p.81] for the case of two levels.

Background.

In this section, we review briefly several theoretical results related to additive Schwarz preconditioners. We introduce them for the sake of completeness.

Lemma 2.1 (fictitious subspace lemma). Let

A ∈ R n A ×n A , B ∈ R n B ×n B be two
This manuscript is for review purposes only. symmetric positive definite matrices. Let R be an operator defined as

R : R n B → R n A v → Rv,
and let R be its transpose. Suppose that the following conditions hold:

1. The operator R is surjective.

2. There exists c u > 0 such that

(Rv) A (Rv) c u v Bv, ∀v ∈ R n B . 3. There exists c l > 0 such that for all v n A ∈ R n A , ∃v n B ∈ R n B |v n A = Rv n B and c l v n B Bv n B (Rv n B ) A (Rv n B ) = v n A Av n A .
Then, the spectrum of the operator

RB -1 R A is contained in the segment [c l , c u ].
Proof. We refer the reader to [START_REF] Dolean | An introduction to domain decomposition methods[END_REF]Lemma 7.4 p.164] or [START_REF]Decomposition and fictitious domains methods for elliptic boundary value problems[END_REF][START_REF] Nepomnyaschikh | Mesh theorems of traces, normalizations of function traces and their inversions[END_REF][START_REF] Griebel | On the abstract theory of additive and multiplicative Schwarz algorithms[END_REF] for a detailed proof.

Lemma 2.2. The operator R i,2 as defined in (1.2) is surjective.

Proof. The proof follows from the definition of R i,2 (1.2).

Lemma 2.3. Let k i,c for i = 1, . . . , L be the minimum number of distinct colors so that span{R i,j } 1 j Ni of the same color are mutually A i -orthogonal. Then, we have

(R i,2 u Bi ) A i (R i,2 u Bi ) (k i,c + 1) Ni j=0 u j R i,j A i R i,j u j , ∀u Bi = (u j ) 0 j Ni ∈ N j=0 R ni,j .
Proof. We refer the reader to [9, Theorem 12 p.93] for a detailed proof.

We note that at level i, the number k i,c is smaller than the maximum number of neighbors over the set of subdomains 1;

N i k i,c max 1 j Ni #V i,j .
Due to the sparse structure of the connectivity graph between the subdomains at level i, the maximum number of neighbors over the set of subdomains 1;

N i is independent of the number of subdomains N i . Then, so is k i,c . Lemma 2.4. Let u Ai ∈ R n A i and u Bi = {u j } 0 j Ni ∈ Ni j=0 R ni,j such that u Ai = R i,2 u Bi .
The additive Schwarz operator without any other restriction on the coarse space S i verifies the following inequality

Ni j=0 u j R i,j A i R i,j u j 2u Ai A i u Ai + (2k i,c + 1) Ni j=1 u j R i,j A i R i,j u j ,
where k i,c is defined in Lemma 2.3.

This manuscript is for review purposes only.

Proof. We refer the reader to [12, Lemma 7.12, p. 175] to view the proof in detail.

Lemma 2.5. Let A, B ∈ R m×m be two symmetric positive semi-definite matrices.

Let ker(A), range(A) denote the null space and the range of A respectively. Let P 0 be an orthogonal projection on range(A). Let τ be a positive real number. Consider the generalized eigenvalue problem,

P 0 BP 0 u k = λ k Au k , (u k , λ k ) ∈ range(A) × R.
Let P τ be an orthogonal projection on the subspace

Z = ker(A) ⊕ span {u k |λ k > τ } ,
then, the following inequality holds:

(2.1) (u -P τ u) B (u -P τ u) τ u Au, ∀u ∈ R m .
Proof. We refer the reader to [3, Lemma 2.4] and [START_REF] Dolean | An introduction to domain decomposition methods[END_REF]Lemma 7.7] for a detailed proof.

2.1. GenEO coarse space. In [START_REF] Spillane | Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps[END_REF][START_REF] Dolean | An introduction to domain decomposition methods[END_REF] the authors present the GenEO coarse space which relies on defining appropriate symmetric positive semi-definite (SPSD) matrices Ãj ∈ R n×n for j = 1, . . . , N . These are the unassembled Neumann matrices, corresponding to the integration on each subdomain of the operator defined in the variational form of the PDE. These matrices are local, i.e., R j,∆ Ãj = 0. Furthermore, they verify the relations

u Ãj u u Au, ∀u ∈ R n , u N j=1 Ãj u k GenEO u Au, ∀u ∈ R n ,
where k GenEO N is the maximum number of subdomains that share an unknown.

2.2.

Local SPSD splitting of an SPD matrix. In [START_REF] Daas | A class of efficient locally constructed preconditioners based on coarse spaces[END_REF], the authors present the local SPSD splitting of an SPD matrix. Given the permutation matrix P j , a local SPSD splitting matrix Ãj of A associated with subdomain j is defined as (2.2)

P j Ãj P j =   R j,I AR j,I R j,I AR j,Γ R j,Γ AR j,I Ãj Γ 0   ,
where Ãj Γ ∈ R γj ×γj satisfies the two following conditions: For all u ∈ R γj ,

• u R j,Γ AR j,I R j,I AR j,I -1 R j,I AR j,Γ u u Ãj γ u • u Ãj Γ u u R j,Γ AR j,Γ -R j,Γ AR j,∆ R j,∆ AR j,∆ -1 R j,∆ AR j,Γ u.
The authors prove that the matrices Ãj defined in such a way verify the following relations:

R j,∆ Ãj = 0, (2.3) u Ãj u ≤ u Au, ∀u ∈ R n , (2.4) u N j=1 Ãj u ku Au, ∀u ∈ R n , (2.5)
This manuscript is for review purposes only.

where k is a number that depends on the local SPSD splitting matrices and can be at most equal to the number of subdomains k N . The authors also show that the local matrices defined in GenEO [START_REF] Spillane | Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps[END_REF][START_REF] Dolean | An introduction to domain decomposition methods[END_REF] can be seen as a local SPSD splitting.

In [START_REF] Daas | A class of efficient locally constructed preconditioners based on coarse spaces[END_REF], the authors highlight that the key idea to construct a coarse space relies on the ability to identify the so-called local SPSD splitting matrices. They present a class of algebraically constructed coarse spaces based on the local SPSD splitting matrices. Moreover, this class can be extended to a larger variety of local SPSD matrices. This extension has the advantage of allowing to construct efficient coarse spaces for a multilevel structure in a practical way. This is discussed in the following section.

Extension of the class of coarse spaces.

In this section we extend the class of coarse spaces presented in [START_REF] Daas | A class of efficient locally constructed preconditioners based on coarse spaces[END_REF]. To do so, we present a class of matrices, that is larger than the class of local SPSD splitting matrices. This will be our main building block in the construction of efficient coarse spaces. Furthermore, this extension can lead to a straightforward construction of hierarchical coarse spaces in a multilevel Schwarz preconditioner setting.

Extension of the class of local SPSD splitting matrices.

Regarding the two-level additive Schwarz method, the authors of [START_REF] Daas | A class of efficient locally constructed preconditioners based on coarse spaces[END_REF] introduced the local SPSD splitting related to a subdomain as defined in (2.2). As it can be seen from the theory presented in that paper, it is not necessary to have the exact matrices R j,I AR j,I , R j,I AR j,Γ , and R j,Γ AR j,I in the definition of the local SPSD splitting in order to build an efficient coarse space. Indeed, the one and only necessary condition is to define for each subdomain j an SPSD matrix Ãj for j = 1, . . . , N such that:

R j,∆ Ãj = 0, u N j=1 Ãj u ku Au, ∀u ∈ R n , (3.1)
where k is a number that depends on the local SPSD matrices Ãj for j = 1, . . . , N .

The first condition means that Ãj has the local SPSD structure associated with subdomain j, i.e., it has the following form:

P j Ãj P j = Ãj I,Γ 0 0 0 ,
where Ãj I,Γ ∈ R nj ×nj . The second condition is associated with the stable decomposition property [START_REF] Xu | Theory of Multilevel Methods[END_REF][START_REF] Dolean | An introduction to domain decomposition methods[END_REF]. Note that with regard to the local SPSD matrices, the authors in [START_REF] Spillane | Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps[END_REF] only use these two conditions. That is to say, with matrices that verify conditions (3.1) the construction of the coarse space is straightforward through the theory presented in either [START_REF] Spillane | Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps[END_REF] or [START_REF] Daas | A class of efficient locally constructed preconditioners based on coarse spaces[END_REF]. To this end, we define in the following the local SPSD (LSPSD) matrix associated with subdomain j as well as the associated local filtering subspace that contributes to the coarse space.

Definition 3.1 (local SPSD matrices). An SPSD matrix Ãi,j ∈ R ni×ni is called local SPSD (LSPSD) with respect to subdomain j if

• R i,j,∆ Ãi,j = 0, • u Ni j=1 Ãi,j u k i u A i u,
where k i > 0.

This manuscript is for review purposes only.

We note that the local SPSD splitting matrices form a subset of the local SPSD matrices.

3.2. Multilevel coarse spaces. This section summarizes the steps to be performed in order to construct the coarse space at level i once we have the LSPSD matrices associated with each subdomain at that level.

Definition 3.2 (coarse space based on LSPSD matrices). Let Ãi,j ∈ R ni×ni for j = 1, . . . , N i be LSPSD matrices. Let D i,j ∈ R ni,j for j = 1, . . . , N i be the partition of unity. Let τ i > 0 be a given number. For a subdomain j ∈ 1; N i , let

G i,j = D i,j R i,j A i R i,j D i,j .
Let Pi,j be the projection on range(R i,j Ãj R i,j ) parallel to ker(R i,j Ãj R i,j ). Let K i,j = ker(R i,j Ãi,j R i,j ). Consider the generalized eigenvalue problem:

Pi,j G i,j Pi,j u i,j,k = λ i,j,k R i,j Ãi,j R i,j u i,j,k , (u i,j,k , λ i,j,k ) ∈ range(R i,j Ãi,j R i,j ) × R. (3.2) Set (3.3) Z i,j = K i,j ⊕ span {u i,j,k |λ i,j,k > τ i } .
Then, the coarse space associated with LSPSD matrices Ãi,j for j = 1, . . . , N i at level i is defined as:

(3.4) S i = Ni j=1 R i,j D i,j Z i,j .
Following notations from section 1, the columns of V i span the coarse space S i . The matrix A i+1 is defined as:

(3.5) A i+1 = V i A i V i .
The local SPSD splitting matrices at level 1 will play an important role in the construction of the LSPSD matrices at subsequent levels. In the following, we present an efficient approach for computing LSPSD matrices for levels greater than 1.

4. Partitioning for levels strictly greater than 1. In this section, we explain how to obtain the partitioning sets Ω i,j,I for (i, j) ∈ 2; L × 1; N i . Once the sets Ω i,j,I for j = 1, . . . , N i are defined at level i, the following elements are readily available: sets Γ i,j , ∆ i,j , and Ω i,j ; restriction operators R i,j,I , R i,j,Γ , R i,j,∆ , and R i,j ; permutation matrices P i,j for j = 1, . . . , N i . The partition of unity is constructed in an algebraic way. The mth diagonal element of D i,j is 1 if m n i,j,I and 0 otherwise. This manuscript is for review purposes only.

D i,1 Z i,1 D i,3 Z i,3                                           (R i,1 D i,1 Z i,1 ) A i (R i,1 D i,1 Z i,1 ) (R i,3 D i,3 Z i,3 ) A i (R i,4 D i,4 Z i,4 )                   Fig. 4.2.
Illustration of the correspondence of indices between the columns of V i (left) and the rows and columns of A i+1 (right). Having no overlap in V i is possible through a non-overlapping partition of unity.

4.2. Heritage from superdomains. Let e i,j be the set of indices of the vectors that span R i,j D i,j Z i,j in the matrix V i for some (i, j) ∈ 1; L -1 × 1; N i , see Figure 4.2. We define Ω i+1,j,I = ∪ k∈Gi,j e i,k , for j = 1, . . . , N i+1 . We denote Ω i+1,j,Γ the subset of 1; n i+1 \Ω i+1,j,I whose elements are at distance 1 from Ω i+1,j,I through the graph of A i+1 . We note that Ω i+1,j,Γ ⊂ p∈Gi,j k∈Vi,p e i,k , where V i,j represents the set of subdomains that are neighbors of subdomain j at level i for j = 1, . . . , N i . The overlapping subdomain j is defined by the set Ω i+1,j = [Ω i+1,j,I , Ω i+1,j,Γ ]. The rest of the sets, restriction, and prolongation operators can be defined as given in section 1.

This manuscript is for review purposes only.

5. LSPSD matrices for levels strictly greater than 1. In [START_REF] Spillane | Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps[END_REF][START_REF] Dolean | An introduction to domain decomposition methods[END_REF][START_REF] Daas | A class of efficient locally constructed preconditioners based on coarse spaces[END_REF], different methods are suggested to obtain local SPSD splitting matrices at level 1. These matrices are used to construct efficient two-level additive Schwarz preconditioners.

Here in this section, we do not discuss the construction of these matrices at level 1. We suppose that we have the local SPSD matrices Ã1,j ∈ R n1×n1 for j = 1, . . . , N 1 . We focus on computing LSPSD matrices Ãi,j ∈ R ni×ni for (i, j) ∈ 2; L × 1; N i . We also suppose that the coarse space S 1 is available, i.e., the matrices V 1 and

A 2 = V 1 A 1 V 1
are known explicitly.

Proposition 5.1. Let i be a fixed level index, and let Ãi,j be an LSPSD of A i , (see Definition 3.1), associated with subdomain j, for j = 1, . . . , N i . Let G i,1 , . . . , G i,Ni+1 be a set of superdomains at level i associated with the partitioning at level i + 1, see subsection 4.1. Let V i be the restriction matrix to the coarse space at level i. Then, the matrix Ãi+1,j which is defined as:

Ãi+1,j = k∈Gi,j V i Ãi,k V i ,
satisfies the conditions in Definition 3.1. That is, Ãi+1,j is LSPSD of A i+1 with respect to subdomain j for j = 1, . . . , N i+1 . Proof. To prove that Ãi+1,j is LSPSD of A i+1 with respect to subdomain j, we have to prove the following:

• R i+1,j,∆ Ãi+1,j = 0 • u Ni+1 j=1 Ãi+1,j u k i+1 u A i+1 u for all u ∈ R ni+1 .
First, note that R i,k Ãi,j = 0 for all non-neighboring subdomains k of subdomain j. This yields Z i,k D i,k R i,k Ãi,j = 0 for these subdomains k. Now, let m ∈ 1; n i+1 \ Ω i+1,j . We will show that the mth row of Ãi+1,j is zero.

Following the partitioning of subdomains at level i + 1, there exists a subdomain Ω p0 such that the mth column of V i is part of R i,p0 D i,p0 Z i,p0 . We denote this column vector by v m . Furthermore, the subdomain p 0 is not a neighbor of any subdomain that is a part of the superdomain G i,j . Hence, v m Ãi,k = 0 for k ∈ G i,j . The mth row of Ãi+1,j is given as v m k∈Gi,j Ãi,k V i . Then, v m k∈Gi,j Ãi,k = 0, and the mth row of Ãi+1,j is zero.

To prove the second condition, we have

u Ni+1 j=1 Ãi+1,j u = u Ni+1 j=1 k∈Gi,j V i Ãi,k V i u.
Since {G i,j } 1 j Ni+1 form a disjoint partitioning of 1; N i , we can write

u Ni+1 j=1 Ãi+1,j u = u Ni k=1 V i Ãi,k V i u, = u V i Ni k=1 Ãi,k V i u.
This manuscript is for review purposes only.

Ãi,k is an LSPSD matrix of A i for k = 1, . . . , N i . Hence, we have

u Ni+1 j=1 Ãi+1,j u k i u V i A i V i u, k i u A i+1 u.
We finish the proof by setting k i+1 = k i .

A 2 = V 1 A 1 V 1 .
Note that column and row indices of A 2 are associated with column indices of V 1 . Four subdomains are used at level 2.

The partitioning at level 2 is related to the superdomain G 1,j = 4(j-1)+1; 4(j-1)+4 for j = 1, . . . , 4. Figure 5.1 (bottom right) represents an LSPSD matrix of A 2 with respect to subdomain 1 at level 2. Theorem 5.2 shows that the third condition of the fictitious subspace lemma Lemma 2.1 holds at level i for i = 1, . . . , L.

Theorem 5.2. Let Ãi,j be an LSPSD of A i associated with subdomain j, for (i, j) ∈ 1; L × 1; N i . Let τ i > 0, Z i,j be the subspace associated with Ãi,j , and P i,j be the projection on Z i,j as defined in Lemma 2.5. Let u i ∈ R ni and let u i,j = D i,j I ni,j -P i,j R i,j u i for (i, j) ∈ 1; L × 1; N i . Let u i,0 be defined as,

u i,0 = V i V i -1 V i   Ni j=1 R i,j D i,j P i,j R i,j u i   .
Let

m i = (2 + (2k i,c + 1)k i τ i ) -1 . Then, u i = Ni j=0 R i,j u i,j , and
(5.1) m i Ni j=0 u i,j R i,j A i R i,j u i,j u i A i u i .
Proof. We have

Ni j=0 R i,j u i,j = V i V i V i -1 V i   Ni j=1 R i,j D i,j P i,j R i,j u i   + Ni j=1 R i,j u i,j
This manuscript is for review purposes only.

                                 
A 2 = V 1 A 1 V 1 , bottom right: Ã2,1 = j∈G 1,1 V 1 Ã1,j V 1 , where G 1,1 = 1, . . . , 4
Since for all y ∈ S i ,

V i V i V i -1 V i y = y, we have Ni j=0 R i,j u i,j = Ni j=1 R i,j D i,j P i,j R i,j u i + Ni j=1 R i,j D i,j I ni,j -P i,j R i,j u i , = Ni j=1 R i,j D i,j R i,j u i , = u i .
To prove the inequality (5.1), we start with the inequality from Lemma 2.4. We have

Ni j=0 u i,j R i,j A i R i,j u i,j 2u i A i u i + (2k i,c + 1) Ni j=1 u i,j R i,j A i R i,j u i,j , (5.2)
where we chose u Bi in Lemma 2.4 to be (u i,j ) j=0,...,Ni and u Ai = u i . In Definition 3.2, This manuscript is for review purposes only.

we defined Z i,j , such that for all w ∈ R ni,j we have

(I ni,j -P i,j )w D i,j R i,j A i R i,j D i,j (I ni,j -P i,j )w τ i w R i,j Ãi,j R i,j w.
Hence, in the special case w = R i,j u i , we can write

(I ni,j -P i,j )R i,j u i D i,j R i,j A i R i,j D i,j (I ni,j -P i,j )R i,j u i τ i (R i,j u i ) R i,j Ãi,j R i,j (R i,j u i ).
Equivalently,

u i,j R i,j A i R i,j u i,j τ i (R i,j u i ) R i,j Ãi,j R i,j (R i,j u i ).
Plugging this inequality in (5.2) gives

Ni j=0 u i,j R i,j A i R i,j u i,j 2u i A i u i + (2k i,c + 1) τ i Ni j=1 (R i,j u i ) R i,j Ãi,j R i,j (R i,j u i ).
Since Ãi,j is local, we have

(R i,j u i ) R i,j Ãi,j R i,j (R i,j u i ) = u i Ãi,j u i , for j = 1, . . . , N i .
By using the fact that Ãi,j is LSPSD of A i for j = 1, . . . , N i , we obtain the following:

Ni j=0 u i,j R i,j A i R i,j u i,j 2u i A i u i + (2k i,c + 1) k i τ i u i A i u i .
Multiplying both sides with m i ends the proof, i.e.,

m i Ni j=0 u i,j R i,j A i R i,j u i,j u i A i u i .
In [START_REF] Daas | A class of efficient locally constructed preconditioners based on coarse spaces[END_REF], the authors presented the minimal subspace that replaces Z i,j (defined in (3.3) and used in Theorem 5.2) that is required to prove Theorem 5.2. The main difference with respect to the subspace that we define in (3.3) is that it is not necessary to include the entire kernel of the LSPSD matrix, K i,j , in Z i,j , see Definition 3.2. Nevertheless, in this work, we include the entire kernel of the LSPSD matrix in the definition of Z i,j . This allows us to ensure that the kernels of Neumann matrices are transferred across the levels, see Theorem 5.4. And in addition, this corresponds to the definition used in GenEO [START_REF] Dolean | An introduction to domain decomposition methods[END_REF]Lemma 7.7] and to its implementation in the HPDDM library [START_REF] Jolivet | Scalable domain decomposition preconditioners for heterogeneous elliptic problems[END_REF].

Theorem 5.3 provides an upper bound on the condition number of the preconditioned matrix M -1 i A i for i = 1, . . . , L.

Theorem 5.3. Let M i be the additive Schwarz preconditioner at level i combined with the coarse space correction induced by S i defined in (3.4). The following inequality holds,

κ M -1 i A i (k i,c + 1) (2 + (2k i,c + 1)k i τ i ) .
This manuscript is for review purposes only.

Proof. Lemma 2.2, Lemma 2.3, and Theorem 5.2 prove that the multilevel preconditioner verifies the conditions in Lemma 2.1 at each level i. Hence, the spectrum of the

preconditioned matrix M -1 i A i is contained in the interval [(2 + (2k i,c + 1)k i τ i ) -1 , k i,c + 1]
. Equivalently, the condition number of the preconditioned matrix at level i verifies the following inequality

κ M -1 i A i (k i,c + 1) (2 + (2k i,c + 1)k i τ i ) .
Proposition 5.1 shows that the constant k i associated with the LSPSD matrices at level i is independent of the number of levels and bounded by the number of subdomains at level 1. Indeed,

k 1 ≥ k i for i = 2, . . . , L.
Furthermore, in the case where the LSPSD matrices at the first level are the Neumann matrices, k i is bounded by the maximum number of subdomains at level 1 that share an unknown.

The constant k i,c for i = 1, . . . , L is the minimum number of distinct colors so that span{R i,j } 1 j Ni of the same color are mutually A i -orthogonal. Both constants k i and k i,c are independent of the number of subdomains for each level i.

The constant τ i can be chosen such that the condition number of the preconditioned system at level i is upper bounded by a prescribed value. Hence, this allows to have a robust convergence of the preconditioned Krylov solver at each level.

Algorithm 5.1 presents the construction of the multilevel additive Schwarz method by using GenEO. The algorithm iterates over the levels. At each level, three main operations are performed. First, the construction of the LSPSD matrices. At level 1, the LSPSD matrices are the Neumann matrices, otherwise, Proposition 5.1 is used to compute them. Once the LSPSD matrix is available, the generalized eigenvalue problem in (3.2) has to be solved concurrently. Given the prescribed upper bound on the condition number, Z i,j can be set. Finally, the coarse space is available and the coarse matrix is assembled.

The following Theorem 5.4, describes how the kernel of Neumann matrices are transferred across the levels.

Theorem 5.4. Suppose that Ã1,j is the Neumann matrix associated with the subdomain Ω 1,j for j ∈ 1; N 1 . For (i, j) ∈ 2; L × 1; N i , let

• Ãi,j be the LSPSD matrices associated with A i,j defined in Proposition 5.1,

• G i-1,j be the corresponding superdomains,

• G 1 i-1,j be the union of subdomains at level 1 which contribute hierarchically to obtain G i-1,j ,

• ÃGi-1,j be the Neumann matrix associated with

G 1 i-1,j (seeing G 1 i-1,j as a subdomain),
• A Gi-1,j be the restriction of A to the subdomain

G 1 i-1,j .
Then, the kernel of ÃGi-1,j is included in the kernel of

i-1 l=1 V l Ãi,j i-1 l=1 V l .
Proof. First, note that for any LSPSD matrix computed as in Proposition 5.1, we have

i-1 l=1 V l Ãi,j i-1 l=1 V l = i-1 l=1 V l i-1 l=1 V l k∈G 1 i,j Ã1,k i-1 l=1 V l i-1 l=1 V l .
Algorithm 5.1 Multilevel GenEO Require: A 1 = A ∈ R n×n SPD, L + 1 number of levels, N i number of subdomains at each level, G i,j sets of superdomains Ensure: preconditioner at each level i, M -1 i with bounded condition number of

M -1 i A i 1: for i = 1, . . . , L do 2:
for each subdomain j = 1, . . . , N i do 3:

A i,j = R i,j A i R i,j (local matrix associated with subdomain j)

4: if i = 1 then 5:
local SPSD Ãi,j is Neumann matrix of subdomain j 6:

else 7: compute local SPSD matrix as Ãi,j = k∈Gi,j V i-1 Ãi-1,k V i-1 8: end if 9:
solve the generalized eigenvalue problem (3.2), set Z i,j as in (3.3)

S i = Ni j=1 D i,j R i,j Z i,j , V i basis of S i 12: coarse matrix A i+1 = V i A i V i , A i+1 ∈ R ni+1×ni+1 13: end for 14: M -1 i = V i A -1 i+1 V i + Ni j=1 R i,j A -1 i,j R i,j
Moreover, due to the fact that ÃGi-1,j and Ã1,k are Neumann matrices, we have u ÃGi-1,j u u

k∈G 1 i,j Ã1,k u k 1 u ÃGi-1,j u.
On one hand, the kernels of Ã1,k for k ∈ G 1 i,j are included, by construction, in the image of V 1 , see Definition 3.2. So is their intersection which is the kernel of k∈G 1 i,j Ã1,k .

On the other hand, the previous two-sided inequality implies that the kernels of ÃGi-1,j and k∈G 1 i,j Ã1,k are identical. Hence, the kernel of ÃGi-1,j is included in the image of QQ , where

Q = i-1 l=1 V l .
Theorem 5.4 proves that the kernel of the Neumann matrix of a union of subdomains at level 1 that hierarchically contribute to form a subdomain at level i is conserved by the construction of the hierarchical coarse spaces. For example in the case of linear elasticity, it is essential to include the rigid body motions in the coarse space in order to have a fast convergence. As these are included in the kernel of the Neumann matrix of the subdomain, the hierarchical coarse space includes them, consequently.

6. Numerical experiments. In this section, the developed theory is validated numerically with FreeFEM [START_REF] Hecht | New development in FreeFem++[END_REF] for finite element discretizations and HPDDM [START_REF] Jolivet | Scalable domain decomposition preconditioners for heterogeneous elliptic problems[END_REF] for domain decomposition methods. We present numerical experiments on two highly challenging problems illustrating the efficiency and practical usage of the proposed method. For both problems, we use N 1 = 2,048 MPI processes (equal to the number of subdomains at level 1), and the domain partitioning is performed using ParMETIS

This manuscript is for review purposes only.

[22], with no control on the alignments of subdomain interfaces. We compare the two-level GenEO preconditioner and its multilevel extension by varying N 2 between 4 and 256. For the two-level method, N 2 corresponds to the number of MPI processes that solve the coarse problem in a distributed fashion using MKL CPARDISO [START_REF]MKL web page[END_REF].

For the multilevel method, N 3 is set to 1, i.e., a three-level method is used. The goal of these numerical experiments is to show that when one switches from a two-level method with an exact coarse solver, to our proposed multilevel method, the number of outer iterations is not impacted. Thus, three levels are sufficient. As an outer solver, since all levels but the coarsest are solved approximately, the flexible GMRES [START_REF] Saad | A Flexible Inner-Outer Preconditioned GMRES Algorithm[END_REF] is used. It is stopped when relative unpreconditioned residuals are lower than 10 -6 . Subdomain matrices {A i,j } 1 i 2,1 j Ni are factorized concurrently using MKL PARDISO, and eigenvalue problems are solved using ARPACK [START_REF] Lehoucq | ARPACK users' guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF]. In both, twoand three-level GenEO, we factorize the local matrices A 1,j for j ∈ 1; N 1 and solve the generalized eigenvalue problems concurrently at the first level. For this reason, we do not take into account the time needed for these two steps which are performed without any communication between MPI processes. We compare the time needed to assemble and factorize A 2 in the two-level approach against the time needed to assemble A 2 and local SPSD matrices Ã2,j for j ∈ 1; N 2 , solve the generalized eigenvalue problems concurrently on the second level, assemble, and factorize the matrix A 3 in the three-level approach. We also compare the time spent in the outer Krylov solver during the solution phase. Readers interested by a comparison of the efficiency of GenEO and multigrid methods such as GAMG [START_REF] Adams | Ultrascalable Implicit Finite Element Analyses in Solid Mechanics with over a Half a Billion Degrees of Freedom[END_REF] are referred to [START_REF] Jolivet | Domain decomposition methods[END_REF].

FreeFEM scripts used to produce the following results are available at the following URL: https://github.com/prj-/aldaas2019multi 1 . The exterior normal vector to the boundary of Ω is denoted n. Γ D is the subset of the boundary of Ω corresponding to x = 0 in 2D and 3D. Γ N is defined as the complementary of Γ D with respect to the boundary of Ω. We discretize the equation using P 2 and P 4 finite elements in the 3D and 2D test cases, respectively. The number of unknowns is 441 × 10 6 and 784 × 10 6 , with approximately 28 and 24 nonzero elements per row in the 3D and 2D cases, respectively. The heterogeneity is due to the jumps in the diffusion coefficient κ, see Figure 6.1, which is modeled using a combination of jumps and channels, cf. the file coefficients.idp from https: //github.com/prj-/aldaas2019multi.

The results in two dimensions are reported in Table 6.1. The number of outer iterations for both two-and three-level GenEO is 32. The size of the level 2 operator is n 2 = 25 × 2,048 = 51,200. In all numerical results, the number of eigenvectors per subdomain, here 25, is fixed. This is because ARPACK cannot a priori compute all eigenpairs below a certain threshold, and an upper bound has to be provided instead.

1 note to reviewers: the repository is now public For both two-and three-level methods, the second column is the time spent in the outer Krylov solver once the preconditioner has been set up. In the last column of the three-level method, the number of inner iterations for solving systems involving A 2 , which is not inverted exactly anymore, is reported. For all tables, this column is an average over all successive outer iterations. Another important numerical property of our method is that, thanks to fully controlled bounds at each level, the number of inner iterations is low, independently of the number of superdomains N 2 . Because this problem is not large enough, it is still tractable by a two-level method, for which HPDDM was highly optimized for.

Thus, there is no performance gain to be expected at this scale. However, one can notice that the construction of the coarse operator(s) scales nicely with N 2 for the three-level method, whereas the performance of the direct solver MKL CPARDISO quickly stagnates because of the finer and finer parallel workload granularity.

The results in three dimensions are reported in The physical domain Ω is a beam of dimensions [0, 10] × [0, 1], extruded for z ∈ [0, 1] in 3D. The Cauchy stress tensor σ(•) is given by Hooke's law: it can be expressed in terms of Young's modulus E and Poisson's ratio ν.

σ ij (u) = 2µε ij (u) i = j, 2µε ii (u) + λdiv(u) i = j,
where

ε ij (u) = 1 2
∂u i ∂x i + ∂u j ∂x j , µ = E 2(1 + ν)
, and λ = Eν 1 -2ν .

The exterior normal vector to the boundary of Ω is denoted n. Γ D is the subset of the boundary of Ω corresponding to x = 0 in 2D and 3D. Γ N is defined as the complementary of Γ D with respect to the boundary of Ω. We discretize (6.1) using the following vectorial finite elements: (P Results in two (resp. three) dimensions are reported in Table 6.3 (resp. Table 6.4).

The number of outer iterations are 73 and 45 respectively. For these test cases, we

This manuscript is for review purposes only. slightly relaxed the criterion for selecting eigenvectors in coarse spaces, which explains why the iteration counts increase. However, the same observations as for the diffusion test cases still hold. The dimension of the level 2 matrix is n 2 = 50×2,048 = 1.02•10 5 , while for the level 3 matrix it is n 3 = 20 × N 2 . This means that 50 (resp. [START_REF] Jones | AMGe Based on Element Agglomeration[END_REF] eigenvectors are kept per level 1 (resp. level 2) subdomains. We observe that the number of iterations of the inner solver increases slowly when increasing the number of subdomains from 4 to 256 in the 2D case and remains almost constant in the 3D case. In terms of runtime, the two-level GenEO is faster than three-level GenEO for these matrices of medium dimensions.

To show the potential of our method at larger scales, a three-dimensional linear elasticity problem of size 593 × 10 6 is now solved on N 1 = 16,384 processes and N 2 = 256 superdomains. With the two-level method, A 2 is assembled and factorized in 40.8 seconds. With the three-level method, this step now takes 35.1 seconds, see Table 6.5. There is a two iterations difference in the iteration count. Not taking into account the preconditioner setup, the problem is solved in 222.5 seconds in the two-level case and 90.1 seconds in the multilevel case. In this test case the cost of applying the two-level preconditioner on a given vector is approximately twice the cost of applying the multilevel variant. At this regime, it is clear that there are important gains for the solution phase. At even greater scales, gains for the setup phase are also expected. Moreover, another interesting fact to note regarding computation time is that the generalized eigenvalue problems solved concurrently at the first level to obtain V 1 actually represents a significant part of the total time of 377.6 seconds (resp. 244.8 seconds) with the two-(resp. three-)level method: 78.2 seconds. This cost can be reduced by taking a larger number of (smaller) subdomains, with the drawback of increasing the size of V 1 and thus A 2 . This drawback represents a clear bottleneck for the two-level method but is alleviated by using the three-level method, making it a good candidate for problems at greater scales.

Conclusion.

In this paper, we reviewed general properties of overlapping Schwarz preconditioners and presented a framework for its multilevel extension. We This manuscript is for review purposes only. 8. Acknowledgments. We would like to thank the anonymous referees for their comments and remarks that helped us improve the clarity of this manuscript. This work was granted access to the HPC resources of TGCC@CEA under the allocation

Fig. 1 .

 1 Fig.1.1. Left: a triangular mesh is used to discretize the unit square. The set of nodes of the mesh is partitioned into 16 disjoint subsets, non-overlapping subdomains, Ω 1,j,I for j = 1, . . . , 16. Right: Illustration of the matrix A 1 whose connectivity graph corresponds to the mesh on the left. The diagonal block j of A 1 corresponds to the non-overlapping subdomain Ω 1,j,I . Each submatrix A 1 (Ω 1,j,I , Ω 1,j,I ) is colored with a distinct color. The same color is used to color the region of the square that contains nodes in Ω 1,j,I .

4. 1 .

 1 Superdomains as unions of several subdomains. In this section, we introduce the notion of a superdomain. It refers to the union of several neighboring subdomains. Let G i,1 , . . . , G i,Ni+1 be disjoint subsets of 1; N i , where Ni+1 j=1 G i,j = 1; N i . We call the union of the subdomains {k ∈ 1; N i : k ∈ G i,j } superdomain j, for j = 1, . . . , N i+1 . Figure4.1 gives an example of how to set superdomains. Though this definition of superdomains may look somehow related to the fine mesh, it is in practice done at the algebraic level, as explained later on. Note that the indices of columns and rows of A i+1 are associated with the vectors contributed by the subdomains at level i in order to build the coarse space S i , see Figure4.2. Hence, defining subdomains on the structure of A i+1 is natural once we have the subsets G i,j , for j = 1, . . . , N i+1 .

Ω 1 1 G 1 , 2 G 1 , 3 G 1 , 4 Fig. 4 . 1 .

 1112131441 Fig. 4.1. Left: 16 subdomains at level 1. Right: 4 superdomains at level 1. G 1,j = 4(j -1) + 1; 4(j -1) + 4 .

Figure 5 .

 5 Figure 5.1 gives an illustration of the LSPSD construction provided by Proposition 5.1.Figure 5.1 (top left) represents the matrix A 1 . The graph of A 1 is partitioned

Figure 5 . 1 (

 51 top left) represents the matrix A 1 . The graph of A 1 is partitioned into 16 subdomains. Each subdomain is represented by a different color.

Figure 5 . 1 (

 51 top right) represents the matrix V 1 whose column vectors form a basis of the coarse space S 1 . Colors of columns of V 1 correspond to those of subdomains in A 1 . Figure 5.1 (bottom left) represents the matrix

Fig. 5 . 1 .

 51 Fig. 5.1. Illustration of the LSPSD construction provided by Proposition 5.1. Top left: the matrix A 1 , top right: V 1 , bottom left: the matrixA 2 = V 1 A 1 V 1 , bottom right: Ã2,1 =

6. 1 .

 1 Diffusion test cases. The scalar diffusion equation with highly heterogeneous coefficient κ is solved in [0, 1] d (d = 2 or 3). The strong formulation of the equation is: -∇ • (κ∇u) = 1 in Ω, u = 0 on Γ D , ∂u ∂n = 0 on Γ N .

1 5 • 10 5 1 • 10 6 1. 7 • 10 6 κFig. 6 . 1 .

 1567661 Fig. 6.1. Variation of the coefficient κ used for the diffusion test case

6. 2 .

 2 Linear elasticity test cases. The system of linear elasticity with highly heterogeneous elastic moduli is solved in 2D and 3D. The strong formulation of the equation is given as: div σ(u) + f = 0 in Ω, u = 0 on Γ D , σ(u) • n = 0 on Γ N . (6.1)

Fig. 6 . 2 .

 62 Fig. 6.2. Variation of the structure coefficients used for the elasticity test casetwo-level GenEO three-level GenEO N 2 CS solve % of nnz A 2 CS solve inner it. % of nnz A 3 4 4.8 52.7 0.18 22.5 179.3 31 43.0 16 3.9 50.3 9.3 124.9 57 17.0 64 4.0 53.1 7.2 71.5 34 4.9 256 4.8 63.2 6.8 71.2 44 1.4

  generalized the local SPSD splitting presented in[START_REF] Daas | A class of efficient locally constructed preconditioners based on coarse spaces[END_REF] to cover a larger set of matrices leading to more flexibility for building robust coarse spaces. Based on local SPSD matrices on the first level, we presented how to compute local SPSD matrices for coarser levels. The multilevel solver based on hierarchical local SPSD matrices is robust and guarantees a bound on the condition number of the preconditioned matrix at each level depending on predefined values. Numerical experiments illustrate the theory and prove the efficiency of the method on challenging problems of large size arising from heterogeneous linear elasticity and diffusion problems with jumps in the coefficients of multiple orders of magnitude.

Table 6

 6 CS solve % of nnz A 2 CS solve inner it. % of nnz A 3

	.2. The number of outer

in two dimensions still hold, and the dimensions of A 2 and A 3 are the same. Once again, it is important to note that the number of outer iterations is the same for both methods.

This manuscript is for review purposes only.

Table 6

 6

	.1 Diffusion 2D test case, comparison between two-and
	three-level GenEO. The percentage of nonzero entries
	in A 1 is 0.3%.
	two-level GenEO		three-level GenEO
	N 2 CS solve % of nnz A 2	CS solve inner it. % of nnz A 3
	4 7.0 20.9		16.9 43.6	17	62.0
	16 5.0 19.8 64 5.1 20.1	0.36	7.7 5.8	26.7 32.7	17 15	28.0 8.9
	256 5.2 24.1		5.3	22.6	14	2.6

Table 6 . 2

 62 Diffusion 3D test case, comparison between two-and three-level GenEO. The percentage of nonzero entries in A 1 is 0.5%.

  2 , P 2 , P 2 ) in 3D and (P 3 , P 3 ) in 2D. The number of unknowns is 146 × 10 6 and 847 × 10 6 , with approximately 82 and 34 nonzero elements per row in the 3D and 2D cases, respectively. The heterogeneity is due to the jumps in E and ν. We consider discontinuous piecewise constant values

for E and ν: (E 1 , ν 1 ) = (2 × 10 11 , 0.25), (E 2 , ν 2 ) = (10 7 , 0.45), see Figure

6

.2.

Table 6 . 3

 63 Elasticity 2D test case, comparison between two-and three-level GenEO. The percentage of nonzero entries in A 1 is 0.4%.

Table 6 . 4

 64 two-level GenEO three-level GenEO N 2 CS solve % of nnz A 2 CS solve inner it. % of nnz A 3 4 28.5 46.9 Elasticity 3D test case, comparison between two-and three-level GenEO. The percentage of nonzero entries in A 1 is 3.3%.

			78.9 296.7	23	43.0
	16 17.3 35.4 64 15.0 33.2	0.38	24.5 124.5 15.4 62.2	23 21	19.0 7.9
	256 13.6 40.7		10.6 50.7	23	2.5
	two-level GenEO	three-level GenEO
	N 2 CS	solve		CS solve inner it.
	256 40.8	222.5		35.1 90.1	11

Table 6 . 5

 65 Elasticity 3D test case, comparison between two-and three-level GenEO

This manuscript is for review purposes only.

A0050607519 made by GENCI. The work of the second author was supported by the NLAFET project as part of European Union's Horizon 2020 research and innovation program under grant 671633.