# A scalable Domain Decomposition method for Saddle Point problems with GenEO coarse spaces

Filippo Brunelli, Guillaume Delay, Frédéric Nataf, Emile Parolin, Pierre-Henri Tournier

<figure><img src="image_1.png" /></figure>

To cite this version:

Filippo Brunelli, Guillaume Delay, Frédéric Nataf, Emile Parolin, Pierre-Henri Tournier. A scalable Domain Decomposition method for Saddle Point problems with GenEO coarse spaces. 2025. hal-05123311

HAL Id: hal-05123311

https://hal.science/hal-05123311v1

Preprint submitted on 20 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# A scalable Domain Decomposition method for Saddle Point problems with GenEO coarse spaces

Filippo Brunelli*<sup>1,2</sup>, Guillaume Delay<sup>1,3</sup>, Frédéric Nataf<sup>1,2</sup>, Emile Parolin<sup>§,1,2</sup>, and Pierre-Henri Tournier¶<sup>1,2</sup>

<sup>1</sup>Laboratory J.L. Lions (UMR 7598), Sorbonne University

<sup>2</sup>Alpines Team, INRIA Paris

<sup>3</sup>Commedia team, INRIA Paris

June 20, 2025

## Abstract

We present an adaptive domain decomposition (DD) preconditioning technique for the solution of saddle point problems with a 2x2 blocks structure. This work utilises the GenEO theory for symmetric positive definite (SPD) problems (Spillane et al., 2014) to tackle the solution of saddle point linear systems with iterative methods. The latters are preferred over direct methods due to the typically large size of these problems. Under the same assumptions of the GenEO method we can develop a robust and scalable preconditioning technique for the saddle point linear system. We present a preconditioner for the primal-primal block and one for the Schur complement matrix, both based on the DD paradigm. Then, these ingredients are used to precondition the 2x2 block saddle point matrix. Numerical tests are performed on a classical systems of partial differential equations (PDEs) discretized with finite elements, for which the assumptions are easy to verify. The results emphasize preconditioner's robustness with respect to the specific domain's partitioning.

<footnote>* filippo.brunelli@sorbonne-universite.fr</footnote>

<footnote>† guillaume delaysorbonne-universite.fr</footnote>

<footnote>‡ frederic.nataf@sorbonne-universite.fr</footnote>

<footnote>§ emile.parolin@inria.fr</footnote>

<footnote>¶ pierre-henri.tournier@sorbonne-universite.fr</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>3</td></tr><tr><td>2</td><td>Preconditioning of the primal problem</td><td>4</td></tr><tr><td>2.1</td><td>Schwarz DD preconditioners</td><td>4</td></tr><tr><td>2.2</td><td>Fictitious Space Lemma (FSL)</td><td>5</td></tr><tr><td>2.3</td><td>GenEO method for Schwarz preconditioners</td><td>6</td></tr><tr><td>2.4</td><td>Bounds for the condition numbers</td><td>7</td></tr><tr><td>3</td><td>Schur complement preconditioning</td><td>8</td></tr><tr><td>3.1</td><td>Approximation of the Schur complement matrix</td><td>9</td></tr><tr><td>3.2</td><td>Two-level preconditioner for the sum of local Schur complement matrices</td><td>12</td></tr><tr><td>3.3</td><td>Preconditioner for the Schur complement matrix</td><td>14</td></tr><tr><td>3.3.1</td><td>Spectral bounds for the Schur complement preconditioned system</td><td>14</td></tr><tr><td>3.3.2</td><td>Application of the preconditioner with nested iterative solution</td><td>15</td></tr><tr><td>3.3.3</td><td>New way to apply the preconditioner using an algebraic formula</td><td>16</td></tr><tr><td>3.3.4</td><td>Practical aspects for the preconditioner application</td><td>16</td></tr><tr><td>4</td><td>Saddle point preconditioning</td><td>17</td></tr><tr><td>5</td><td>Numerical experiments</td><td>20</td></tr><tr><td>5.1</td><td>Heterogeneous linear elasticity problem</td><td>20</td></tr><tr><td>5.2</td><td>Hardware, softwares and implementation details</td><td>21</td></tr><tr><td>5.3</td><td>Weak scalability experiment for heterogeneous beam</td><td>23</td></tr><tr><td>5.3.1</td><td>Experiment using the method based on the Woodbury formula</td><td>23</td></tr><tr><td>5.3.2</td><td>Experiment using the method based on nested iterative solution</td><td>24</td></tr><tr><td>5.4</td><td>Comparison with norm based preconditioner</td><td>25</td></tr><tr><td>6</td><td>Conclusion and outlooks</td><td>26</td></tr><tr><td></td><td>References</td><td>26</td></tr></table>
# 1 Introduction

Saddle point problems arise in a multitude of applications throughout scientific computing such as solid and fluid mechanics, computational electromagnetism, inverse problems and optimization. In particular, discretization of partial differential equations (PDEs) using mixed finite elements and several numerical methods for constrained optimization problems quickly became popular sources of large-scale saddle point linear systems. For a comprehensive list of applications and contexts where saddle point problems arise, see [4] and references therein.

Solving large saddle point linear systems typically proves to be a difficult challenge. Indeed, due to the large size of the problems and their often poor spectral properties, both direct and iterative method struggles. The formers have large memory requirements, while the latter are affected by indefiniteness and ill-conditioning, which are usual characteristics of problems from this class. However, when working with sparse saddle point systems, which is the case in most of the fields above, iterative methods are preferred over direct ones. Consequently, the derivation of suitable preconditioners becomes a crucial task in order to improve the convergence of such methods. The choice of a preconditioner for a saddle point system is known to be highly problem-dependent. This is why in this work we aim to present a method to address the widest possible range of problems.

Iterative methods for large sparse linear systems can benefit from the increasing availability of massively parallel architectures. Accordingly, domain decomposition (DD) methods became one of the dominant paradigms to combine preconditioning and parallelism, alongside with multi-grid methods. The interest in them comes from the possibility of defining efficient and robust preconditioners that are naturally well-fitted to modern parallel machines.

For symmetric positive definite (SPD) systems, two-level Schwarz DD methods using adaptive spectral coarse spaces have been successfully developed, see [8, 16, 22], book [6] and references therein. For specific systems of equations with saddle point structure, efficient DD preconditioners have been designed, e.g [19, 20], book [23] and references therein. More recently, a method to tackle the solution of generic saddle point linear systems combining different DD strategies is proposed in [15]. However, a nested iterative solution is needed at every application of the proposed preconditioner and the iteration count of the inner solution step is not under theoretical control, see Section 3.3.2. In the present work we provide a new way to apply the preconditioner developed in [15] without the need of a nested iterative solution. In this way we can benefit from the theoretical control that we have on the iteration count for the outer preconditioner.

Let us consider the linear system

$$Au = b, \qquad (1.1)$$

where

$$A := \begin{bmatrix} A & B^T \\ B & -C \end{bmatrix}, \quad u := \begin{bmatrix} x \\ y \end{bmatrix}, \quad \text{and} \quad b := \begin{bmatrix} f \\ g \end{bmatrix}.$$

Here, $A \in \mathbb{R}^{n \times n}$ is the primal-primal block, $B \in \mathbb{R}^{m \times n}$, with $m < n$, is the constraint matrix and $C \in \mathbb{R}^{m \times m}$ is the dual-dual block. Unknowns $x \in \mathbb{R}^n$ and $y \in \mathbb{R}^m$ are the primal and dual variables respectively. Moreover $f \in \mathbb{R}^n$ and $g \in \mathbb{R}^m$ are the right hand sides for the primal and dual variables respectively. Let us assume that $A$ is SPD, that $C$ is symmetric positive semidefinite (SPSD) and that matrices $B$ and $C$ fulfill the assumption $\mathrm{Ker}(B^T) \cap \mathrm{Ker}(C) = \{0\}$. In particular, the case where $C = 0$ is allowed if $B$ is full rank. These are necessary and sufficient conditions for the saddle point matrix $A$ to be symmetric and nonsingular, see [4] for more general properties. Moreover, let us assume that all matrices are sparse. Indeed, sparsity allows to derive efficient DD preconditioners. If $A$ is nonsingular, we have the following factorization of
A,

$$A = \begin{bmatrix} A & B^T \\ B & -C \end{bmatrix} = \begin{bmatrix} I & 0 \\ BA^{-1} & I \end{bmatrix} \begin{bmatrix} A & 0 \\ 0 & -S \end{bmatrix} \begin{bmatrix} I & A^{-1}B^T \\ 0 & I \end{bmatrix}, \quad (1.2)$$

where we denote by $S := C+BA^{-1}B^T$ the Schur complement matrix of $A$. Matrix $S$ is symmetric and also positive definite under the previous assumptions but it is dense and typically very ill-conditioned. From factorization (1.2) it is evident that $A$ is nonsingular if and only if $S$ is also nonsingular. Finally, matrix $A$ is also indefinite with $n$ positive and $m$ negative eigenvalues.

The iterative resolution of linear systems like (1.1) can be performed in several ways. One approach is to split it into steps and alternate updates of the unknowns $x$ and $y$. This consists in alternating solutions of linear systems for matrices $A$ and $S$. An algorithm of this kind belongs from the class of segregated iterations algorithms or Schur complement approaches. These algorithms are usually derived from some factorizations of $A$, like (1.2). If the problems for $A$ and $S$ are large as well, preconditioning must be introduced to solve them iteratively. A different approach is solving the original problem (1.1) with an iterative method. In this case, block preconditioners are used to correct the indefiniteness of $A$ or to provide a more favorable eigenvalue distribution. They are typically derived from decompositions of $A$ as well, by replacing the blocks $A$ and $S$ with some approximations.

Both ways require preconditioning for the primal-primal block $A$ and for the Schur complement $S$. Working with segregated iterations allows to analyze the individual performances of preconditioners for $A$ and $S$ introduced at every step. On the other hand, one big disadvantage is that convergence can be assessed only by looking at the individual steps. Instead, an iterative method on the original saddle point system allows to assess convergence by looking at the residual on the global unknowns vector.

In [15], a two-level additive Schwarz preconditioner for the primal-primal block $A$ is presented, followed by the derivation of a two-level DD preconditioner for the Schur complement matrix $S$. A segregated iteration algorithm is considered to individually check the performances of the presented preconditioners. Then, a block preconditioning technique is introduced as well. Following the same argument of [15], Section 2 and 3 are devoted to define suitable preconditioners for $A$ and $S$. We will give the details on the argument and present our contribution at the end of Section 3. Then, in Section 4 the preconditioner developed in the previous sections will be used to build a block preconditioner for (1.1). Finally, is Section 5 some numerical results will be presented while conclusions and some further directions of research will be outlined in Section 6.

# 2 Preconditioning of the primal problem

In this section we aim to introduce a two-level DD preconditioner for the primal-primal block matrix $A$. We will start by presenting the algebraic setting to define a Schwarz DD preconditioner and some well-known choices of two-level Schwarz preconditioners. Then we will briefly present the GenEO method to define a suitable coarse space for SPD problems [22] and we will end by reporting some spectral results for the preconditioned matrices

## 2.1 Schwarz DD preconditioners

Let us denote $\mathcal{N}$ the global set of degrees of freedom (dofs) for the primal variable, with $n := \#\mathcal{N}$. Let us decompose $\mathcal{N}$ into $N$ overlapping subdomains, denoted by $N_i$, with $n_i := \#N_i$. This can be done with a graph partitioning algorithm on the connectivity graph of matrix $A$. For $1 \le i \le N$, let us define the extension with zero matrix $R_i^T : \mathbb{R}^{n_i} \to \mathbb{R}^n$, the associated restriction
matrix $R_i$. Moreover let us define a family of diagonal matrices $D_i \in \mathbb{R}^{n_i \times n_i}$, for $1 \le i \le N$, called partition of unity (PoU), such that

$$I_n = \sum_{i=1}^{N} R_i^T D_i R_i.$$

We aim to approximate the inverse of A with the resolution of local problems on the subspaces defined by matrices $R_i$ and on a subspace $V_0$ of $\mathbb{R}^n$, called coarse space, for which a definition will be given later. There exist several possibilities to define a two-level Schwarz preconditioner, see Chapter 2 of [23] for a presentation of the general theory about Schwarz DD methods. One possibility is the two-level additive Schwarz (AS) preconditioner, which reads

$$M_{AS,2}^{-1} := R_0^T (R_0 A R_0^T)^{-1} R_0 + \sum_{i=1}^{N} R_i^T (R_i A R_i^T)^{-1} R_i \quad (2.1) \\ = R_0^T (R_0 A R_0^T)^{-1} R_0 + M_{AS,1}^{-1}.$$

Let $n_0 := \dim V_0$ and let $R_0^T : \mathbb{R}^{n_0} \to \mathbb{R}^n$ be the extension with zeros matrix from $V_0$ to $\mathbb{R}^n$. It is defined such that its columns form a basis for the coarse space $V_0$. The first term in (2.1) corresponds the second level, while $M_{AS,1}^{-1}$ is the one-level AS preconditioner. See Chapter 4 in [6] for an explanation of why the resolution of a coarse problem is fundamental to achieve faster convergence.

Other popular formulae to integrate the coarse correction to the one level preconditioner are the hybrid Schwarz method (HSM) and the multiplicative method, see Chapter 2 in [23]. In particular, the HSM formulation reads

$$M_{HSM}^{-1} := R_0^T (R_0 A R_0^T)^{-1} R_0 + (I - P_0) M_{AS,1}^{-1} (I - P_0^T), \quad (2.2)$$

where $P_0 := R_0^T (R_0 A R_0^T)^{-1} R_0 A$ is the A-orthogonal projection from $\mathbb{R}^n$ on $V_0$.

## 2.2 Fictitious Space Lemma (FSL)

In this section we will introduce a very useful tool to prove spectral bounds for the preconditioned matrix within the DD framework. This Lemma is frequently encountered in the DD literature, see [17] for the original paper and [10] for a more modern version of the statement. It can be also related to the general theory of Schwarz methods from [23]. It is used to have a result on the approximation quality that one gets using an operator defined on the decomposition of the original space as an approximation of some operator. Let us state the Lemma.

**Lemma 2.1 (FSL).** Let $H$, $H_D$ be two Hilbert spaces with some associated scalar products $(\cdot, \cdot)$, $(\cdot, \cdot)_D$. Let $a: H \times H \to \mathbb{R}$ and $b: H_D \times H_D \to \mathbb{R}$ be two symmetric, continuous, coercive bilinear forms. Moreover let $A: H \to H'$ and $B: H_D \to H'_D$ be differential operators defined as

$$\langle Au, v \rangle := a(u, v), \quad \forall u, v \in H, \quad \langle Bu_D, v_D \rangle := b(u_D, v_D), \quad \forall u_D, v_D \in H_D.$$

Finally, let $\mathcal{R}: H_D \to H$ be a continuous linear operator. Assume that the following assumptions hold

1. $\mathcal{R}$ is surjective.

2. (Continuity). $\exists c_R > 0$ such that $\forall u_D \in H_D$,

$$a(\mathcal{R}u_D, \mathcal{R}u_D) \le c_R b(u_D, u_D).$$
3. (Stable decomposition). $\exists c_T > 0$ such that $\forall u \in H$, $\exists u_D \in H_D$ with $Ru_D = u$ and

$$c_T b(u_D, u_D) \leq a(u, u).$$

Let us define the adjoint operator $R^* : H \to H_D$ by $(Ru_D, u) = (u_D, R^*u)_D$, $\forall u_D \in H_D$, $\forall u \in H$. Then the following spectral estimate holds $\forall u \in H$

$$c_T a(u, u) \leq a(RB^{-1}R^*Au, u) \leq c_R a(u, u). \qquad (2.3)$$

## 2.3 GenEO method for Schwarz preconditioners

The GenEO method [22] consists in a way of designing the space $V_0$ and the related extension matrix $R_0^T$ so that the condition number of the preconditioned system is independent of the specific decomposition and in particular from the global number of subdomains. The underlying idea is based on the solution of some generalized eigenvalue problems (GEPs) locally on every subdomain. Indeed, we look for $(z_{ik}, \lambda_{ik}) \in \mathbb{R}^{n_i} \times \mathbb{R}$ such that

$$D_i R_i A R_i^T D_i z_{ik} = \lambda_{ik} A'_i z_{ik}. \qquad (2.4)$$

In this way we are able to identify the modes responsible for the slow convergence of the one-level preconditioner which is purely local, see Chapter 4 in [6]. Matrices $A'_i \in \mathbb{R}^{n_i \times n_i}$ fulfill the following assumption.

**Assumption 2.2.** There exist a family of SPSD matrices $(A'_i)_{1 \leq i \leq N}$ that defines a stable SPSD splitting of A. Namely there exists a fixed constant $k_A > 0$ such that $\forall x \in \mathbb{R}^n$

$$(\sum_{i=1}^{N} R_i^T A'_i R_i x, x) \leq k_A (Ax, x).$$

This assumption is fundamental for the GenEO method. Indeed, the definition of such matrices $(A'_i)_{1 \leq i \leq N}$ is the key step to show the stable decomposition assumption of Lemma 2.1. In the finite element context, $(A'_i)_{1 \leq i \leq N}$ are called Neumann matrices and they are typically easy to define, see Section 5.2, or Lemma 7.13 in [6]. For more details on the definition of stable SPSD splitting, see [1]. Let us consider $\tau_A > 0$ user-defined parameter and let us define

$$Z_i := \text{Span}\{z_{ik} \text{ solutions of (2.4)} | \lambda_{ik} > \tau_A\}, \quad 1 \leq i \leq N.$$

The coarse space $V_0$ is defined as

$$V_0 := \bigoplus_{i=1}^{N} R_i^T D_i Z_i, \qquad (2.5)$$

where the abuse of notation $R_i^T D_i Z_i$ means applying the PoU and the extension to every vector in the basis of $Z_i$.

Let us define some recurrently appearing parameters. Let

$$k_0 := \max_{1 \leq i \leq N} \#\{1 \leq j \leq N | R_i A R_j^T \neq 0\}, \quad k_0 > 1,$$

be the maximum number of interactions that one subdomains has with his neighbors and

$$k_1 := \max_{k \in \mathcal{N}} \#\{1 \leq i \leq N | k \in \mathcal{N}_i\}, \quad k_1 \geq 1, \qquad (2.6)$$

be the maximum multiplicity of intersections between subdomains. These two parameters depend neither on the global number of subdomains $N$ nor on the mesh parameter. Typically, in the finite element contexts the parameter from Assumption 2.2 takes the value $k_A = k_1$.
## 2.4 Bounds for the condition numbers

Let us state one definition that we will use throughout the paper.

**Definition 2.3.** Two SPD operators $P_1, P_2$ on $\mathbb{R}^k$ are said to be spectrally equivalent with some constants $c_1, c_2$ if

$$ (P_1x, x) \le (P_2x, x) \le (P_1x, x), \quad \forall x \in \mathbb{R}^k. $$

Considering the two-level AS preconditioner (2.1), thanks to Lemma 2.1, it can be proved that $M_{AS,2}^{-1}$ and $A^{-1}$ are spectrally equivalent with some constants that depend only on the parameters defined in the previous section. Following what is done in Chapter 7 of [6], one can set the two-level AS preconditioner (2.1) into the FSL framework and obtain that

$$ \lambda_{m,AS}(Ax, x) \le (AM_{AS,2}^{-1}Ax, x) \le \lambda_{M,AS}(Ax, x), \quad \forall x \in \mathbb{R}^n, $$

with

$$ \lambda_{m,AS} := \frac{1}{2 + (2k_0 + 1)k_1\tau_A}, \quad \lambda_{M,AS} := 2k_0. $$

From this, the spectral equivalence between $A^{-1}$ and $M_{AS,2}^{-1}$ follows straightforwardly. Indeed, it is enough to consider $x = A^{-1}\hat{x}$ to obtain (after renaming all the $\hat{x}$ to $x$)

$$ \lambda_{m,AS}(A^{-1}x, x) \le (M_{AS,2}^{-1}x, x) \le \lambda_{M,AS}(A^{-1}x, x), \quad \forall x \in \mathbb{R}^n. \quad (2.7) $$

**Lemma 2.4.** Let $P_1, P_2$ be two SPD operators on $\mathbb{R}^k$. If $\exists c > 0$ such that

$$ (P_2x, x) \le c(P_1x, x), \quad \forall x \in \mathbb{R}^k, $$

Then $P_1^{-1}P_2$ has real positive eigenvalues bounded from above by $c$.

Recalling that $A$ and $M_{AS,2}$ are SPD, (2.7) and Lemma 2.4 imply that every eigenvalue of $M_{AS,2}^{-1}A$ is real, positive and moreover

$$ \lambda_{m,AS} \le \lambda(M_{AS,2}^{-1}A) \le \lambda_{M,AS}. \quad (2.8) $$

**Definition 2.5.** An operator $P^{-1}$ is said to be a robust preconditioner for an SPD operator $T$ if $P^{-1}$ is spectrally equivalent to $T^{-1}$ with constants that are independent of mesh parameters or problem coefficients.

**Definition 2.6.** An operator $P^{-1}$ is said to be a scalable DD preconditioner for an SPD operator $T$ if $P^{-1}$ is spectrally equivalent to $T^{-1}$ with constants that are independent of the number of subdomains $N$.

The two-level AS preconditioner using the GenEO method is a robust and scalable preconditioner for $A$. Indeed, constants $\lambda_{m,AS}, \lambda_{M,AS}$ depend only on local features of the decomposition and are independent of the mesh parameter $h$ and from the problem coefficients.

**Remark 2.7.** In practical applications, scalability may also refer to the computational time needed to have a solution. In Definition 2.6, the focus is on the eigenvalue distribution, which affects the number of iterations required by an iterative solver to converge. However, in order to get a competitive method also the setup time to build the preconditioner must be taken into consideration. In particular, it is important to ensure that the coarse space is not very large, so that the coarse problem is not too large to be factorized and solved with a direct method. Moreover, it is also important that the coarse space is made with extensions by zero of local vectors. This allows us to build the coarse space problem in an efficient way. In particular, the latter assumption is verified for the GenEO coarse space, see (2.5).
**Remark 2.8.** It is worth to notice that taking a smaller value for parameter $\tau_A$ means on one side to have a stricter bound for the smallest eigenvalue, hence possibly faster convergence, but on a larger coarse space on the other side. The choice of $\tau_A$ must take both into consideration.

Alternatively, formulation (2.2) can be used to incorporate the GenEO coarse space in the preconditioner. Using Lemma 2.1 we can prove some spectral bounds leading to the spectral equivalence between $M_{HSM}^{-1}$ and $A^{-1}$, see Chapter 7.5 of [6]. In particular we have that the eigenvalues of $M_{HSM}^{-1}A$ are such that

$$ \lambda_{m,HS} \le \lambda(M_{HSM}^{-1}A) \le \lambda_{M,HS}, \quad \lambda_{m,HS} := \frac{1}{1 + k_1\tau_A}, \quad \lambda_{M,HS} := k_0. \quad (2.9) $$

Interestingly, HSM preconditioner provides a sharper bound for the smallest eigenvalue. It is expected to bring faster convergence but there is one more synchronization needed due to the presence of the projector $P_0$.

## 3 Schur complement preconditioning

In this section we aim to derive a robust and scalable preconditioner for the Schur complement matrix $S$. This proves to be delicate since it means finding an approximation of $S$ that is cheap to invert and spectrally equivalent with some suitable constants.

Usually, information about the specific problem are exploited in order to approximate the Schur complement with something simpler, see Section 10.1.3 in [4] for an extensive survey of approximations of $S$. More recently, the theory of "natural" or norm-based block-diagonal preconditioners has been unified in an application-independent way in [21]. In this framework, an approximation for the Schur complement matrix $S$ is given by the Gram matrix representing some norm on the dual discrete functional space. If a discrete inf-sup condition is proved using such norm, then the corresponding Gram matrix is spectrally equivalent to the Schur complement. One example is the mass matrix that is spectrally equivalent to $S$ in the case of linear elasticity or steady Stokes discretized using inf-sup stable finite elements spaces. Unfortunately, of practical interest are the constants entering such spectral equivalence. They depend on the coercivity and inf-sup condition constants associated to the problem, which in turn depend on the coefficients and parameters of the problem at hand. It will be shown in Section 5.4 that for linear elasticity with highly heterogeneous coefficients, using the mass matrix as a preconditioner requires more than 500 iterations for the solution to converge. One attempt to lower this iteration count is to consider the rescaled mass matrix, which instead proves to be not robust.

In order to develop a method that can be used for general problems, something more sophisticated than the mass matrix must be derived. In particular, the goal is to find an operator that is spectrally equivalent to $S$, with some constants under control. To achieve this goal, we will exploit the spectral equivalences between $A$ and the Schwarz DD preconditioners defined in the previous section.

In the following, we are going to define some operators that are spectrally equivalent to $S$. We will then prove some conditions of the form (2.7) from which a result on the condition number follows directly. The following argument was first developed in [15] and here we give more details about it. Finally, we conclude by presenting a new way, with respect to [15], to apply the preconditioner developed for the Schur complement matrix.
## 3.1 Approximation of the Schur complement matrix

A first approximation of $S$ is obtained by replacing $A^{-1}$ with $M_{AS,2}^{-1}$ from (2.1) in the definition of the Schur complement, that is

$$ P_S := C + B M_{AS,2}^{-1} B^T. \qquad (3.1) $$

This is strictly related to constraint preconditioning (CP), see [3] or Chapter 10.2 in [4] for a survey. The underlying idea of CP is to approximate $A$ by replacing only the primal-primal block $A$ with some SPD approximation $P_A$, that is

$$ P_c := \begin{bmatrix} P_A & B^T \\ B & -C \end{bmatrix}. \qquad (3.2) $$

The application of $P_c^{-1}$ involves the solution of a linear system for the matrix $C + B P_A^{-1} B^T$, hence $P_A$ is chosen so that such linear system is easy to solve, e.g. $P_A = \text{diag}(A)$. However, this is usually not satisfactory in term of approximation, so more sophisticated $P_A$ are considered. But then, an approximation of $C + B P_A^{-1} B^T$ that is easier to use is needed. With $P_S$ from (3.1) we are in the latter case.

Let us denote $\mathcal{M}$ the global set of degrees of freedom for the dual variable, with $m := \#\mathcal{M}$. We claim that operator $P_S$ is spectrally equivalent to $S$. Indeed, thanks to (2.7) and the fact that $0 < \lambda_{m,AS} < 1 < \lambda_{M,AS}$, we can show that

$$ \frac{1}{\lambda_{M,AS}}(P_{Sy}, y) \le (Sy, y) \le \frac{1}{\lambda_{m,AS}}(P_{Sy}, y), \quad \forall y \in \mathbb{R}^m, $$

that from now on we will shorten with

$$ \frac{1}{\lambda_{M,AS}} P_S \le S \le \frac{1}{\lambda_{m,AS}} P_S. \qquad (3.3) $$

This can be easily proved by noticing that (2.7) holds $\forall x \in \mathbb{R}^n$, so in particular for $x = B^T y$, for some $y \in \mathbb{R}^m$. Writing $P_S$ explicitly we get

$$ P_S = C + B R_0^T (R_0 A R_0^T)^{-1} R_0 B^T + \sum_{i=1}^{N} B R_i^T (R_i A R_i^T)^{-1} R_i B^T. $$

To apply $P_S$ to some vector $g \in \mathbb{R}^m$, we need the whole $g$ since $B^T$ is a global matrix. Instead, we would like to perform computations locally on every core using only a portion of the vector. This can be done by taking advantage of the sparsity of $B^T$. Indeed, the restriction operations after $B^T$ suggest that the global vector is not needed. Only the part that is nonzero after the restriction is. Let us assume that there exists a family of restriction matrices $\tilde{R}_i : \mathbb{R}^m \to \mathbb{R}^{m_i}$, for $1 \le i \le N$ and with $m_i$ defined later, such that

$$ R_i B^T (I - \tilde{R}_i^T \tilde{R}_i) = 0, \quad 1 \le i \le N. \qquad (3.4) $$

If condition (3.4) is fulfilled, only the component projected by $\tilde{R}_i^T \tilde{R}_i$ will survive the application of $P_S$. In this way by defining some local operators $\tilde{B}_i := \tilde{R}_i B R_i^T \in \mathbb{R}^{m_i \times n_i}$, we can rewrite $P_S$ as

$$ P_S = C + B R_0^T (R_0 A R_0^T)^{-1} R_0 B^T + \sum_{i=1}^{N} \tilde{R}_i^T \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T \tilde{R}_i, $$
where we have now the sum of $N$ operators localized over some subsets of $\mathcal{M}$. As a result, we commuted the order of the operations represented by $B^T$ and $(R_i)_{1\le i\le N}$.

Let us denote the k-th component of vector $y \in \mathbb{R}^m$ with $[y]_k$. In [15] it is proposed to define $\tilde{R}_i$ as the restriction matrix from the global dofs set $\mathcal{M}$ to

$$ \mathcal{M}_i := \{k \in \mathcal{M} \mid [BR_i^T x_i]_k \neq 0, \quad \forall x_i \in \mathbb{R}^{n_i}\}, \quad 1 \le i \le N. \qquad (3.5) $$

With this definition, every $\mathcal{M}_i$ defines the smallest subset of dofs over which a vector in the $\mathcal{I}m(BR_i^T)$ is non zero. Which also means that the family $(\mathcal{M}_i)_{1\le i\le N}$ forms the smallest decomposition of $\mathcal{M}$ that fulfills condition (3.4). Let us denote $m_i := \#\mathcal{M}_i$. Moreover, let us state the following.

**Lemma 3.1.** Given the $(\mathcal{M}_i)_{1\le i\le N}$ defined in (3.5) and the associated restriction matrices $(\tilde{R}_i)_{1\le i\le N}$, we have

$$ \bigcup_{i=1}^{N} \mathcal{M}_i = \{k \in \mathcal{M} \mid [Bx]_k \neq 0, \quad \forall x \in \mathbb{R}^n\} \qquad (3.6) $$

*Proof.* The left-to-right inclusion is straightforward to prove, being $\mathcal{I}m(BR_i^T) \subseteq \mathcal{I}m(B)$.

In order to show the opposite inclusion, let us consider $y \in \mathcal{I}m(B)$, so that there exists $x \in \mathbb{R}^n$ such that $y = Bx$. Then using (3.4) we have

$$ y = Bx = \sum_{i=1}^{N} BR_i^T D_i R_i x = \sum_{i=1}^{N} \tilde{R}_i^T \tilde{R}_i BR_i^T D_i R_i x. $$

This means that we can define some local vectors $y_i := \tilde{R}_i BR_i^T D_i R_i x \in \mathbb{R}^{m_i}$ so that

$$ y = \sum_{i=1}^{N} \tilde{R}_i^T y_i. $$

If $[y]_k \neq 0$ for some $k$, then there exists at least one $i \in [1, N]$ for which $[\tilde{R}_i^T y_i]_k \neq 0$. Since $\tilde{R}_i^T y_i \in \mathcal{I}m(BR_i^T)$, it means that $k \in \mathcal{M}_i$ for at least one $i \in [1, N]$ which shows the left inclusion and thus the equality. $\square$

**Lemma 3.2.** If $B$ is full rank, there exists a partition of unity $(\tilde{D}_i)_{1\le i\le N}$ for $\mathbb{R}^m$ such that

$$ I_m = \sum_{i=1}^{N} \tilde{R}_i^T \tilde{D}_i \tilde{R}_i. \qquad (3.7) $$

*Proof.* Given the family of subsets $(\mathcal{M}_i)_{1\le i\le N}$ from (3.5), the associated restriction matrices $(\tilde{R}_i)_{1\le i\le N}$ and extension with zero matrices, we can define a partition of unity for $\mathbb{R}^m$ such that (3.7) is fulfilled if the following holds

$$ \bigcup_{i=1}^{N} \mathcal{M}_i = \mathcal{M}. $$

If $B$ is full-rank, $\mathcal{I}m(B) = \mathbb{R}^m$ so the left hand side of (3.6) is exacty $\mathcal{M}$. So this condition is a direct consequence of Lemma 3.1. $\square$
If $B$ is full rank, family $(M_i)_{1 \le i \le N}$ defines a decomposition of $M$. This means that a PoU like (3.7) can be defined. The case of $B$ is not full rank will be addressed in Section 3.2.

Let us assume that matrices $(\tilde{R}_i)_{1 \le i \le N}$ define a proper decomposition of $M$, meaning that a PoU like (3.7) can be defined. Moreover, let us make the following assumption.

**Assumption 3.3.** For $C \neq 0$, there exists a stable SPSD splitting. Namely there exist some matrices $(\tilde{C}_i)_{1 \le i \le N}$ and a fixed constant $k_C > 0$ such that

$$C \le \sum_{i=1}^{N} \tilde{R}_i^T \tilde{C}_i \tilde{R}_i \le k_C C.$$

If $C = 0$ this assumption is not needed. If $C$ comes from the finite element discretization of some SPSD bilinear form, then usually the splitting is easy to find and the idea to define such matrices is the same as for the Neumann matrices $A_i^{\text{Neu}}$ from Assumption 3.3. In this case $k_C = \tilde{k}_1$, where $\tilde{k}_1$ is the equivalent of (2.6) but for the decomposition defined by $M_i$, for $1 \le i \le N$.

Thanks to Assumption 3.3, $C$ can be replaced with the sum of local operators. This allows the definition of a new operator

$$M_S := BR_0^T (R_0 A R_0^T)^{-1} R_0 B^T + \sum_{i=1}^{N} \tilde{R}_i^T (\tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T) \tilde{R}_i, \quad (3.8)$$

which is spectrally equivalent to $S$ as well with

$$\frac{1}{\lambda_{M,AS}} \min\{1, \frac{1}{k_C}\} M_S \le S \le \frac{1}{\lambda_{m,AS}} M_S. \quad (3.9)$$

This can be showed combining (3.3) and Assumption 3.3.

Let us denote the matrices in brackets from (3.8) as

$$\tilde{S}_i := \tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T, \quad 1 \le i \le N.$$

They can be seen as local Schur complement matrices and a-priori they are only semi definite operators. Let us now use the following notations,

$$S_{1:N} := \sum_{i=1}^{N} \tilde{R}_i^T (\tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T) \tilde{R}_i = \sum_{i=1}^{N} \tilde{R}_i^T \tilde{S}_i \tilde{R}_i, \quad (3.10)$$

$$S_0 := BR_0^T (R_0 A R_0^T)^{-1} R_0 B^T = BR_0^T A_0^{-1} R_0 B^T,$$

so that

$$M_S = S_0 + S_{1:N}. \quad (3.11)$$

Operator $S_0$ is only semidefinite since it contains a restriction followed by an extension with zeros. Moreover it is a low-rank matrix, having the same rank as the coarse problem matrix $A_0 := R_0 A R_0^T$, which is $n_0 \ll n, m$. On the other hand, $S_{1:N}$ is the sum of the local Schur complement matrices and is always invertible if $S$ is.

**Remark 3.4.** To see why $S_{1:N}$ is invertible we can retrace backwards the argument that led to the definition of $S_{1:N}$. Indeed, using Assumption 3.3 we notice that $S_{1:N}$ is spectrally equivalent to

$$C + B \sum_{i=1}^{N} R_i^T (R_i A R_i^T)^{-1} R_i B^T = C + B M_{AS,1}^{-1} B^T,$$
that is what one obtaines when replacing $A^{-1}$ in the Schur complement definition with the one-level preconditioner from (2.1). Finally, if A is positive definite so will be $M_{AS,1}^{-1}$, therefore $S_{1:N}$ as well.

Let us assume for the moment to have a spectrally equivalent approximation of $S_{1:N}$, that we denote by $M_{S_{1:N}}$. Using this approximation, let us define one last spectrally equivalent operator for $S$, that is

$$N_S := S_0 + M_{S_{1:N}}. \qquad (3.12)$$

This operator was proposed in [15] as well. We want to use $N_S^{-1}$ as a preconditioner for $S$. We must provide also a way to apply $N_S^{-1}$ to some right hand side $r \in \mathbb{R}^m$. This will be discussed in Section 3.3.

## 3.2 Two-level preconditioner for the sum of local Schur complement matrices

Let us first define a robust and scalable preconditioner for $S_{1:N}$. We will need it in section 3.3. To this aim, it is possible to borrow the same underlying idea of the Neumann-Neumann (NN) preconditioner. It was originally defined to solve the substructuring problem on the interfaces's skeleton for nonoverlapping DD methods, see Chapter 6 in [6] and references therein. We aim to extend the same idea to find a preconditioner for the sum of general semidefinite positive operators, such as $S_{1:N}$. Assuming that all $\tilde{S}_i$ are definite we can define the one-level NN preconditioner for $S_{1:N}$ as

$$M_{S_{1:N},1}^{-1} := \sum_{i=1}^{N} \tilde{R}_i^T \tilde{D}_i \tilde{S}_i^{-1} \tilde{D}_i \tilde{R}_i. \qquad (3.13)$$

However, in general some $\tilde{S}_i$ may have a nontrivial kernel. This can be the case even for the original NN preconditioner. To take this into account, the Balancing Neumann-Neumann (BNN) variant is introduced, see [7] or Chapter 7.8.2 in [6]. The latter involves also the resolution of a problem on the product space of the kernels spaces. We aim to extend the idea of BNN to $S_{1:N}$ by considering the product space of the Ker($\tilde{S}_i$), which may be nonempty.

It is important to stress that the BNN preconditioner still represents a one-level method. In order to add a second level to it, a coarse space can be derived using the GenEO method, see Chapter 7.8.3 of [6]. To identify which modes are responsible for the slow convergence using the BNN preconditioner, some GEPs are solved on every subdomain $M_i$. By setting the BNN operator into the FSL framework, one notices that in this case the continuity constant may be very large. For A, conversely, it was the stable decomposition constant that could be very small. By defining a GenEO coarse space, Lemma 2.1 gives some useful spectral bounds.

The same analysis can be done for $S_{1:N}$. In order to understand which GEPs have to be solved we can look at preconditioner (3.13), since it leads to the same GEPs of the adapted-BNN version. In [15], operator (3.13) is set into the FSL framework. On every subdomain $M_i$, we must look for $(w_{ik}, \mu_{ik}) \in \mathbb{R}^{m_i} \times \mathbb{R}$ such that

$$\tilde{D}_i \tilde{R}_i S_{1:N} \tilde{R}_i^T \tilde{D}_i w_{ik} = \mu_{ik} \tilde{S}_i w_{ik}. \qquad (3.14)$$

Let us consider $\tau_{S_{1:N}} > 0$ a user-defined parameter and let us define the following spaces of local vectors

$$Z_i := \text{Span}\{w_{ik} \text{ solution of (3.14)} \mid \mu_{ik} > \tau_{S_{1:N}}\}, \quad 1 \le i \le N.$$

From this definition we can notice that the Ker($\tilde{S}_i$) $\subseteq W_i$, so the BNN-GenEO version will have the same formulation of the BNN preconditioner. With the same abuse of notation used in (2.5)
we can define the following global coarse space

$$ \tilde{W}_0 := \bigoplus_{i=1}^{N} \tilde{R}_i^T \tilde{D}_i W_i. $$

Let $Z_{S_{1:N}}$ be the matrix whose columns span the coarse space $\tilde{W}_0$. Moreover, let $\tilde{P}_{S_{1:N}}$ be the $S_{1:N}$-orthogonal projection from $\mathbb{R}^m$ on $\tilde{W}_0$ defined as

$$ \tilde{P}_{S_{1:N}} := Z_{S_{1:N}} (Z_{S_{1:N}}^T S_{1:N} Z_{S_{1:N}})^{-1} Z_{S_{1:N}}^T S_{1:N}. $$

Finally, the BNN-GenEO preconditioner for $S_{1:N}$ is denoted by $M_{S_{1:N},2}^{-1}$, or simply by $M_{S_{1:N}}^{-1}$, and it reads

$$ M_{S_{1:N}}^{-1} := Z_{S_{1:N}} (Z_{S_{1:N}}^T S_{1:N} Z_{S_{1:N}})^{-1} Z_{S_{1:N}}^T S_{1:N} \\ + (I - \tilde{P}_{S_{1:N}}) \sum_{i=1}^{N} \tilde{R}_i^T \tilde{D}_i \tilde{S}_i^\dagger \tilde{D}_i \tilde{R}_i (I - \tilde{P}_{S_{1:N}}^T), \quad (3.15) $$

where $\tilde{S}_i^\dagger$ is a pseudoinverse for $\tilde{S}_i$, for $1 \le i \le N$. Those are needed since a-priori operators $\tilde{S}_i$ from (3.10) may be singular. Adapting the proof from Chapter 7.8.3 of [6] that makes use of Lemma 2.1, it can be showed that

$$ (S_{1:N} y, y) \le (S_{1:N} M_{S_{1:N}}^{-1} S_{1:N} y, y) \le \alpha (S_{1:N} y, y), \quad \forall y \in \mathbb{R}^m, \quad (3.16) $$

where

$$ \alpha := \max\{1, \tilde{k}_0 \tau_{S_{1:N}}\}, $$

$$ \tilde{k}_0 := \max_{1 \le i \le N} \#\{1 \le j \le N \mid \tilde{D}_i \tilde{R}_i S_{1:N} \tilde{R}_j^T \tilde{D}_j \neq 0\}, \quad \tilde{k}_0 > 1. $$

The latter is the maximum number of interactions through the $S_{1:N}$-scalar product that one subdomain $M_i$ has with the others, plus one. Typically, $\tilde{k}_0$ represents the number of neighbors of $M_i$ plus the number of neighbors of the neighbors, plus one.

From (3.16) one can obtain

$$ (S_{1:N}^{-1} y, y) \le (M_{S_{1:N}}^{-1} y, y) \le \alpha (S_{1:N}^{-1} y, y), \quad \forall y \in \mathbb{R}^m, \quad (3.17) $$

and thanks to Lemma 2.4, the following result on the eigenvalues of $M_{S_{1:N}}^{-1} S_{1:N}$ holds

$$ 1 \le \lambda(M_{S_{1:N}}^{-1} S_{1:N}) \le \alpha. \quad (3.18) $$

**Remark 3.5.** It is worth to notice that in this case, parameter $\tau_{S_{1:N}}$, hence the GenEO method acts on the largest eigenvalues instead of the smallest as for A. This is related to the fact that in the BNN case is the continuity constant from the FSL to be problematic. However, as for A, a smaller value of $\tau_{S_{1:N}}$ leads on one side to stricter bound on the largest eigenvalue but a larger coarse space on the other side.

**Constraint matrix not full rank.** The case of B not full rank is frequenlty encountered in saddle point problems. One example is the Stokes problem discretized using mixed finite elements with stabilization terms, e.g. the stabilized P1/P1 elements. When B is not full rank the image of B has not full dimension, so the left hand side of (3.6) is strictly included in $M$.
This means that there exist some $y \in \mathbb{R}^m$ that can not be represented as sums of extended local vectors living inside subsets $M_i$. This might seem a problem to definition of a PoU on the dual decomposition. However, there is a simple way to circumvent the problem without changing the method.

Let us consider any family of subsets $(M_i)_{1 \le i \le N}$ and the associated family of restriction matrices $(\tilde{R}_i)_{1 \le i \le N}$ for which (3.4) holds true. If $B$ is not full rank, $B^T$ has a nontrivial kernel. So there exist some $y \in \mathbb{R}^m$ for which

$$R_i B^T y = 0, \quad \forall i \in [1, N],$$

which thanks to (3.4) implies that

$$R_i B^T \tilde{R}_i \tilde{R}_i y = \tilde{B}_i^T \tilde{R}_i y = 0, \quad \forall i \in [1, N]. \qquad (3.19)$$

From (3.19) one notices that if $B$ is not full rank, matrices $\tilde{B}_i^T$ has a larger kernel, for some $i \in [1, N]$. Recalling that

$$\tilde{S}_i = \tilde{C}_i + \tilde{B}_i (R_i A R_i^T)^{-1} \tilde{B}_i^T, \quad 1 \le i \le N,$$

it may be that $\tilde{S}_i$ has a larger kernel as well. As explained before in this same section, the Ker($\tilde{S}_i$) is automatically included inside $\tilde{W}_0$. To conclude, if subsets $(M_i)_{1 \le i \le N}$ and the associated restriction matrices $(\tilde{R}_i)_{1 \le i \le N}$ are such that (3.4) holds true, the same algorithm can be used for general $B$ since the size of $\tilde{W}_0$ will automatically adapt by including the additional elements in the kernels.

## 3.3 Preconditioner for the Schur complement matrix

Let us recall that we aim to use $N_S = S_0 + M_{S_{1:N}}$, defined in (3.12), as a preconditioner for $S$. We will first state the result on the spectrum of $N_S^{-1}S$, then we will provide some ways to apply $N_S^{-1}$, in order to use it as a preconditioner. The one outlined in Section 3.3.2 is the one proposed in [15]. Instead, the one presented in Section 3.3.3 is a new way introduced in this work to extend the method to general problems.

### 3.3.1 Spectral bounds for the Schur complement preconditioned system

Let us state a result on the spectrum of the preconditioned system $N_S^{-1}S$. The spectral equivalence between $N_S$ and $S$ is the statement of the following lemma.

**Lemma 3.6.** Operator $N_S$ is spectrally equivalent to $S$ with constants

$$\lambda_{m,N_S} := \frac{1}{2k_0} \min\{1, \frac{1}{k_C}\}, \quad \lambda_{M,N_S} := (2 + (2k_0 + 1)k_1\tau_A) \max\{1, \tilde{k}_0\tau_{S_{1:N}}\},$$

so that every eigenvalue of $N_S^{-1}S$ is real, positive and such that

$$\lambda_{m,N_S} \le \lambda(N_S^{-1}S) \le \lambda_{M,N_S}. \qquad (3.20)$$

*Proof.* Operator $S_{1:N}$ is definite, hence $M_{S_{1:N}}^{-1}$ as well thanks to (3.17). Let us call $M_{S_{1:N}}$ the inverse of $M_{S_{1:N}}^{-1}$. First of all, we must show the spectral equivalence between $M_{S_{1:N}}$ and $S_{1:N}$. It can be proved by using Lemma 2.1 with

$$H = H_D = \mathbb{R}^m, \quad a(y, w) := (S_{1:N}^{-1} y, w), \quad b(y, w) := (M_{S_{1:N}}^{-1} y, w), \quad \mathcal{R}(y) = y.$$
In this context, $\mathcal{R}$ is surjective being the identity, while using (3.17) it can be proved that continuity holds with constant $c_R = 1$ and stable decomposition holds with constant $c_T = \alpha^{-1}$. The result given by the FSL is

$$ \frac{1}{\alpha}(S_{1:N}^{-1}y, y) \leq (S_{1:N}^{-1}M_{S_{1:N}}S_{1:N}^{-1}y, y) \leq (S_{1:N}^{-1}y, y), \quad \forall y \in \mathbb{R}^m. $$

By taking $y = S_{1:N}w$, after renaming all the $w$, one gets

$$ \frac{1}{\alpha}(S_{1:N}y, y) \leq (M_{S_{1:N}}y, y) \leq (S_{1:N}y, y), \quad \forall y \in \mathbb{R}^m. \quad (3.21) $$

Combining (3.21) and (3.9) one gets the spectral equivalence between $N_S$ and $S$ that reads

$$ \lambda_{m,N_S} N_S \leq S \leq \lambda_{M,N_S} N_S. $$

Then (3.20) is a consequence of Lemma 2.4. $\square$

Constants $\lambda_{m,N_S}$, $\lambda_{M,N_S} N_S$ only depend on parameters that are fixed with respect to the discretization, problem coefficients or global number of subdomains. This means that $N_S^{-1}$ is a robust and scalable preconditioner for $S$.

**Remark 3.7.** In agreement with the previous sections, smaller values for both parameters $\tau_A$ and $\tau_{S_{1:N}}$ leads to a sharper bound on the largest eigenvalue, alongside with larger coarse spaces $V_0$ and $W_0$.

### 3.3.2 Application of the preconditioner with nested iterative solution

The result from Lemma 3.6 is of no use unless a way to apply $N_S^{-1}$ to some right hand side $r \in \mathbb{R}^m$ is provided as well. In [15] it is proposed to solve a problem like

$$ N_S y = r, \quad (3.22) $$

every time $N_S^{-1}r$ has to be performed, for some right hand side $r \in \mathbb{R}^m$. A direct method can not be used since we do not possess $M_{S_{1:N}}$ explicitly. On the other hand, solving (3.22) with an iterative method would suffer of slow convergence due to the same poor conditioning affecting $S$, being $N_S$ spectrally equivalent to $S$. So preconditioning must be introduced.

In [15], it is proposed to use $M_{S_{1:N}}^{-1}$ as a preconditioner every time a problem like (3.22) has to be solved iteratively. Since we do not know $N_S$ explicitly, problem (3.22) is replaced in practice by problem

$$ M_{S_{1:N}}^{-1} N_S y = (M_{S_{1:N}}^{-1} S_0 + I_m)y = M_{S_{1:N}}^{-1} r, \quad (3.23) $$

solved with a nonpreconditioned Krylov method. However, this approach is not totally satisfactory for two main reasons:

* it requires an inner solution step at every iteration for $N_S^{-1}S$, increasing the result time,

* the number of iterations needed to solve (3.23) is not under theoretical control. Namely, there are no results to guarantee the theoretical scalability and robustness of $M_{S_{1:N}}^{-1}$ as a preconditioner for $N_S$.

However, this work remains of interest because it brings experimental evidence of controlled number of iterations for (3.23) in the case of heterogeneous linear elasticity, see [15] or Section 5.3 where the same test case is considered.

In the following section we will provide a new way to apply $N_S^{-1}$ that does not imply nested iterative solutions. It allows to directly exploit the result from Lemma 3.6 at the expense of building and factorizing a dense matrix.
### 3.3.3 New way to apply the preconditioner using an algebraic formula

In order to have a method that is application-independent, we must rely on some purely algebraic trick in order to apply $N_S^{-1}$. In this way we do not need more specific informations about the differential operators underlying matrices. The tool used is the Sherman-Morrison-Woodbury formula, called Woodbury formula in the following. It is an exact algebraic result particularly useful to compute the inverse a matrix $M$ plus a low-rank correction. Let us state the result, followed by one way to prove it.

**Proposition 3.8** (Woodbury formula, [25, 12]). Given $M \in \mathbb{R}^{n \times n}$, $K \in \mathbb{R}^{k \times k}$, with $k \le n$, invertible matrices, $U \in \mathbb{R}^{n \times k}$, $V \in \mathbb{R}^{k \times n}$ we have

$$ (M + UKV)^{-1} = M^{-1} - M^{-1}U(K^{-1} + VM^{-1}U)^{-1}VM^{-1}. \quad (3.24) $$

*Proof.* Matrix $(M + UKV)^{-1}$ can be seen as the solution of the matrix equation

$$ (M + UKV)X = I_n. \quad (3.25) $$

Now we can introduce the matrix variable $Y := KVX$. In this way (3.25) can be rewritten in a mixed formulation as

$$ \begin{bmatrix} M & U \\ V & -K^{-1} \end{bmatrix} \begin{bmatrix} X \\ Y \end{bmatrix} = \begin{bmatrix} I_n \\ 0_m \end{bmatrix}, \quad (3.26) $$

from which (3.25) can be obtained with an elimination of $Y$. On the other hand, if we perform an elimination of $X$ instead of $Y$ in (3.26), from the first equation we get $X = M^{-1}(I_n - UY)$. If we replace this expression in the second line we get

$$ VM^{-1}(I_n - UY) + K^{-1}Y = 0_m. $$

So now $Y$ is given by $Y = (K^{-1} + VM^{-1}U)^{-1}VM^{-1}$. If we replace this back into the expression for $X$ we get

$$ X = M^{-1}(I_n - U(K^{-1} + VM^{-1}U)^{-1}VM^{-1}), $$

which, combined with (3.25), proves exactly the Woodbury formula (3.24). $\square$

Let us recall that $N_S$ from (3.12) is the sum of a matrix $S_0$ of rank $n_0 \ll n$ and a full rank operator $M_{S_{1:N}}$, whose inverse is defined in (3.15). Let us also recall that $A_0 = R_0AR_0^T$. If we replace $M^{-1}$, $K$, $U$ and $V$ in (3.24) with $M_{S_{1:N}}^{-1}$, $A_0^{-1}$, $BR_0^T$ and $R_0B^T$ respectively, we can use the formula to obtain an explicit expression for the inverse of $N_S$. Indeed, one gets

$$ N_S^{-1} = M_{S_{1:N}}^{-1} - M_{S_{1:N}}^{-1}BR_0^T(A_0 + R_0B^T M_{S_{1:N}}^{-1}BR_0^T)^{-1}R_0B^T M_{S_{1:N}}^{-1}. \quad (3.27) $$

This allows us to apply $N_S^{-1}$ without any inner solution step and moreover, the result from Lemma 3.6 becomes meaningful in practice. Indeed, the iteration count for the preconditioned iterative solution for $N_S^{-1}S$ is theoretically controlled thanks the result in Lemma 3.6. This comes at the expense of applying the right hand side of (3.27), which will be addressed in the next section.

### 3.3.4 Practical aspects for the preconditioner application

Some considerations about the use of (3.27) in practice must be made. To perform a matrix-vector product between the right hand side of (3.27) and a vector $r \in \mathbb{R}^m$, we have to first perform the product $M_{S_{1:N}}^{-1}r$, which is almost fully parallelizable. The coarse problem contained
in $M_{S_{1:N}}^{-1}$ is solved using a direct method. Then, the product with $R_0 B^T$ can be done fully in parallel. The critical point is the solution of the problem represented by

$$E_0 := A_0 + R_0 B^T M_{S_{1:N}}^{-1} B R_0^T. \qquad (3.28)$$

This matrix has rank $n_0 \ll n$, but it is dense due to the solution of the coarse problem on $\tilde{W}_0$, contained in $M_{S_{1:N}}^{-1}$. In order to apply $E_0^{-1}$, one must build $E_0$ and factorize it as a dense matrix. The time needed for both the building phase and the factorization are reported among the numerical results.

The building phase can be optimized by exploiting the particular structures of the matrices involved in the computation. The actual implementation is not optimized. It consists in applying the operator $B^T M_{S_{1:N}}^{-1} B$ to every column of $R_0^T$ sequentially. Instead, the factorization time depends only on the size of $E_0$, which depends on the problem at hand. At scale, this step can become the bottleneck point of the presented method, preventing us to use it for problems with very large coarse spaces. Indeed, if the coarse space $V_0$ is large, the factorization of (3.28) becomes prohibitive in terms of memory and time.

If coarse space $V_0$ is not large the presented method is of practical use, the factorization of $E_0$ being affordable. This is also the regime where classical GenEO method works best, even though GenEO only requires the factorization of sparse matrices, so in principle can be used for larger problems.

## 4 Saddle point preconditioning

In this section we will combine the DD preconditioners derived in the previous sections with block preconditioning. The latter is one of the two main paradigms to solve large saddle point problems iteratively. The other one are segregated iteration methods (or Schur complement approaches). As briefly explained in the introduction, both paradigms need some approximations of operators $A$ and $S$. In this work we presented an approximation of the Schur complement $S$ (which serves as a preconditioner as well) using the DD paradigm. In order to study the individual performance of $N_S^{-1}$ as a preconditioner for $S$ we will use Algorithm 1. It allows us to divide the solution and solve a problem for the Schur complement individually, at Step 3.

**Remark 4.1.** Algorithm 1 is based on factorization (1.2). It coincides with the elimination of the primal unknown $x$, which is then replaced in the second equation. A problem for $S$ appears and it can be solved for the dual unknown $y$. The algorithm ends with an update of the primal variable to find the final value of $x$.

<caption>**Algorithm 1** Saddle point solver with segregated iterations</caption>

**Input:** $[f \ g] \in \mathbb{R}^{n+m}$. **Output:** $[x \ y]$ solution to (1.1).
1. Solve $A\tilde{x} = f$ with a Krylov method using $M_{HSM}^{-1}$ from (2.2) as a preconditioner.
2. Compute $b_g := g - B\tilde{x}$.
3. Solve $Sy = -b_g$ with a Krylov method with $N_S^{-1}$ from (3.27) as a preconditioner.
4. Compute $b_f := f - B^T y$.
5. Solve $Ax = b_f$ with a Krylov method using $M_{HSM}^{-1}$ as a preconditioner.
Block preconditioning means to define an approximation of A that has a similar block structure. Typically, this is done by replacing A and S with some approximations (or preconditioners). Then, the original problem is solved all-at-once using some preconditioned iterative method e.g. a Krylov methods. With respect to Algorithm 1, the choice of a block preconditioner brings two main advantages.

* Convergence is assessed by looking at the global relative residual, instead of at the residuals on the single variable independently.

* The product with S, that is instead needed in step 3 of Algorithm 1 can be avoided. This means to avoid solving a problem for matrix A at every application of S.

We are interested in using the DD approximations that we derived in the previous sections in this context. The goal is to obtain a robust and scalable global preconditioner for A.

In the literature, there exist several choices for block preconditioners. The most classical ones are block diagonal preconditioners, block triangular preconditioners, Uzawa preconditioners and block approximate factorization preconditioners. References for all these choices are reported in the classical reference [4] or more recently in [18], alongside a careful and unified spectral analysis. Among the above choices, let us consider the block approximate factorization preconditioners, for which one representative is

$$ \mathcal{P} = \begin{bmatrix} M_{HSM} & \\ B & -N_S \end{bmatrix} \begin{bmatrix} I & M_{HSM}^{-1}B^T \\ & I \end{bmatrix}, \qquad (4.1) $$

which is obtained replacing A and S in (1.2) with $M_{HSM}$ from (2.2) and $N_S$ from (3.27) respectively. The application of $\mathcal{P}$ amounts to a block backward-forward substitution process where only the applications of $M_{HSM}^{-1}$ and $N_S^{-1}$ have to be performed.

In [18], given two generic approximations for A and S and the extremal spectral values of the corresponding preconditioned matrices, a region of the complex plane where the eigenvalues of $\mathcal{P}^{-1}A$ are included is outlined. In the same work, this analysis is carried out for all the block preconditioner choices listed above. We aim to use the result from Theorem 4.3 and Table 1 in [18] for $\mathcal{P}$ defined in (4.1). First of all, let us define the extremal spectral values for the preconditioned matrices borrowing the notation from [18], namely

$$ \underline{\mu} := \lambda_{\min}(M_{HSM}^{-1}A), \quad \bar{\mu} := \lambda_{\max}(M_{HSM}^{-1}A), $$

$$ \underline{\nu} := \lambda_{\min}(N_S^{-1}S), \quad \bar{\nu} := \lambda_{\max}(N_S^{-1}S). $$

Let us assume that $\bar{\mu} < 2$, which is equivalent to assume that $\rho_A < 1$, with $\rho_A = \max\{1-\underline{\mu}, \bar{\mu}-1\}$ spectral radius of $I-M_{HSM}^{-1}A$. This assumption is made in [18] to carry out the proofs for spectral bounds. We can always ensure that it is satisfied by operating a rescaling of the preconditioner for A. Indeed, we know from (2.9) that $\bar{\mu} \le k_0$, so rescaled preconditioner $2k_0^{-1}M_{HSM}^{-1}$ fulfils the assumption. Let us define

$$ \underline{\xi} := \min\{\underline{\mu}, \underline{\nu}\}, \quad \bar{\xi} := \max\{\bar{\mu}, \bar{\nu}\}, $$

$$ \underline{\chi} := \frac{1}{2}(\underline{\mu} + \underline{\nu}(1 - \rho_A^2)), \quad \bar{\chi} := \frac{1}{2}(\bar{\mu} + \bar{\nu}), $$

$$ \delta^2 := \begin{cases} \bar{\nu} \max\{\frac{4}{27}, \bar{\mu}(1-\bar{\mu})^2\}, & \text{if } \underline{\mu} \le \frac{1}{3} \le \bar{\mu}, \\ \bar{\nu} \max\{\underline{\mu}(1-\underline{\mu})^2, \bar{\mu}(1-\bar{\mu})^2\}, & \text{otherwise,} \end{cases} $$

$$ \zeta := \frac{\bar{\nu}}{1 + \bar{\nu}(2 - \bar{\mu})}. $$
The preconditioned operator $\mathcal{P}^{-1}\mathcal{A}$ is nonsymmetric. Moreover, $\mathcal{A}$ and $\mathcal{P}$ are indefinite, so a-priori the preconditioned matrix could have both real and complex eigenvalues. From [18], we have that the real eigenvalues $\lambda_r$ of $\mathcal{P}^{-1}\mathcal{A}$ are such that

$$ \xi \le \lambda_r \le \bar{\xi}, \qquad (4.2) $$

whereas complex eigenvalues $\lambda_c$ are such that

$$ \underline{\chi} \le \Re(\lambda_c) \le \bar{\chi}, \qquad (4.3) $$

$$ |\mathfrak{F}(\lambda_c)| \le \delta, \qquad (4.4) $$

and

$$ |\lambda_c - \zeta| \le \zeta. \qquad (4.5) $$

Bounds (4.2) - (4.5), identify a region of the complex plane where the eigenvalues of $\mathcal{P}^{-1}\mathcal{A}$ are included. In particular, the extension of such region depends on the extremal spectral values of $M_{HSM}^{-1}\mathcal{A}$ and $N_S^{-1}\mathcal{S}$. In particular, the real eigenvalues lies on the real axis inside the interval $[\xi, \bar{\xi}]$, while the complex ones are included inside the square $[\chi, \bar{\chi}] \times [-\delta, \delta]$ and the circle centered in $(\zeta, 0)$ with radius $\zeta$. The eigenvalue distribution of $\mathcal{P}^{-1}\mathcal{A}$ for the heterogenous linear elasticity problem presented in Section 5.1 is reported in Figure 1, alongside with bounds (4.2) - (4.5).

<figure><img src="image_2.png" /><figcaption>Figure 1: Eigenvalue distribution of $\mathcal{P}^{-1}\mathcal{A}$ for the 3D heterogeneous linear elasticity problem from Section 5.1. Eigenvalues and bounds (4.2) - (4.5) in the complex plane are depicted.</figcaption></figure>

From the DD theory outlined in Sections 2 and 3, parameters $\bar{\mu}, \bar{\mu}, \bar{\nu}, \bar{\nu}$ can be bounded by quantities that are fixed with respect to the number of subdomains or problem parameters. This means that the spreading of the eigenvalue cloud of $\mathcal{P}^{-1}\mathcal{A}$ is controlled. Even if it is well-known that there is no direct connection between spectrum and GMRES convergence speed, see [9], it is preferred to have a more clustered eigenvalue cloud, since it typically brings a lower number of iterations. So in this context, we can say that $\mathcal{P}^{-1}$ is a robust and scalable preconditioner for $\mathcal{A}$, even though the bounds could be not very sharp.

**Remark 4.2.** The analysis performed in [18] shows also how the particular choice of $\mathcal{P}$ affects how much eigenvalues are spread on the complex plane. In particular, block preconditioner $\mathcal{P}$ defined in (4.1) brings a less spread region for the eigenvalues of $\mathcal{P}^{-1}\mathcal{A}$ with respect to other choices e.g. block diagonal or block triangular preconditioners. The latter are easier to apply, but more spread eigenvalues cloud typically results in a higher number of iterations.
**Remark 4.3.** If a rescaling on the preconditioner for A is introduced, it may degenerate the iteration count. Indeed, rescaling the preconditioner means to move the whole eigenvalue cloud. However, since we are interested in scalability, being in a framework where theory applies and some bounds are available is more appealing than looking for the configuration that brings the lowest iteration count.

# 5 Numerical experiments

In this section we aim to show numerically the scalability of preconditioners derived in the previous sections in terms of number of iterations for the iterative solver to convergence. We will report results from weak scalability experiments. Namely, at every test we increase the number of subdomains and the size of the global problem, so that the size of local problems remains almost constant.

Alongside the iterations count we will report also the computational times. Both the setup times to build the entire DD framework and the preconditioners and the solution times will be reported to compare the method presented in this work and the one from [15]. In order to have a competitive method, both setup and solution time must be taken into consideration.

We will test preconditioners on a difficult benchmark case from the class of linear saddle point problems discretized with the finite elements method (FEM). First of all, we will present the discrete formulations of the problem. Then, we will provide some details regarding the implementation and how the DD framework applies to FEM. Finally, complete numerical results will be presented. In particular, a comparison with the original method from [15] is performed.

## 5.1 Heterogeneous linear elasticity problem

The benchmark problem is the 3D elasticity problem for nearly-incompressible materials. In particular we are interested in a beam with square section made of alternating layers of rubber and steel, see Figure 2. Due to the big variation between the mechanical properties of these two materials, this is a typical example of challenging problem, for which the development of a scalable preconditioner is difficult.

The mechanical properties of a solid are characterized by its Young modulus $E$ and its Poisson ratio $\nu$, or alternatively by its Lamé coefficients defined as

$$ \lambda := \frac{E\nu}{(1 + \nu)(1 - 2\nu)} \quad \text{and} \quad \mu := \frac{E}{2(1 + \nu)}. $$

The elasticity problem reads a variational formulation with saddle point structure where $u$ represents the displacement and $p$ the pressure field on the solid. The first equation represents the balance of forces acting on the solid, while the second equation represents the nearly-incompressibility constraint. An Homogeneous Dirichlet boundary condition is imposed on the end face of the beam in order to model the binding to a fixed wall. Let us denote by $\Gamma_D$ this part of the boundary. A stable conforming discretization with continuous P2/P1 Taylor-Hood elements is used, see [5]. Let us denote by $V_h := \mathbb{P}_2^3 \cap (H_{\Gamma_D}^1)^3$, where $H_{\Gamma_D}^1$ is the space of $H_1$ functions with null trace on $\Gamma_D$, and $W_h := \mathbb{P}_1$. The problem reads:

Find $(u_h, p_h) \in V_h \times W_h$ such that

$$ \begin{aligned} \int_{\Omega} 2\mu \underline{\epsilon}(u_h) : \underline{\epsilon}(v_h) - \int_{\Omega} p_h \mathrm{div}(v_h) &= \int_{\Omega} f \cdot v_h, & \forall v_h \in V_h, \\ - \int_{\Omega} q_h \mathrm{div}(u_h) - \int_{\Omega} \frac{1}{\lambda} p_h q_h &= 0, & \forall q_h \in W_h, \end{aligned} $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Rectangular beam composed by 10 alternating layer of rubber and steel. On the left: heterogeneous coefficient distribution. On the right: automatic partition of the beam into 16 subdomains realized with *Metis*.</figcaption></figure>

where $\underline{\underline{\xi}}$ denotes the linearized strain tensor and $f$ is the weight force vector. The algebraic formulation of this variational problem has a saddle point structure and reads

$$ \begin{bmatrix} A & B^T \\ B & -C \end{bmatrix} \begin{bmatrix} u \\ p \end{bmatrix} = \begin{bmatrix} f \\ 0 \end{bmatrix}. \qquad (5.1) $$

Coercivity of the strain tensor integral term over $V_h$ implies that matrix $A$ is SPD, while the inf-sup condition for the divergence term over $V_h, W_h$ implies that matrix $B^T$ has an empty kernel, namely $B$ is full rank.

We will consider 10 alternating layers of steel and rubber for which material properties are given by $(E_1, \nu_1) = (2 \times 10^9, 0.35)$ and $(E_2, \nu_2) = (10^7, 0.4999)$ respectively. The subdivision in layers is performed along one axis and we consider an automatic domain decomposition that does not resolve the discontinuities in the coefficients distribution, see Figure 2.

## 5.2 Hardware, softwares and implementation details

The problem is discretized and solved using the parallel finite element software *FreeFEM* [11]. In particular, the algorithm developed in this work is implemented on top of the *ffddm* framework already available in *FreeFEM*. Indeed, *ffddm* is a set of parallel *FreeFEM* scripts that allow the user to implement Schwarz DD methods, see [24] for the documentation. In particular, the GenEO method is available in *ffddm* as a possible choice for the coarse space for SPD problems. To obtain the results from [15], the framework was extended to the saddle point linear systems case. As usual in domain decomposition, every subdomain is assigned to a MPI process. No multithreading is used.

Numerical results are obtained on the French GENCI supercomputer Joliot-Curie, on the AMD Rome partition composed by 2286 nodes of AMD Epyc processors clocked at 2.6 Ghz. Every node features a total of 128 cores divided into two AMD CPUs of 64 cores each. Every cores has a RAM memory availability of 1.8Gb. We assigned multiple cores to every MPI process in order to benefit of more memory. This is done to match the local problems size of the test performed in [15] with the memory availability of the machine used.

Using *FreeFEM*, a mesh $T_h$ for the computational domain $\Omega$ is defined. It consists in a triangulation of $\Omega$ made by simplicial triangular elements. Subdomains are obtained first by defining $N$ disjoint subsets $\Omega'_i$ of elements, for $1 \le i \le N$, using the automatic graph partitioner *Metis* [13]. Let us stress that the partition is generic, namely it does not resolve the discontinuities in the coefficients as one can see in Figure 2. To create the overlap regions, one layer of elements
is added to every $\Omega'_i$. Let us denote the extended subdomains by $\Omega_i$, for $1 \le i \le N$. The family of $(\Omega_i)_{1 \le i \le N}$ defines a FEM overlapping decomposition of $\Omega$. Let us consider $V_h$ discrete space for the primal variable and let us denote by $\mathcal{N}$ the global set of dofs associated to $V_h$, with $n := \#\mathcal{N}$. Let us denote by $\{\varphi_j\}_{1 \le j \le n}$ the basis functions for $V_h$ and let us define

$$ \mathcal{N}_i := \{k \in \mathcal{N} \mid \mathrm{meas}(\mathrm{supp}(\varphi_k) \cap \Omega_i) > 0\}, \quad 1 \le i \le N. $$

Subset $\mathcal{N}_i$ contains the degrees of freedom inside $\Omega_i$ and on $\partial\Omega_i$. Let us define $n_i := \#\mathcal{N}_i$.

The definitions of $(\mathcal{M}_i)_{1 \le i \le N}$ is problem specific since it depends on $B$. But as pointed out in Section 3.2, it is enough to ensure that the dual decomposition fulfills (3.4). First of all, one additional layer of elements is added to every $\Omega_i$. The resulting subdomains are denoted with $\tilde{\Omega}_i$. Finally, considering $\{\psi_j\}_{1 \le j \le m}$ basis functions for $W_h$, let us define

$$ \mathcal{M}_i := \{k \in \mathcal{M} \mid \mathrm{meas}(\mathrm{supp}(\psi_k) \cap \tilde{\Omega}_i) > 0\}, \quad 1 \le i \le N. $$

Let us define $m_i := \#\mathcal{M}_i$. Restriction matrices $\tilde{R}_i$ are consequently defined. The same definition would work even if matrix $B$ is not full rank. The application of $B^T$ represents the derivative of piecewise linear functions over $\Omega \setminus \tilde{\Omega}_i$, which is identically zero when restricted to $\Omega_i$ since $\Omega_i \subset \tilde{\Omega}_i$. Hence, condition (3.4) is fulfilled.

Matrices $A'_i \in \mathbb{R}^{n_i \times n_i}$ from Assumption (2.2) and $\tilde{C}_i \in \mathbb{R}^{m_i \times m_i}$ from Assumption (3.3) are obtained by evaluating the bilinear forms underlying $A$ and $C$ on $\Omega_i$ and $\tilde{\Omega}_i$ respectively, without any additional boundary conditions on the subdomains boundaries. In the FEM context, these matrices are also called Neumann matrices.

Generalized eigenvalue problems are solved locally with Arpack [14] using the shift-and-invert transformation. This requires the factorization of the Neumann matrices $A'_i$ and the local Schur complement matrices $\tilde{S}_i$ for problems (2.4) and (3.14) respectively. In both cases, the sparse solver MUMPS [2] is used. While matrices $A'_i$ are sparse, $\tilde{S}_i$ are dense matrices so their factorization is performed by considering the augmented sparse saddle point matrices

$$ \tilde{S}_{i, \mathrm{aug}} := - \begin{bmatrix} R_i A R_i^T & \tilde{B}_i^T \\ \tilde{B}_i & -\tilde{C}_i \end{bmatrix} \in \mathbb{R}^{(n_i+m_i) \times (n_i+m_i)}. \quad (5.2) $$

Indeed, solving the linear system $\tilde{S}_i y_i = g_i$ is equivalent to solve

$$ - \begin{bmatrix} R_i A R_i^T & \tilde{B}_i^T \\ \tilde{B}_i & -\tilde{C}_i \end{bmatrix} \begin{bmatrix} x_i \\ y_i \end{bmatrix} = \begin{bmatrix} 0 \\ g_i \end{bmatrix}, $$

which can be efficiently done with MUMPS, being the matrix sparse.

**Remark 5.1.** While GEPs (2.4) can be solved independently on every subdomain, GEPs (3.14) need neighbour-to-neighbour communications. Indeed, the left hand side needs the operators $\tilde{S}_i$ from the neighbouring subdomains in order to perform the matrix vector product.

When applying preconditioners $M_{HSM}^{-1}$ and $M_{S_{1:N}}^{-1}$, both local and coarse problems are solved with MUMPS, being all sparse. In particular, local problems are factorized and solved by every MPI process independently, while coarse problems represented by $R_0 A R_0^T$ and $Z_{S_{1:N}}^T S_{1:N} Z_{S_{1:N}}$ are factorized and solved by one rank. If these problems become too large to be solved only by one rank, ffddm offers the possibility to use a few cores so that coarse problems can be solved with a direct method in a distributed way.

The chosen parameters for the GenEO method are $(\tau_A, \tau_{S_{1:N}}) = (10, 3.33)$, as in [15]. There is no rule to choose these parameters a-priori. One can select the combination that guarantees a
reasonable number of iterations for a small case. Indeed, scalability ensures that the iterations count is constant while increasing the number of subdomains, hence the global problem size. One always has to consider that a larger parameter typically means less iterations but larger coarse spaces.

In order to solve linear systems, GMRES with right preconditioning is used. The tolerance for convergence is set at tol = 10⁻⁵. Using right-preconditioning, convergence is assessed by looking at the real relative residuals and not at the preconditioned ones. No restarting is considered. GMRES is used to solve problems at Steps 1, 3 and 5 in Algorithm 1 and also to solve Au = b using P⁻¹ from (4.1) as a preconditioner.

**Remark 5.2.** All the bounds on the spectra from Sections 2 and 3 are given for the left precondi-
tioned matrices. From the spectral equivalencies one can simply derive the same spectral bounds
for the right preconditioned matrices.

**5.3 Weak scalability experiment for heterogeneous beam**

In this section we present a weak scalability experiment realized for the heterogeneous linear
elasticity problem on a beam, introduced in Section 5.1. During the experiment the mesh re-
finement is kept fixed, while the beam lenght is doubled at every test, keeping the same number
of alternating layers. In this way, doubling the number of subdomains (and of MPI processes)
accordingly, the local size of subproblems is kept almost constant. We will compare iteration
counts and timings for the same experiment using both methods to apply $N_s^{-1}$, from Section
3.3.2 and Section 3.3.3.

**5.3.1 Experiment using the method based on the Woodbury formula**

Results on the number of iterations are reported in Table 1. The first column reports the number of subdomains (and of MPI processes) used for every test. The second column reports the total number of degrees of freedom, namely the size of the original problem (1.1). Then we report the coarse space sizes for $V_0$ and $\tilde{W}_0$ respectively. Finally, the number of GMRES iterations for problems at Step 1 and Step 3 of Algorithm 1 and for the original problem with block preconditioner $\mathcal{P}$ from (4.1).

We observe a large but constant number of iterations for *A* throughout the experiment. By considering a smaller value for *τ*<sub>*A*</sub> this iteration count can be decreased and made even more constant, at the expense of having larger *V*<sub>0</sub>. On the other hand, the number of iterations observed for *N*<sub>*S*</sub><sup>-1</sup>*S* is very constant. This confirms the theoretical scalability given by Lemma 3.6. For what concerns the use of a preconditioned iterative method on the global problem *A*, we observe very large iteration counts. In this case, the rescaling of the preconditioner for *A* has been performed in order to be inside the theoretical framework of [18]. Such procedure may bring worse performances in some cases. But on the other hand, it guarantees that the eigenvalues are included in a region whose extension does not depend on the subdomain numbers. Indeed, even in this case we observe that the iteration count does not vary too much.
<table><caption>Table 1: Weak scalability experiment for the heterogeneous linear elasticity problem using the method based on the Woodbury Formula. Iteration counts.</caption><thead><tr><th>N (or P)</th><th>n+m</th><th>dimV<sub>0</sub></th><th>dimW&#770;<sub>0</sub></th><th>It. A</th><th>It. S</th><th>It. A</th></tr></thead><tbody><tr><td>16</td><td>1,058,188</td><td>726</td><td>79</td><td>46</td><td>17</td><td>84</td></tr><tr><td>32</td><td>2,110,941</td><td>1473</td><td>285</td><td>48</td><td>17</td><td>89</td></tr><tr><td>64</td><td>4,216,392</td><td>3088</td><td>914</td><td>49</td><td>18</td><td>90</td></tr><tr><td>128</td><td>8,427,373</td><td>6185</td><td>1726</td><td>52</td><td>18</td><td>96</td></tr><tr><td>256</td><td>16,849,298</td><td>12551</td><td>3476</td><td>58</td><td>19</td><td>106</td></tr></tbody></table>

<table><caption>Table 2: Weak scalability experiment for the heterogeneous linear elasticity problem using the method based on the Woodbury Formula. Timings.</caption><thead><tr><th>Build. E<sub>0</sub> (s)</th><th>Fact. E<sub>0</sub> (s)</th><th>WF (s)</th><th>Total setup (s)</th><th>GMRES S (s)</th><th>GMRES A (s)</th></tr></thead><tbody><tr><td>585.94</td><td>0.0735</td><td>27.25</td><td>944.6</td><td>167.4</td><td>335.2</td></tr><tr><td>1226.58</td><td>0.3490</td><td>28.26</td><td>1634.6</td><td>176.5</td><td>360.2</td></tr><tr><td>1674.45</td><td>2.0300</td><td>19.22</td><td>2441.5</td><td>150.5</td><td>238.7</td></tr><tr><td>3628.20</td><td>11.33</td><td>21.23</td><td>3875.0</td><td>134.5</td><td>273.4</td></tr><tr><td>7700.55</td><td>60.70</td><td>23.63</td><td>8052.1</td><td>156.2</td><td>329.6</td></tr></tbody></table>

On the other hand, relevant timings are reported in Table 2. First of all, we can notice that while the ratio between the iteration number for $N_S^{-1}S$ and $P^{-1}A$ is around 5, the ratio between solution times for $N_S^{-1}S$ and $P^{-1}A$ is around 2. Indeed, even though the latter requires more iterations, a problem with $A$ has to be solved iteratively every time $S$ has to be applied, increasing the solution time. Then, one can notice that the time to build matrix $E_0$ is the largest portion of the total setup time. In particular it is around 62% of the total setup for the smaller tests to pass to the 94% for the largest ones. We can notice that this time almost doubles at every test. This is because since the size of $V_0$ doubles as well. A more optimised building of $E_0$ will be the addressed in the future in order to decrease the time for setup. Moreover, factorization time of $E_0$ rapidly increases $E_0$ is a dense matrix. Finally, once $E_0$ is factorized, the time needed to apply (3.27) to some right hand side is remains constant, as expected.

### 5.3.2 Experiment using the method based on nested iterative solution

In this section we will perform the same weak scalability experiment using the method from [15]. So here $N_S^{-1}$ will be applied by solving a problem like (3.23). Numerical results are reported in Tables 3 and 4.

The sizes of coarse spaces and the iterations for matrix $A$ are reported again for a matter of clarity, even though they are exactly the same as in the previous experiment. Indeed, the way $N_S^{-1}$ is implemented does not affect the iteration count for $M_{HSM}^{-1}A$ nor the GEPs. The 7-th column of Table 3 contains the average number of iterations needed to solve (3.23) at every application of $N_S^{-1}$. We can notice that this value remains constant throughout the experiment. This allows us to say that for this specific problem $M_{S_{1,N}}^{-1}$ is a scalable preconditioner for $N_S$. The number of iterations to solve $N_S^{-1}S$ is the same in the two experiments. Indeed, the result in Lemma 3.6 holds true for both methods. Also the iteration count for $P^{-1}A$ is the same between the two experiments.

Finally, in Table 4 timings are reported. The total setup time is lower in this experiment. In this case, the largest portion of the setup is the solution of GEPs, which can be done fully in
parallel. Indeed, the setup time is scalable throughout the experiment, while it was not using the Woodbury formula. On the other hand, the GMRES time is higher in this case due to the inner iterative solution to perform. In particular, this is evident when solving the original problem with block preconditioning (last column of Table 4) where the iterations required are more.

<table><caption>Table 3: Weak scalability experiment for the heterogeneous linear elasticity problem using the original method from [15]. Iteration counts and timings.</caption><thead><tr><th>N (or P)</th><th>n+m</th><th>dimV<sub>0</sub></th><th>dimW<sub>0</sub></th><th>It. A</th><th>It. S</th><th>It. N<sub>S</sub></th><th>It. A</th></tr></thead><tbody><tr><td>16</td><td>1,058,188</td><td>726</td><td>79</td><td>46</td><td>17</td><td>13</td><td>84</td></tr><tr><td>32</td><td>2,110,941</td><td>1473</td><td>285</td><td>48</td><td>18</td><td>14.4</td><td>89</td></tr><tr><td>64</td><td>4,216,392</td><td>3088</td><td>914</td><td>49</td><td>18</td><td>15.5</td><td>91</td></tr><tr><td>128</td><td>8,427,373</td><td>6185</td><td>1726</td><td>52</td><td>19</td><td>16.1</td><td>97</td></tr><tr><td>256</td><td>16,849,298</td><td>12551</td><td>3476</td><td>58</td><td>19</td><td>15.9</td><td>107</td></tr></tbody></table>

<table><caption>Table 4: Weak scalability experiment for the heterogeneous linear elasticity problem using the original method from [15]. Iteration counts and timings.</caption><thead><tr><th>Total setup (s)</th><th>GMRES S (s)</th><th>GMRES A (s)</th></tr></thead><tbody><tr><td>370.8</td><td>290.1</td><td>2301.7</td></tr><tr><td>406.5</td><td>367.5</td><td>2887.1</td></tr><tr><td>364.9</td><td>304.3</td><td>1928.0</td></tr><tr><td>317.1</td><td>278.1</td><td>2288.4</td></tr><tr><td>343.4</td><td>299.2</td><td>2583.3</td></tr></tbody></table>

## 5.4 Comparison with norm based preconditioner

In this section we want to compare the performances of $N_S^{-1}$ as a preconditioner for $S$ with norm based preconditioner defined presented in [21]. In the latter, the Gram matrix corresponding to the norm on the dual space used to prove the inf-sup condition is considered as a preconditioner for $S$. In the case of heterogeneous elasticity one can prove an inf-sup condition using the $L^2$ discrete norm, [5]. The corresponding Gram matrix is the mass matrix. We will perform a robustness study where we will change the material coefficients of the elasticity problem by considering $\nu_2 = 0.49, 0.499, 0.4999$. The closer $\nu_2$ is to 0.5, the higher is the jump in the coefficients which means more difficult problem.

<figure><img src="image_4.png" /><figcaption>Figure 3: Residual plots from multiple log-derived data files.</figcaption></figure>
In the most left picture in Figure 3 the number of iterations to solve the problem with $S$ using the mass matrix $M$ as a preconditioner for increasing $\nu_2$ are reported. One can observe that the iteration count is very large, but it does not vary much for $\nu_2$ increasing. This means that the mass matrix is a robust preconditioner in this case, but it is not suitable due to the large number of iterations.

The central picture in Figure 3 represents the iteration count using the rescaled mass matrix $\lambda^{-1}M$ as a preconditioner. This is the most natural choice being equal to matric $C$. Rescaling is an attempt to take into account the constants entering the spectral equivalence between $M$ and $S$. We observe a lower number of iterations for $\nu_2 = 0.49$ but then that the iteration number increases with $\nu_2$, which means that the rescaled mass matrix is not a robust preconditioner.

Finally, if we consider $N_S^{-1}$ as preconditioner for $S$, we observe small iteration counts which remain fixed when the problem becomes difficult. This shows the robustness of our preconditioner with respect to the problem physical coefficient and moreover shows that for highly heterogenous problems a more sophisticated preconditioner than the (rescaled) mass matrix must be used.

# 6 Conclusion and outlooks

In this work we proposed a method to precondition the Schur complement matrix in a robust and scalable way for general saddle point linear systems. Such preconditioner can be used to solve saddle point systems iteratively, either with a segregated iteration algorithm where a problem for the Schur complement is solved or to build a block preconditioner for the original system.

To develop such preconditioner, some assumptions must be satisfied. We ask the diagonal blocks in the original system to be spectrally equivalent to sums of symmetric positive semi-definite matrices. This is easily verified e.g. in the FEM context. The presented method generalizes the one proposed in [15] and differs from it in the way the preconditioner is applied. The starting point of the argument is a two-level Schwarz method for the primal unknowns block.

For the original method, there is experimental evidence on the scalability of the inner solution needed to apply the preconditioner for the elasticity problem, but no theoretical results supporting it. Instead, the method presented in this work allows to take advantage of the spectral bounds we have for the Schur complement matrix to define a scalable and robust approximation. This, at the cost of building and factorizing a dense matrix which prevents to solve larger problems efficiently.

The perspective of an efficient and application-independent preconditioner that can be used to solve general saddle point problems with only few extra assumptions is fascinating. Further work is needed to make this method more efficient. First-of-all, a more optimised way to build the dense matrix must be considered. This will make the method of practical interest to tackle larger and more complex linear systems. Also interesting would be the introduction of multilevel methods to reduce the size of the dense matrix to factorize. Finally, how to decrease the iteration count when using block preconditioners by taking advantage of the DD context could be addressed as well.

# References

[1] H. Al Daas and L. Grigori. "A class of efficient locally constructed preconditioners based on coarse spaces". In: *SIAM J. Matrix Anal. Appl.* 40.1 (2019), pp. 66–91. doi: 10.1137/18M1194365.
[2] P. R. Amestoy, I. S. Duff, J.-Y. L'Excellent, and J. Koster. "A fully asynchronous multi-frontal solver using distributed dynamic scheduling". In: *SIAM J. Matrix Anal. Appl.* 23.1 (2001), pp. 15-41. DOI: 10.1137/S0895479899358194.

[3] O. Axelsson and M. Neytcheva. "Preconditioning methods for linear systems arising in constrained optimization problems". In: vol. 10. 1-2. Dedicated to the 60th birthday of Raytcho Lazarov. 2003, pp. 3-31. DOI: 10.1002/nla.310.

[4] M. Benzi, G. H. Golub, and J. Liesen. "Numerical solution of saddle point problems". In: *Acta Numer.* 14 (2005), pp. 1-137. DOI: 10.1017/S0962492904000212.

[5] D. Boffi, F. Brezzi, and M. Fortin. *Mixed Finite Element Methods and Applications*. Springer Series in Computational Mathematics. Springer Berlin Heidelberg, 2013.

[6] V. Dolean, P. Jolivet, and F. Nataf. *An Introduction to Domain Decomposition Methods: algorithms, theory and parallel implementation*. SIAM, 2015.

[7] M. Dryja and O. B. Widlund. "Schwarz methods of Neumann-Neumann type for three-dimensional elliptic finite element problems". In: *Comm. Pure Appl. Math.* 48.2 (1995), pp. 121-155. DOI: 10.1002/cpa.3160480203.

[8] Y. Efendiev, J. Galvis, R. Lazarov, and J. Willems. "Robust domain decomposition preconditioners for abstract symmetric positive definite bilinear forms". In: *ESAIM Math. Model. Numer. Anal.* 46.5 (2012), pp. 1175-1199. DOI: 10.1051/m2an/2011073.

[9] A. Greenbaum, V. Pták, and Z. Strakoš. "Any nonincreasing convergence curve is possible for GMRES". In: *SIAM J. Matrix Anal. Appl.* 17.3 (1996), pp. 465-469. DOI: 10.1137/S0895479894275030.

[10] M. Griebel and P. Oswald. "On the abstract theory of additive and multiplicative Schwarz algorithms". In: *Numer. Math.* 70.2 (1995), pp. 163-180. DOI: 10.1007/s002110050115.

[11] F. Hecht. "New development in FreeFem++". In: *J. Numer. Math.* 20.3-4 (2012), pp. 251-265.

[12] H. V. Henderson and S. R. Searle. "On Deriving the Inverse of a Sum of Matrices". In: *SIAM Review* 23.1 (1981), pp. 53-60.

[13] G. Karypis and V. Kumar. "METIS—A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes and Computing Fill-Reducing Ordering of Sparse Matrices". In: (1997).

[14] R. B. Lehoucq, D. C. Sorensen, and C. Yang. *ARPACK Users' Guide*. Society for Industrial and Applied Mathematics, 1998. DOI: 10.1137/1.9780898719628.

[15] F. Nataf and P.-H. Tournier. "A GenEO Domain Decomposition method for Saddle Point problems". In: *Comptes Rendus. Mécanique* 351.S1 (2023), pp. 1-18.

[16] F. Nataf, H. Xiang, V. Dolean, and N. Spillane. "A coarse space construction based on local Dirichlet-to-Neumann maps". In: *SIAM J. Sci. Comput.* 33.4 (2011), pp. 1623-1642. DOI: 10.1137/100796376.

[17] S. V. Nepomnyaschikh. "Mesh theorems on traces, normalizations of function traces and their inversion". In: *Soviet J. Numer. Anal. Math. Modelling* 6.3 (1991), pp. 223-242.

[18] Y. Notay. "A new analysis of block preconditioners for saddle point problems". In: *SIAM J. Matrix Anal. Appl.* 35.1 (2014), pp. 143-173.

[19] J. Pasciak and J. Zhao. "Overlapping Schwarz methods in H(curl) on polyhedral domains". In: *Journal of Numerical Mathematics* 10 (Jan. 2002). DOI: 10.1515/JNMA.2002.221.
[20] L. Pavarino and O. Widlund. "Balancing Neumann-Neumann Methods for Incompressible Stokes Equations". In: *Communications on Pure and Applied Mathematics* 55 (Mar. 2002), pp. 302-335. DOI: 10.1002/cpa.10020.

[21] J. Pestana and A. J. Wathen. "Natural preconditioning and iterative methods for saddle point systems". In: *SIAM Rev*. 57.1 (2015), pp. 71-91. DOI: 10.1137/130934921.

[22] N. Spillane, V. Dolean, H. Patrice, F. Nataf, C. Pechstein, et al. "Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps". In: *Numerische Mathematik* 126.4 (2014), pp. 741-770. DOI: 10.1007/s00211-013-0576-y.

[23] A. Toselli and O. Widlund. *Domain Decomposition Methods - Algorithms and Theory*. Vol. 34. Springer Series in Computational Mathematics. Springer, 2005.

[24] P.-H. Tournier, P. Jolivet, and F. Nataf. *FFDDM: FreeFem Domain Decomposition Method*. https://doc.freefem.org/documentation/ffddm/index.html. 2019.

[25] M. A. Woodbury. *Inverting modified matrices*. Statistical Research Group, Memo. Rep. no. 42, Princeton University, Princeton, NJ, 1950, p. 4.