# Schur complement based preconditioners for twofold and block tridiagonal
saddle point problems

Mingchao Cai Corresponding author, Department of Mathematics, Morgan State
University, Baltimore, MD 21251, USA, E-mail<EMAIL_ADDRESS>Guoliang Ju
School of Mechatronics Engineering, Harbin Institute of Technology, Harbin,
Heilongjiang 150001, China, and National Supercomputing Center in
Shenzhen,Shenzhen, Guangdong 518055, China, E-mail<EMAIL_ADDRESS>Jingzhi LI
Department of Mathematics, Southern University of Science and Technology,
Shenzhen, Guangdong 518055, China, E-mail<EMAIL_ADDRESS>
###### Abstract

In this paper, two types of Schur complement based preconditioners are studied
for twofold and block tridiagonal saddle point problems. One is based on the
nested (or recursive) Schur complement, the other is based on an additive type
Schur complement after permuting the original saddle point systems. We analyze
different preconditioners incorporating the exact Schur complements. It is
shown that some of them will lead to positive stable preconditioned systems.
Our discussion is instructive for devising various exact and inexact
preconditioners, as well as iterative solvers for many twofold and block
tridiagonal saddle point problems. Numerical experiments for a 3-field
formulation of Biot model is provided to justify our theoretical analysis.

Keywords Saddle point problems; Schur complement; block tridiagonal systems.

##  1 Introduction

Many application problems will lead to twofold and/or block tridiagonal saddle
point linear systems. Important examples include mixed formulations of the
Biot model [9, 15, 19, 24], the coupling of fluid flow with porous media flow
[6, 14], hybrid discontinuous Galerkin approximation of Stokes problem [11],
liquid crystal problem [3, 26] and optimization problems [16, 19, 21, 25, 27].
Some of these problems (or after permutations) will lead to a twofold saddle
point problem [3, 14, 16, 17, 27, 30] (or the so-called double saddle point
problem) of the following form.

$\mathcal{A}{\bm{x}}=\left[\begin{array}[]{ccc}A_{1}&B^{T}_{1}&0\\\\[2.84526pt]
C_{1}&-A_{2}&B^{T}_{2}\\\\[2.84526pt] 0&C_{2}&A_{3}\\\
\end{array}\right]\left[\begin{array}[]{c}{\bm{x}}_{1}\\\\[2.84526pt]
{\bm{x}}_{2}\\\\[2.84526pt]
{\bm{x}}_{3}\end{array}\right]=\left[\begin{array}[]{c}{\bm{f}}_{1}\\\\[2.84526pt]
{\bm{f}}_{2}\\\\[2.84526pt] {\bm{f}}_{3}\end{array}\right].$ (1)

A negative sign in front of $A_{2}$ is just for the ease of notation. After
simple permutations, the system matrix of (1) can be rewritten into the
following form.

$\mathcal{A}=\left[\begin{array}[]{ccc}A_{1}&0&B^{T}_{1}\\\\[2.84526pt]
0&A_{3}&C_{2}\\\\[2.84526pt] C_{1}&B_{2}^{T}&-A_{2}\\\ \end{array}\right].$
(2)

We call the system matrix in (2) is permutation-equivalent to that in (1).
Without causing confusion, we continue to use the notation $\mathcal{A}$ for
the permuted matrix (2). The linear system in (2) arises naturally from the
domain decomposition methods [22, 29]. In this work, we only assume that
$A_{1}$ is invertible and the global system matrix $\mathcal{A}$ is
invertible. Many special cases can be cast into the above forms of twofold
saddle point systems. For example, (a) $A_{2}=0$; (b) $A_{3}=0$; (c) both
$A_{2}=0$ and $A_{3}=0$. Our discussions will try to cover all these special
cases.

The above $3$-by-$3$ block linear problems (1) and (2) can be naturally
extended to the $n$-tuple cases. For example, when the system matrix in (1) is
extended to the $n$-tuple case, it is the block tridiagonal systems discussed
in [27]. When the system matrix in (2) is extended to the $n$-tuple case, it
corresponds to the linear system resulting from the domain decomposition
method for elliptic problem with $n-1$ subdomains. In many references, these
linear systems are assumed to be symmetric. No matter symmetric or not,
$\mathcal{A}$ is generally indefinite. For solving such a system in large-
scale computations, Krylov subspace methods with preconditioners are usually
applied. The analysis in [23, 18] indicates that one should employ Schur
complement based preconditioners [4] and then resort to replacing each exact
Schur complement with an approximate one. If a Krylov subspace method has
Galerkin property, the number of iterations depends on the degree of the
minimal polynomial that the preconditioned matrix $\mathcal{T}$ satisfies [23,
18, 28]. More precisely, the number of iterations depends on the dimension of
the Krylov subspace
$\mbox{Span}\\{{\bm{r}},\mathcal{T}{\bm{r}},\mathcal{T}^{2}{\bm{r}},\mathcal{T}^{3}{\bm{r}},\ldots\\}$.
For a classical $2$-by-$2$ saddle point problem, by employing block-diagonal
or block-triangular preconditioners with the $(2,2)$ block being the Schur
complement, it has been shown that the corresponding preconditioned system
satisfies a polynomial of the degree of $4$ or $2$ [23, 18]. Therefore,
devising preconditioners based on Schur complements is of vital importance.

In this work, we mainly discuss two approaches to design preconditioners for
$3$-by-$3$ block systems of block tridiagonal form and (2), as well as their
generalizations to $n$-tuple cases. One approach is based on the nested (or
recursive) Schur complement [27]. Unlike that in [27], we have a much general
assumption: $A_{1}$, $A_{2}$ and $A_{3}$ can be non-symmetric or negative
definite. Furthermore, when generalized to $n$-tuple block tridiagonal cases,
we investigate different preconditioners and emphasize the effects of the sign
in front of each Schur complement. We derive the polynomial equations of the
preconditioned systems and highlight the possible advantages and disadvantages
of each preconditioner. It is shown that choosing proper sign in front of each
Schur complement will result in a positive stable preconditioned system [10].
Another approach is based on the permuted matrix (2). One can firstly treat
the first 2-by-2 block as a whole (although maybe indefinite) and cast the
system into the saddle pointed form as that discussed in [18]. We then
investigate preconditioners based on an additive type Schur complement. If
block-triangular preconditioners are considered, the additive type Schur
complement based preconditioner will lead to a preconditioned system which
satisfies a polynomial of the degree 2 or 4, while the preconditioners based
on the nested Schur complement satisfies a polynomial with the degree 3 or 6.
When generalized to $n$-tuple cases, the degree of the polynomial that a
nested Schur complement based preconditioner satisfies is much higher than the
degree of the polynomial that an additive Schur complement based
preconditioner satisfies. Therefore, from a purely algebraic point of view, if
the GMRES method is applied, the additive type preconditioner is more favored.
When the form of (2) is generalized to $n$-tuple cases, it corresponds to the
linear system resulting from domain decomposition methods. Some results from
the literature of domain decomposition methods are referred to.

We started with a twofold saddle point problem and generalize our theory to
$n$-tuple block tridiagonal saddle point problems. It is shown that properly
selected sign in front of Schur complement based preconditioners will result
in positive stable preconditioned system [10]. Numerical experiments for a
3-field formulation of Biot model are provided to show that positive stable
preconditioners outperform other preconditioners. It is noteworthy that we do
not assume the positive definiteness of each $A_{i}$, the linear systems
studied in this work are actually not limited to saddle point problems. In the
concluding remarks, the potential advantages and disadvantages of the two
types of preconditioners will be discussed and highlighted, particularly, when
the corresponding inexact versions are applied.

The outline of the remainder of this paper is as follows. In section 2, we
briefly recall the classic saddle point problem and its Schur complement, and
introduces the twofold saddle point problem and the form of Schur complement
firstly, then we construct and analyze the block-triangular and block-diagonal
preconditioners based on Schur complement for twofold saddle point problems.
Furthermore, we extend these results to $n$-tuple saddle point problem in
Section 3. Some additive Schur complement based preconditioners are
constructed and analyzed in Section 4 for twofold saddle point problems.
Generalization to $n$-tuple cases are provided in Section 5. In Section 6,
numerical experiments for a 3-field Biot model is provided to justify the
advantages of using positive stable preconditioners. Finally, concluding
remarks are given in Section 7.

##  2 Nested Schur complement based preconditioners for twofold saddle point
systems

The matrix form of a classic saddle point problem [4, 18] reads as

$\mathcal{A}=\left[\begin{array}[]{ccc}A&B^{T}\\\ C&-D\end{array}\right].$

Following the notation in [27], the Schur complement of $\mathcal{A}$ is
defined as

$\mbox{Schur}(\mathcal{A})=-D-CA^{-1}B^{T}.$

In many references, $A$ is assumed to be symmetric positive definite, $D$ is
assumed to be symmetric and semi-positive definite, and $C$ is assumed to be
$B$. In this paper, we only make some assumptions that can guarantee the
invertibility of $\mathcal{A}$. We assume that $A$ and the Schur complement
$\mbox{Schur}(\mathcal{A})$ are invertible. For the twofold saddle point
problem (1), following the notations in [27], we denote

$S_{1}=A_{1},\quad S_{2}=A_{2}+C_{1}S_{1}^{-1}B_{1}^{T},$

and denote the nested Schur complement as

$S_{3}=A_{3}+C_{2}S_{2}^{-1}B_{2}^{T}.$

In [27], each $S_{i}$ is assumed to be symmetric positive definite (SPD). From
the expressions, if $\mathcal{A}$ is symmetric and all $A_{i}~{}i=1,\ldots,3$
are symmetric positive definite, then $S_{i}$ are SPD. In this work, we only
assume that each $S_{i}$ is invertible. Such assumptions are sufficient enough
to guarantee the invertibility of the global system.

We study both block-triangular and block-diagonal preconditioners for the
system matrix (1). For block-triangular preconditioners, we focus on a lower
triangular type with left preconditioning because an upper triangular one with
right preconditioning can be discussed in a similar way [4, 18]. We consider
the following preconditioner:

$\mathcal{P}_{1}=\left[\begin{array}[]{ccc}A_{1}&0&0\\\\[2.84526pt]
C_{1}&-S_{2}&0\\\\[2.84526pt] 0&C_{2}&S_{3}\end{array}\right].$ (3)

######  Theorem 1.

Assume that all $S_{i},i=1,2,3$, are invertible, then, for the preconditioner
(3),

$\mathcal{P}_{1}^{-1}\mathcal{A}=\left[\begin{array}[]{ccc}I&A_{1}^{-1}B_{1}^{T}&0\\\\[2.84526pt]
0&I&-S^{-1}_{2}B_{2}^{T}\\\\[2.84526pt] 0&0&I\end{array}\right],$

and $\mathcal{P}_{1}^{-1}\mathcal{A}$ satisfies the polynomial equation
$(\lambda-1)^{3}=0$.

###### Proof.

The block $LDU$ decomposition of the system matrix is

$\mathcal{A}=\left[\begin{array}[]{ccc}I&0&0\\\\[2.84526pt]
C_{1}A_{1}^{-1}&I&0\\\\[2.84526pt]
0&-C_{2}S_{2}^{-1}&I\end{array}\right]\left[\begin{array}[]{ccc}A_{1}&0&0\\\\[2.84526pt]
0&-S_{2}&0\\\\[2.84526pt]
0&0&S_{3}\end{array}\right]\left[\begin{array}[]{ccc}I&A_{1}^{-1}B_{1}^{T}&0\\\\[2.84526pt]
0&I&-S^{-1}_{2}B_{2}^{T}\\\\[2.84526pt] 0&0&I\end{array}\right].$ (4)

The multiplication of the $LD$ part leads to

$\mathcal{A}=\left[\begin{array}[]{ccc}A_{1}&0&0\\\\[2.84526pt]
C_{1}&-S_{2}&0\\\\[2.84526pt]
0&C_{2}&S_{3}\end{array}\right]\left[\begin{array}[]{ccc}I&A_{1}^{-1}B_{1}^{T}&0\\\\[2.84526pt]
0&I&-S^{-1}_{2}B_{2}^{T}\\\\[2.84526pt] 0&0&I\end{array}\right].$

Therefore, the conclusions of this theorem are obvious. ∎

Other preconditioners which will make the corresponding preconditioned systems
have eigenvalues $1$ or $-1$ are

$\mathcal{P}_{2}=\left[\begin{array}[]{ccc}A_{1}&0&0\\\ C_{1}&S_{2}&0\\\
0&-C_{2}&S_{3}\\\
\end{array}\right],\quad\mathcal{P}_{3}=\left[\begin{array}[]{ccc}A_{1}&0&0\\\
C_{1}&S_{2}&0\\\ 0&-C_{2}&-S_{3}\\\
\end{array}\right],\quad\mathcal{P}_{4}=\left[\begin{array}[]{ccc}A_{2}&0&0\\\\[2.84526pt]
C_{1}&-S_{2}&0\\\\[2.84526pt] 0&C_{2}&-S_{3}\end{array}\right].$ (5)

Remark 1. For preconditioners listed in (5), we have the following
conclusions.

$\mathcal{P}_{2}^{-1}\mathcal{A}=\left[\begin{array}[]{ccc}I&A_{1}^{-1}B_{1}^{T}&0\\\
0&-I&S^{-1}_{2}B_{2}^{T}\\\ 0&0&I\\\ \end{array}\right],$

and $\mathcal{P}_{2}^{-1}\mathcal{A}$ satisfies the polynomial equation
$(\lambda-1)^{2}(\lambda+1)=0$.

$\mathcal{P}_{3}^{-1}\mathcal{A}=\left[\begin{array}[]{ccc}I&A_{1}^{-1}B_{1}^{T}&0\\\
0&-I&S^{-1}_{2}B_{2}^{T}\\\ 0&0&-I\\\ \end{array}\right],$

and $\mathcal{P}_{3}^{-1}\mathcal{A}$ satisfies the polynomial equation
$(\lambda-1)(\lambda+1)^{2}=0$.

$\mathcal{P}_{4}^{-1}\mathcal{A}=\left[\begin{array}[]{ccc}I&A_{1}^{-1}B_{1}^{T}&0\\\\[2.84526pt]
0&I&-S^{-1}_{2}B_{2}^{T}\\\\[2.84526pt] 0&0&-I\end{array}\right],$

and $\mathcal{P}_{4}^{-1}\mathcal{A}$ satisfies the polynomial equation
$(\lambda-1)^{2}(\lambda+1)=0$.

######  Theorem 2.

Let the assumptions in Theorem 1 hold. We further assume that $A_{2}=0$ and
$A_{3}=0$. If the block-diagonal preconditioner is

$\mathcal{P}_{D_{1}}=\left[\begin{array}[]{ccc}A_{1}&0&0\\\ 0&S_{2}&0\\\
0&0&S_{3}\end{array}\right],$ (6)

then,

$\mathcal{T}_{1}=\mathcal{P}_{D_{1}}^{-1}\mathcal{A}=\left[\begin{array}[]{ccc}I&A_{1}^{-1}B^{T}_{1}&0\\\\[2.84526pt]
S_{2}^{-1}C_{1}&0&S_{2}^{-1}B^{T}_{2}\\\\[2.84526pt]
0&S_{3}^{-1}C_{2}&0\end{array}\right],$ (7)

which satisfies the polynomial equation
$(\lambda-1)(\lambda^{2}-\lambda-1)(\lambda^{3}-\lambda^{2}-2\lambda+1)=0$.

###### Proof.

By a simple calculation, one can derive that (7) holds true. Then, we see that

$\mathcal{T}_{1}^{2}=\left[\begin{array}[]{ccc}I+A_{1}^{-1}B^{T}_{1}S_{2}^{-1}C_{1}&A_{1}^{-1}B^{T}_{1}&A_{1}^{-1}B^{T}_{1}S_{2}^{-1}B^{T}_{2}\\\\[5.69054pt]
S_{2}^{-1}C_{1}&I+S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2}&0\\\\[5.69054pt]
S_{3}^{-1}C_{2}S_{2}^{-1}C_{1}&0&I\end{array}\right],$

$\mathcal{T}_{1}^{3}=\left[\begin{array}[]{ccc}I+2A_{1}^{-1}B^{T}_{1}S_{2}^{-1}C_{1}&A_{1}^{-1}B^{T}_{1}(2I+S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2})&A_{1}^{-1}B^{T}_{1}S_{2}^{-1}B^{T}_{2}\\\\[5.69054pt]
(2I+S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2})S_{2}^{-1}C_{1}&I&2S_{2}^{-1}B^{T}_{2}\\\\[5.69054pt]
S_{3}^{-1}C_{2}S_{2}^{-1}C_{1}&2S_{3}^{-1}C_{2}&0\end{array}\right],$

and

$\mathcal{T}_{1}-I=\left[\begin{array}[]{ccc}0&A_{1}^{-1}B^{T}_{1}&0\\\\[5.69054pt]
S_{2}^{-1}C_{1}&-I&S_{2}^{-1}B^{T}_{2}\\\\[5.69054pt] 0&S_{3}^{-1}C_{2}&-I\\\
\end{array}\right],$
$(\mathcal{T}_{1}^{2}-\mathcal{T}_{1}-I)=\left[\begin{array}[]{ccc}-I+A_{1}^{-1}B^{T}_{1}S_{2}^{-1}C_{1}&0&A_{1}^{-1}B^{T}_{1}S_{2}^{-1}B^{T}_{2}\\\\[5.69054pt]
0&S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2}&-S_{2}^{-1}B^{T}_{2}\\\\[5.69054pt]
S_{3}^{-1}C_{2}S_{2}^{-1}C_{1}&-S_{3}^{-1}C_{2}&0\end{array}\right],$
$\mathcal{T}_{1}^{3}-\mathcal{T}_{1}^{2}-2\mathcal{T}_{1}+I=\left[\begin{array}[]{ccc}-I+A_{1}^{-1}B^{T}_{1}S_{2}^{-1}C_{1}&A_{1}^{-1}B^{T}_{1}(-I+S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2})&0\\\\[5.69054pt]
(-I+S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2})S_{2}^{-1}C_{1}&I-S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2}&0\\\\[5.69054pt]
0&0&0\end{array}\right].$

Combining all the above together, one can verify that

$\begin{array}[]{ll}(\mathcal{T}_{1}^{2}-\mathcal{T}_{1}-I)(\mathcal{T}_{1}-I)=(\mathcal{T}_{1}-I)(\mathcal{T}_{1}^{2}-\mathcal{T}_{1}-I)=\\\\[11.38109pt]
\qquad\qquad\qquad\left[\begin{array}[]{ccc}0&A_{1}^{-1}B^{T}_{1}S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2}&-A_{1}^{-1}B^{T}_{1}S_{2}^{-1}B^{T}_{2}\\\\[5.69054pt]
S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2}S_{2}^{-1}C_{1}&-2S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2}&2S_{2}^{-1}B^{T}_{2}\\\\[5.69054pt]
-S_{3}^{-1}C_{2}S_{2}^{-1}C_{1}&2S_{3}^{-1}C_{2}&-I\end{array}\right].\end{array}$

The $(i,j)$ block of
$(\mathcal{T}_{1}-I)(\mathcal{T}_{1}^{2}-\mathcal{T}_{1}-I)(\mathcal{T}_{1}^{3}-\mathcal{T}_{1}^{2}-2\mathcal{T}_{1}+I)$
is denoted as $A_{ij}$. One can verify that

$\begin{array}[]{ll}A_{11}=A_{1}^{-1}B^{T}_{1}S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2}(-I+S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2})S_{2}^{-1}C_{1}=0,\\\\[5.69054pt]
A_{12}=A_{1}^{-1}B^{T}_{1}S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2}(I-S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2})=0,\\\\[5.69054pt]
A_{21}=S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2}S_{2}^{-1}C_{1}(-I+A_{1}^{-1}B^{T}_{1}S_{2}^{-1}C_{1})\\\\[5.69054pt]
\qquad\quad-2S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2}(-I+S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2})S_{2}^{-1}C_{1}=0,\\\\[5.69054pt]
A_{22}=S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2}S_{2}^{-1}C_{1}A_{1}^{-1}B^{T}_{1}(-I+S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2})\\\\[5.69054pt]
\qquad\quad-2S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2}(I-S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2})=0,\\\\[5.69054pt]
A_{31}=-S_{3}^{-1}C_{2}S_{2}^{-1}C_{1}(-I+A_{1}^{-1}B^{T}_{1}S_{2}^{-1}C_{1})\\\\[5.69054pt]
\qquad\quad+2S_{3}^{-1}C_{2}(-I+S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2})S_{2}^{-1}C_{1}=0,\\\\[5.69054pt]
A_{32}=-S_{3}^{-1}C_{2}S_{2}^{-1}C_{1}A_{1}^{-1}B^{T}_{1}(-I+S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2})\\\\[5.69054pt]
\qquad\quad+2S_{3}^{-1}C_{2}(I-S_{2}^{-1}B^{T}_{2}S_{3}^{-1}C_{2})=0,\\\\[5.69054pt]
A_{13}=0,\quad A_{23}=0,\quad A_{33}=0.\end{array}$

We therefore conclude that

$(\mathcal{T}_{1}-I)(\mathcal{T}_{1}^{2}-\mathcal{T}_{1}-I)(\mathcal{T}_{1}^{3}-\mathcal{T}_{1}^{2}-2\mathcal{T}_{1}+I)=0.$
(8)

∎

Remark 2. We note that (8) can be factorized into distinct linear factors.
Correspondingly, the roots of the above polynomial are:

$-1.2470,~{}-0.6180,~{}0.4450,~{}1,~{}1.6180,~{}1.8019.$

It follows that the eigenvalues of $\mathcal{T}_{1}$ must be among these
values. Theorem 2.3 in [27] shows that if $\mathcal{A}$ is a symmetric matrix,
then the eigenvalues of the preconditioned system $\mathcal{T}_{1}$ are

$\sigma_{p}(\mathcal{T}_{1})=\left\\{2\cos\left(\frac{2i-1}{2j+1}\pi\right)\colon
j=1,\ldots,3,\ i=1,\ldots,j\right\\},$

Correspondingly, the condition number is

$\kappa(\mathcal{T}_{1})\leq\frac{\cos\left(\frac{\pi}{7}\right)}{\sin\left(\frac{\pi}{14}\right)}\approx
4.05.$

Our conclusion for eigenvalues in Theorem 2 is consistent with that in [27].
However, the assumption that $\mathcal{A}$ is symmetric is removed in our
analysis.

We also consider the following Schur complement based block-diagonal
preconditioners:

$\mathcal{P}_{D_{2}}=\left[\begin{array}[]{ccc}A_{1}&0&0\\\ 0&S_{2}&0\\\
0&0&-S_{3}\end{array}\right],\quad\mathcal{P}_{D_{3}}=\left[\begin{array}[]{ccc}A_{1}&0&0\\\
0&-S_{2}&0\\\
0&0&S_{3}\end{array}\right],\quad\mathcal{P}_{D_{4}}=\left[\begin{array}[]{ccc}A_{1}&0&0\\\
0&-S_{2}&0\\\ 0&0&-S_{3}\end{array}\right].$ (9)

Furthermore, if $A_{2}=0$ and $A_{3}=0$, we have the following conclusions: If
the preconditioner is $\mathcal{P}_{D_{2}}$, then

$\mathcal{T}_{2}=\mathcal{P}_{D_{2}}^{-1}\mathcal{A}=\left[\begin{array}[]{ccc}I&A_{1}^{-1}B^{T}_{1}&0\\\
S_{2}^{-1}C_{1}&0&S_{2}^{-1}B^{T}_{2}\\\ 0&-S_{3}^{-1}C_{2}&0\\\
\end{array}\right],$

satisfies the polynomial equation
$(\lambda-1)(\lambda^{2}-\lambda-1)(\lambda^{3}-\lambda^{2}-1)=0$; Its roots
are

$-0.618,~{}~{}-0.2328+0.7926i,~{}~{}-0.2328-0.7926i,~{}~{}1,~{}~{}1.4656,~{}~{}1.618;$

If the preconditioner is $\mathcal{P}_{D_{3}}$, then

$\mathcal{T}_{3}=\mathcal{P}_{D_{3}}^{-1}\mathcal{A}=\left[\begin{array}[]{ccc}I&A_{1}^{-1}B^{T}_{1}&0\\\
-S_{2}^{-1}C_{1}&0&-S_{2}^{-1}B^{T}_{2}\\\ 0&S_{3}^{-1}C_{2}&0\\\
\end{array}\right],$

which satisfies the polynomial equation
$(\lambda-1)(\lambda^{2}-\lambda+1)(\lambda^{3}-\lambda^{2}+2\lambda-1)=0$;
Correspondingly, the roots are

$0.5698,~{}~{}1,~{}~{}0.5+0.8660i,~{}~{}0.5-0.8660i,~{}~{}0.2151+1.3071i,~{}~{}0.2151-1.3071i;$

If the preconditioner is $\mathcal{P}_{D_{4}}$, then

$\mathcal{T}_{4}=\mathcal{P}_{D_{4}}^{-1}\mathcal{A}=\left[\begin{array}[]{ccc}I&A_{1}^{-1}B^{T}_{1}&0\\\
-S_{2}^{-1}C_{1}&0&-S_{2}^{-1}B^{T}_{2}\\\ 0&-S_{3}^{-1}C_{2}&0\\\
\end{array}\right],$

which satisfies the polynomial equation
$(\lambda-1)(\lambda^{2}-\lambda+1)(\lambda^{3}-\lambda^{2}+1)=0$; The roots
are

$-0.7549,~{}~{}1,~{}~{}0.5+0.8660i,~{}~{}0.5-0.8660i,~{}~{}0.8774+0.7449i,~{}~{}0.8774-0.7449i.$

We highlight here that all eigenvalues of $\mathcal{T}_{3}$ have positive real
parts. We call $\mathcal{P}_{D_{3}}$ a positive stable preconditioner for
$\mathcal{A}$.

##  3 Nested Schur complement based preconditioners for $n$-tuple block
tridiagonal problems

Now, we will discuss preconditioners for the $n$-tuple block tridiagonal
problem of the following form.

$\mathcal{A}_{n}=\left[\begin{array}[]{cccc}A_{1}&B_{1}^{T}&&\\\
C_{1}&-A_{2}&\ddots&\\\ &\ddots&\ddots&B_{n-1}^{T}\\\\[4.30554pt]
&&C_{n-1}&(-1)^{n-1}A_{n}\end{array}\right].$ (10)

Similar to the $3$-by-$3$ block case, we define

$S_{i+1}=A_{i+1}+C_{i}S^{-1}_{i}B_{i}^{T},\quad\text{for}\quad
i=1,2,\ldots,n-1,$

with the initial setting $S_{1}=A_{1}$. The block tridiagonal system can be
written recursively as

$\mathcal{A}_{i+1}=\left[\begin{array}[]{cccc}&\mathcal{A}_{i}&\mathbb{B}_{i}^{T}\\\
&\mathbb{C}_{i}&(-1)^{i}A_{i+1}\end{array}\right].$

Then,

$\begin{array}[]{ll}\mbox{Schur}(\mathcal{A}_{i+1})&=(-1)^{i}A_{i+1}-\mathbb{C}_{i}\mathcal{A}_{i}^{-1}\mathbb{B}_{i}^{T}\\\\[5.69054pt]
&=(-1)^{i}(A_{i+1}+C_{i}S_{i}^{-1}B_{i}^{T}).\end{array}$

Moreover,
$S_{i+1}=(-1)^{i}\mbox{Schur}(\mathcal{A}_{i+1})=A_{i+1}+C_{i}S_{i}^{-1}B_{i}^{T}$
with $S_{1}=A_{1}$.

Similar to the factorization (4), $\mathcal{A}_{n}$ can be factorized into

$\mathcal{A}_{n}=\mathcal{L}_{n}\mathcal{D}_{n}\mathcal{U}_{n},$ (11)

where

$\mathcal{D}_{n}=\left[\begin{array}[]{cccc}A_{1}&0&&\\\ 0&-S_{2}&\ddots&\\\
&\ddots&\ddots&0\\\\[4.30554pt] &&0&(-1)^{n-1}S_{n}\end{array}\right],$ (12)

and

$\begin{array}[]{ll}\mathcal{L}_{n}=\left[\begin{array}[]{cccc}I&0&&\\\
C_{1}A_{1}^{-1}&I&\ddots&\\\ &\ddots&\ddots&0\\\\[4.30554pt]
&&(-1)^{n-2}C_{n-1}S^{-1}_{n-1}&I\end{array}\right]\quad\mathcal{U}_{n}=\left[\begin{array}[]{cccc}I&A_{1}^{-1}B_{1}^{T}&&\\\
0&I&\ddots&\\\ &\ddots&\ddots&(-1)^{n-2}S^{-1}_{n-1}B_{n-1}^{T}\\\\[4.30554pt]
&&0&I\end{array}\right].\end{array}$

######  Theorem 3.

Assume that all $S_{i},i=1,2,3,\ldots n$, are invertible. Let the
preconditioner be

$\mathcal{P}_{n}=\left[\begin{array}[]{cccc}A_{1}&\\\ C_{1}&-S_{2}\\\
&\ddots&\ddots\\\ &&C_{n-1}&(-1)^{n-1}S_{n}\end{array}\right],$

then

$\mathcal{P}_{n}^{-1}\mathcal{A}_{n}=\left[\begin{array}[]{cccc}I&A_{1}^{-1}B_{1}^{T}&&\\\
0&I&\ddots&\\\ &\ddots&\ddots&(-1)^{n-2}S^{-1}_{n-1}B_{n-1}^{T}\\\\[4.30554pt]
&&0&I\end{array}\right],$

and $\mathcal{P}_{n}^{-1}\mathcal{A}_{n}$ satisfies the polynomial
$(\lambda-1)^{n}=0$. Furthermore, if $\mathcal{A}_{n}$ is symmetric and
$A_{i},~{}i=1,2,\ldots,n$ are symmetric positive definite, let the
preconditioner be $\mathcal{D}_{n}$, then all eigenvalues of
$\mathcal{D}_{n}^{-1}\mathcal{A}_{n}$ are positive real.

###### Proof.

Combining $\mathcal{L}_{n}\mathcal{D}_{n}$ in (11), the first part of the
conclusion is obvious.

For the second part, noting that if $\mathcal{A}_{n}$ is symmetric and
$A_{i},~{}i=1,2,\ldots,n$, are symmetric positive definite, then the matrix

$\mathcal{S}_{n}=\left[\begin{array}[]{cccc}S_{1}&&&\\\ &S_{2}&&\\\
&&\ddots&\\\ &&&S_{n}\end{array}\right]$ (13)

is symmetric positive definite [27]. We therefore can define
$\mathcal{S}_{n}^{\frac{1}{2}}$, in which the $i$-th diagonal block is
$S_{i}^{\frac{1}{2}}$. To prove the eigenvalues of
$\mathcal{D}_{n}^{-1}\mathcal{A}_{n}$ are positive real, we apply a similar
transform to $\mathcal{D}_{n}^{-1}\mathcal{A}_{n}$, that is, left multiplying
by $\mathcal{S}_{n}^{\frac{1}{2}}$ and right multiplying by
$\mathcal{S}_{n}^{-\frac{1}{2}}$ to $\mathcal{D}_{n}^{-1}\mathcal{A}_{n}$ , we
have

$\displaystyle\mathcal{S}_{n}^{\frac{1}{2}}\mathcal{D}_{n}^{-1}\mathcal{A}_{n}\mathcal{S}_{n}^{-\frac{1}{2}}$
$\displaystyle=\left[\begin{array}[]{cccc}I&S_{1}^{-\frac{1}{2}}B_{1}^{T}S_{2}^{-\frac{1}{2}}&&\\\
-S_{2}^{-\frac{1}{2}}B_{1}S_{1}^{-\frac{1}{2}}&S_{2}^{-\frac{1}{2}}A_{2}S_{2}^{-\frac{1}{2}}&\ddots&\\\
&\ddots&\ddots&(-1)^{n-2}S_{n-1}^{-\frac{1}{2}}B_{n-1}^{T}S_{n}^{-\frac{1}{2}}\\\\[4.30554pt]
&&(-1)^{n-1}S_{n}^{-\frac{1}{2}}B_{n-1}S_{n-1}^{-\frac{1}{2}}&S_{n}^{-\frac{1}{2}}A_{n}S_{n}^{-\frac{1}{2}}\end{array}\right].$
(18)

We note that the Hermitian part of
$\mathcal{S}_{n}^{\frac{1}{2}}\mathcal{D}_{n}^{-1}\mathcal{A}_{n}\mathcal{S}_{n}^{-\frac{1}{2}}$
is

$\displaystyle\frac{1}{2}\left[(\mathcal{S}_{n}^{\frac{1}{2}}\mathcal{D}_{n}^{-1}\mathcal{A}_{n}\mathcal{S}_{n}^{-\frac{1}{2}})+(\mathcal{S}_{n}^{\frac{1}{2}}\mathcal{D}_{n}^{-1}\mathcal{A}_{n}\mathcal{S}_{n}^{-\frac{1}{2}})^{T}\right]$
$\displaystyle=\left[\begin{array}[]{cccc}I&0&&\\\
0&S_{2}^{-\frac{1}{2}}A_{2}S_{2}^{-\frac{1}{2}}&\ddots&\\\
&\ddots&\ddots&0\\\\[4.30554pt]
&&0&S_{n}^{-\frac{1}{2}}A_{n}S_{n}^{-\frac{1}{2}}\end{array}\right].$ (23)

is symmetric positive definite. Therefore, we conclude that all eigenvalues of
$\mathcal{D}_{n}^{-1}\mathcal{A}_{n}$ are positive real.

∎

Remark 3. In [27], $\mathcal{M}_{n}$ is taken as a preconditioner of the block
tridiagonal system $\mathcal{A}_{n}$. The advantage is that $\mathcal{M}_{n}$
is SPD. However, the eigenvalues of $\mathcal{M}_{n}^{-1}\mathcal{A}_{n}$ may
not be positive real. In comparison, $\mathcal{D}_{n}$ is a positive stable
preconditioner for $\mathcal{A}_{n}$.

######  Theorem 4.

Assume that $A_{i}=0$ for $i=2,3,\ldots,n$ in $\mathcal{A}_{n}$, if the Schur
complement based preconditioner is $\mathcal{D}_{n}$ in (12), then
$\mathcal{T}_{n}=\mathcal{D}_{n}^{-1}\mathcal{A}_{n}$ satisfies a polynomial
equation:

$\prod^{n}_{i=0}\bar{p}_{i}=0,$ (24)

where $\bar{p}_{0}(\lambda)=1,~{}\bar{p}_{1}(\lambda)=\lambda-1$ and

$\bar{p}_{i}=\lambda\bar{p}_{i-1}+\bar{p}_{i-2},~{}i=2,3,4,\ldots,n.$

Moreover, all the roots of $\bar{p}_{i},~{}i=1,2,\ldots,n$, lie in the right
half-plane.

###### Proof.

To verify (24), we let

$\begin{array}[]{ll}\mathcal{C}_{i}=S^{-1}_{i+1}C_{i},\quad\mathcal{B}^{T}_{i}=S^{-1}_{i}B_{i}^{T},\quad\text{for
}\,i=1,2,\ldots,n-1.\end{array}$

Because $A_{i}=0$ for $i=2,\ldots,n$, we see that
$\mathcal{C}_{i}\mathcal{B}^{T}_{i}=I$ and the preconditioned system
$\mathcal{T}_{n}$ degenerates to

$\mathcal{T}_{n}=\mathcal{D}_{n}^{-1}\mathcal{A}_{n}=\left[\begin{array}[]{cccc}I&\mathcal{B}^{T}_{1}&&\\\
-\mathcal{C}_{1}&0&\ddots&\\\
&\ddots&\ddots&(-1)^{n-2}\mathcal{B}^{T}_{n-1}\\\
&&(-1)^{n-1}\mathcal{C}_{n-1}&0\end{array}\right].$ (25)

Let $\lambda$ be an eigenvalue and ${\bf x}$ be the corresponding eigenvector,
we have $\mathcal{T}_{n}\mathbf{x}=\lambda\,\mathbf{x}$. More clearly,

$\displaystyle x_{1}+\mathcal{B}^{T}_{1}x_{2}$ $\displaystyle=\lambda x_{1},$
$\displaystyle\mathcal{C}_{1}x_{1}+\mathcal{B}^{T}_{2}x_{3}$
$\displaystyle=-\lambda x_{2},$ $\displaystyle\vdots$ (26)
$\displaystyle\mathcal{C}_{n-2}x_{n-2}+\mathcal{B}^{T}_{n-1}x_{n}$
$\displaystyle=(-1)^{n-2}\lambda x_{n-1},$
$\displaystyle\mathcal{C}_{n-1}x_{n-1}$ $\displaystyle=(-1)^{n-1}\lambda
x_{n}.$

By using the first equation of (26),

$\mathcal{B}^{T}_{1}x_{2}=\bar{p}_{1}(\lambda)\,x_{1},\quad\text{where}\quad\bar{p}_{1}(\lambda)=\lambda-1.$

We note that $\lambda_{11}=1$, the root of $\bar{p}_{1}(\lambda)$, is an
eigenvalue by setting $x_{2}=x_{3}=\ldots=x_{n}=0$. If
$\lambda\neq\lambda_{11}$, from the second equation and replacing $x_{1}$ by
$\frac{1}{\bar{p}_{1}(\lambda)}\mathcal{B}^{T}_{1}x_{2}$, we then obtain

$\mathcal{B}^{T}_{2}x_{3}=-\mathcal{C}_{1}x_{1}-\lambda\,x_{2}=-\frac{1}{\bar{p}_{1}(\lambda)}\,\mathcal{C}_{1}\mathcal{B}_{1}^{T}x_{2}-\lambda\,x_{2}=-\frac{1}{\bar{p}_{1}(\lambda)}\,x_{2}-\lambda\,x_{2}=-R_{2}(\lambda)\,x_{2},$

where

$R_{2}(\lambda)=\lambda+\frac{1}{\bar{p}_{1}(\lambda)}=\frac{\bar{p}_{2}(\lambda)}{\bar{p}_{1}(\lambda)}\quad\text{with}\quad\bar{p}_{2}(\lambda)=\lambda\bar{p}_{1}(\lambda)+1=\lambda\bar{p}_{1}(\lambda)+\bar{p}_{0}(\lambda).$

By setting $x_{3}=\ldots=x_{n}=0$, we see that the two roots of the polynomial
$\bar{p}_{2}(\lambda)$, denoted as $\lambda_{21}$ and $\lambda_{22}$, are
eigenvalues. Repeating this procedure leads to

$\displaystyle\mathcal{B}^{T}_{j}x_{j+1}=(-1)^{j-1}R_{j}(\lambda)\,x_{j},\text{
for }j=2,\ldots,n-1,\text{ and }0=(-1)^{n-1}R_{n}(\lambda)\,x_{n}\text{ with
}R_{j}(\lambda)=\frac{\bar{p}_{j}(\lambda)}{\bar{p}_{j-1}(\lambda)},$

where the polynomials $\bar{p}_{j}(\lambda)$ are recursively given by

$\bar{p}_{0}(\lambda)=1,\quad\bar{p}_{1}(\lambda)=\lambda-1,\quad\bar{p}_{i+1}(\lambda)=\lambda\,\bar{p}_{i}(\lambda)+\bar{p}_{i-1}(\lambda)\quad\text{for}\,i\geq
1.$ (27)

Therefore, the eigenvalues of $\mathcal{T}_{n}$ are the roots of the
polynomials
$\bar{p}_{1}\left(\lambda\right),\bar{p}_{2}(\lambda),\ldots,\bar{p}_{n}\left(\lambda\right)$.

For proving the second part of the conclusion, we note that all eigenvalues of
$\mathcal{D}^{-1}_{n}\mathcal{A}_{n}$ are the roots of
$\bar{p}_{i},~{}i=1,2,\ldots,n$. We are going to show that all the roots of
polynomial $\bar{p}_{i},i=1,2,\ldots,n$, lie on the right half-plane by using
the Routh-Hurwitz stability criterion [5, 13].

The stability of a polynomial can be determined by the Routh-Hurwitz stability
criterion which suggests using a tabular method [5]: For a general $k$th-
degree polynomial

$p_{k}(\lambda)=a_{k}\lambda^{k}+a_{k-1}\lambda^{k-1}+\cdots+a_{1}\lambda+a_{0},$

its Routh table is a matrix of $k+1$ rows, which has the following structure:

$\begin{array}[]{lllllllllllllll}\lambda^{k}&r_{01}&r_{02}&r_{03}&\cdots\\\\[2.84526pt]
\lambda^{k-1}&r_{11}&r_{12}&r_{13}&\cdots\\\\[2.84526pt]
\lambda^{k-2}&r_{21}&r_{22}&r_{23}&\cdots\\\\[2.84526pt]
\vdots&\vdots&\vdots&\vdots&\ddots\end{array}.$ (28)

In (28), we add a column $\lambda^{i},~{}i=k,k-1,\ldots 0$, as the row marker.
The first two rows, which come from the coefficients directly, are ordered as

$(r_{0j})=(a_{k},a_{k-2},a_{k-4},\cdots),\quad(r_{1j})=(a_{k-1},a_{k-3},a_{k-5},\cdots).$
(29)

Other entries are computed as follows:

$\displaystyle
r_{ij}=\frac{-1}{r_{i-1,1}}\mbox{det}\begin{bmatrix}r_{i-2,1}&r_{i-2,j+1}\\\
r_{i-1,1}&r_{i-1,j+1}\end{bmatrix},\quad i=2,3,4,\cdots.$

When the table is completed, the number of sign changes in the first column is
the number of non-negative roots [5]. Note that the entries of the Routh table
are uniquely determined once the first two rows are fixed, i.e., the Routh
table is uniquely dependent on the coefficients of $p_{k}$.

For the polynomial $\bar{p}_{k}$ generated by (27), its key coefficients (for
the leading terms and the constant) are

$\bar{a}_{k}=1,\quad\bar{a}_{k-1}=-1,\quad\bar{a}_{k-2}=k-1,\bar{a}_{k-3}=-(k-2),\quad\bar{a}_{0}=(-1)^{k}.$
(30)

Such a conclusion can be proved by the method of mathematical induction. The
proof is omitted here. From (30), in general, we have.

$\bar{p}_{k}=\lambda^{k}-\lambda^{k-1}+\bar{a}_{k-2}\lambda^{k-2}+\bar{a}_{k-3}\lambda^{k-3}+\bar{a}_{k-4}\lambda^{k-4}\ldots+\bar{a}_{2}\lambda^{2}+\bar{a}_{1}\lambda+\bar{a}_{0}.$

We are trying to write out the expression of $\bar{p}_{k-1}$ in terms of the
coefficients of $\bar{p}_{k}$. It is notable that the expression of
$\bar{p}_{k-1}$ will have small differences for the case $n$ odd and the case
$n$ even:

$\begin{array}[]{l}\displaystyle\bar{p}_{k-1}=\lambda^{k-1}-\lambda^{k-2}-\bar{a}_{k-3}\lambda^{k-3}-(\bar{a}_{k-4}+\bar{a}_{k-5})\lambda^{k-4}-\bar{a}_{k-5}\lambda^{k-5}-(\bar{a}_{k-6}+\bar{a}_{k-7})\lambda^{k-6}\\\
\displaystyle\qquad\quad\ldots-(\bar{a}_{5}+\bar{a}_{4})\lambda^{5}-\bar{a}_{4}\lambda^{4}-(\bar{a}_{3}+\bar{a}_{2})\lambda^{3}-\bar{a}_{2}\lambda^{2}-(\bar{a}_{1}+\bar{a}_{0})\lambda-\bar{a}_{0},\qquad\mbox{if
$k$ is odd},\\\\[5.69054pt]
\displaystyle\bar{p}_{k-1}=\lambda^{k-1}-\lambda^{k-2}-\bar{a}_{k-3}\lambda^{k-3}-(\bar{a}_{k-4}+\bar{a}_{k-5})\lambda^{k-4}-\bar{a}_{k-5}\lambda^{k-5}-(\bar{a}_{k-6}+\bar{a}_{k-7})\lambda^{k-6}\\\
\displaystyle\qquad\quad\ldots-\bar{a}_{5}\lambda^{5}-(\bar{a}_{4}+\bar{a}_{3})\lambda^{4}-\bar{a}_{3}\lambda^{3}-(\bar{a}_{2}+\bar{a}_{1})\lambda^{2}-\bar{a}_{1}\lambda-1,\qquad\qquad\quad~{}~{}\mbox{if
$k$ is even}.\end{array}$ (31)

For ease of presentation, let us assume that $k$ is even ($k-1$ and $k+1$ are
odd. The case $k$ is odd can be proved similarly by slightly changing the
notations in the following arguments). We denote the matrix $R_{T_{k}}$ as the
Routh table of polynomial $\bar{p}_{k}$. The first two rows of $R_{T_{k}}$ are

$\begin{array}[]{rrrrrrrrrrrrr}\lambda^{k}&\bar{a}_{k}=1&\bar{a}_{k-2}&\bar{a}_{k-4}&\ldots&\bar{a}_{2}&\bar{a}_{0}=1\\\
\lambda^{k-1}&\bar{a}_{k-1}=-1&\bar{a}_{k-3}&\bar{a}_{k-5}&\ldots&\bar{a}_{1}&0,\\\
\end{array}$ (32)

then, the first two rows of $R_{T_{k-1}}$ are

$\begin{array}[]{rrrrrrrrrrrrr}\lambda^{k-1}&1&-\bar{a}_{k-3}&-\bar{a}_{k-5}&\ldots&-\bar{a}_{1}\\\
\lambda^{k-2}&-1&-(\bar{a}_{k-4}+\bar{a}_{k-5})&-(\bar{a}_{k-6}+\bar{a}_{k-7})&\ldots&-1;\\\
\end{array}$ (33)

the first two rows of $R_{T_{k+1}}$ are

$\begin{array}[]{rrrrrrrrrrrrrr}\lambda^{k+1}&1&\bar{a}_{k-2}-\bar{a}_{k-1}&\bar{a}_{k-4}-\bar{a}_{k-3}&\ldots&\bar{a}_{2}-\bar{a}_{3}&\bar{a}_{0}-\bar{a}_{1}\\\
\lambda^{k}&-\bar{a}_{k}=-1&-\bar{a}_{k-2}&-\bar{a}_{k-4}&\ldots&-\bar{a}_{2}&-\bar{a}_{0},\\\
\end{array}$ (34)

and $R_{T_{k+1}}$ satisfy the following relationship.

$\begin{array}[]{rrrrrrrrrrrrrr}\lambda^{k+1}&1&\bar{a}_{k-2}-\bar{a}_{k-1}&\bar{a}_{k-4}-\bar{a}_{k-3}&\ldots&\bar{a}_{2}-\bar{a}_{3}&\bar{a}_{0}-\bar{a}_{1}\\\
\lambda^{k}\\\ \vdots&&&-R_{T_{k}}\\\ \lambda^{1}\\\ \lambda^{0}\end{array}.$
(35)

We actually only need to verify that the derivation from (32) to (34) is
correct. To check this, by using (27), one obtains

$\begin{array}[]{lll}&\bar{p}_{k+1}=\lambda\bar{p}_{k}+\bar{p}_{k-1}\\\
&=\lambda(\lambda^{k}-\lambda^{k-1}+\bar{a}_{k-2}\lambda^{k-2}+\bar{a}_{k-3}\lambda^{k-3}+\bar{a}_{k-4}\lambda^{k-4}\ldots+\bar{a}_{2}\lambda^{2}+\bar{a}_{1}\lambda+\bar{a}_{0})\\\
&+\lambda^{k-1}-\lambda^{k-2}-\bar{a}_{k-3}\lambda^{k-3}-(\bar{a}_{k-4}+\bar{a}_{k-5})\lambda^{k-4}-\bar{a}_{k-5}\lambda^{k-5}-(\bar{a}_{k-6}+\bar{a}_{k-7})\lambda^{k-6}\ldots\\\
&-(\bar{a}_{5}+\bar{a}_{4})\lambda^{5}-\bar{a}_{4}\lambda^{4}-(\bar{a}_{3}+\bar{a}_{2})\lambda^{3}-\bar{a}_{2}\lambda^{2}-(\bar{a}_{1}+\bar{a}_{0})\lambda-\bar{a}_{0}\\\
&=\lambda^{k+1}-\lambda^{k}+(\bar{a}_{k-2}+1)\lambda^{k-1}+(\bar{a}_{k-3}-1)\lambda^{k-2}+(\bar{a}_{k-4}-\bar{a}_{k-3})\lambda^{k-3}-\bar{a}_{k-4}\lambda^{k-4}\\\
&\quad+(\bar{a}_{k-6}-\bar{a}_{k-5})\lambda^{k-5}-\bar{a}_{k-6}\lambda^{k-6}\ldots+(\bar{a}_{4}-\bar{a}_{3})\lambda^{2}-\bar{a}_{3}\lambda^{3}+(\bar{a}_{1}-\bar{a}_{2})\lambda^{2}-\bar{a}_{1}\lambda+1.\end{array}$

As $k$ is even, the expression of $\bar{p}_{k+1}$ satisfies the equation in
(31), and the first two rows of the Routh table of $\bar{p}_{k+1}$ are

$\begin{array}[]{rrrrrrrrrrrrrr}\lambda^{k+1}&1&\bar{a}_{k-2}-\bar{a}_{k-1}&\bar{a}_{k-4}-\bar{a}_{k-3}&\ldots&\bar{a}_{2}-\bar{a}_{3}&\bar{a}_{0}-\bar{a}_{1}\\\
\lambda^{k}&-\bar{a}_{k}=-1&-\bar{a}_{k-2}&-\bar{a}_{k-4}&\ldots&-\bar{a}_{2}&-\bar{a}_{0}.\\\
\end{array}$

Thus, (34) is verified. As the entries of the Routh table are uniquely
determined by the first two rows, the rest of the verifications are
straightforward: for verifying (35), one can directly check the entries row-
by-row by using (29); Actually, one only needs to check the $3$-rd row:

$r_{21}=\frac{-1}{-1}\mbox{det}\begin{bmatrix}1&\bar{a}_{k-2}-\bar{a}_{k-1}\\\
-1&-\bar{a}_{k-2}\end{bmatrix}=-\bar{a}_{k-1},\quad
r_{22}=\frac{-1}{-1}\mbox{det}\begin{bmatrix}1&\bar{a}_{k-4}-\bar{a}_{k-3}\\\
-1&-\bar{a}_{k-4}\end{bmatrix}=-\bar{a}_{k-3},...,$
$r_{2j}=\frac{-1}{-1}\mbox{det}\begin{bmatrix}1&\bar{a}_{k-2j}-\bar{a}_{k-2j+1}\\\
-1&-\bar{a}_{k-2j}\end{bmatrix}=-\bar{a}_{k-2j+1},\quad\mbox{for a general
$j$}.$

By comparing these entries with the $2$nd row of (32), one see clearly that
(35) is verified.

We conclude that the entries of the first row of $R_{T_{k}}$ change values
from $1$ to $-1$ alternatingly. Thus, according to the Routh-Hurwitz stability
criterion [5, 13] and the fact that there is no zero root, the number of roots
of $\bar{p}_{k}$ which are lying in the right half-plane equals $k$.
Therefore, all the roots of $\bar{p}_{i},~{}i=1,2,\ldots,n$, lie in the right
half-plane. ∎

Remark 4. We highlight here that in the assumption of Theorem 3.2, we do not
require that $\mathcal{A}_{n}$ is symmetric. Moreover, under the same
assumption of Theorem 3.2, if $\mathcal{M}_{n}$ in (13) is the preconditioner,
then $\mathcal{M}_{n}^{-1}\mathcal{A}_{n}$ satisfies a polynomial of the form
$\prod^{n}_{i=0}\tilde{p}_{i}$, where
$\tilde{p}_{0}(\lambda)=1,~{}\tilde{p}_{1}(\lambda)=\lambda-1$ and

$\tilde{p}_{i}=\lambda\tilde{p}_{i-1}-\tilde{p}_{i-2},~{}i=2,3,4,\ldots,n.$

The proof of this conclusion is the same as that provided in [27], which is
for the symmetric cases. However, $\mathcal{M}_{n}$ is not a positive stable
preconditioner. In comparison, we show that $\mathcal{D}_{n}$ is a positive
stable preconditioner for $\mathcal{A}_{n}$ no matter it is symmetric or not.

Remark 5. We comment here that if $\mathcal{A}_{n}$ is symmetric, the proof
for verifying $\mathcal{D}_{n}$ is positive stable can be carried out similar
to that for Theorem 3.1. More clearly, we left multiply by
$\mathcal{S}_{n}^{\frac{1}{2}}$ and right multiply by
$\mathcal{S}_{n}^{-\frac{1}{2}}$ to $\mathcal{D}_{n}^{-1}\mathcal{A}_{n}$, and
denote
$\mathcal{H}_{n}=\mathcal{S}_{n}^{\frac{1}{2}}\mathcal{D}_{n}^{-1}\mathcal{A}_{n}\mathcal{S}_{n}^{-\frac{1}{2}}$.
Then, we have

$\displaystyle\mathcal{H}_{n}$
$\displaystyle=\left[\begin{array}[]{cccc}I&S_{1}^{-\frac{1}{2}}B_{1}^{T}S_{2}^{-\frac{1}{2}}&&\\\
-S_{2}^{-\frac{1}{2}}B_{1}S_{1}^{-\frac{1}{2}}&0&\ddots&\\\
&\ddots&\ddots&(-1)^{n-2}S_{n-1}^{-\frac{1}{2}}B_{n-1}^{T}S_{n}^{-\frac{1}{2}}\\\\[4.30554pt]
&&(-1)^{n-1}S_{n}^{-\frac{1}{2}}B_{n-1}S_{n-1}^{-\frac{1}{2}}&0\end{array}\right].$
(40)

We note that the Hermitian part of
$\mathcal{S}_{n}^{\frac{1}{2}}\mathcal{D}_{n}^{-1}\mathcal{A}_{n}\mathcal{S}_{n}^{-\frac{1}{2}}$
is

$\displaystyle\frac{1}{2}\left[(\mathcal{S}_{n}^{\frac{1}{2}}\mathcal{D}_{n}^{-1}\mathcal{A}_{n}\mathcal{S}_{n}^{-\frac{1}{2}})+(\mathcal{S}_{n}^{\frac{1}{2}}\mathcal{D}_{n}^{-1}\mathcal{A}_{n}\mathcal{S}_{n}^{-\frac{1}{2}})^{T}\right]$
$\displaystyle=\left[\begin{array}[]{cccc}I&0&&\\\ 0&0&\ddots&\\\
&\ddots&\ddots&0\\\\[4.30554pt] &&0&0\end{array}\right].$ (45)

Let us denote $\lambda$ as an eigenvalue of $\mathcal{H}_{n}$, $\bar{\lambda}$
as its conjugate, and ${\bf x}$ as the corresponding eigenvector. As
$\mathcal{H}_{n}$ is invertible, $\lambda$ is nonzero and $\mathbf{x}$ is not
a zero vector. We have

$\mathbf{x}^{T}\mathcal{H}_{n}\mathbf{x}=\lambda||\mathbf{x}||^{2},\quad\mathbf{x}^{T}\mathcal{H}_{n}^{T}\mathbf{x}=\bar{\lambda}||\mathbf{x}||^{2}.$

Summing the two equations together and noting from the expression of
$\mathcal{H}_{n}$, we have

$\displaystyle\mathbf{x}^{T}\left(\mathcal{H}_{n}+\mathcal{H}_{n}^{T}\right)\mathbf{x}$
$\displaystyle=2\mbox{Re}(\lambda)||\mathbf{x}||^{2},$ $\displaystyle
2||x_{1}||^{2}$ $\displaystyle=2\mbox{Re}(\lambda)||\mathbf{x}||^{2}.$ (46)

Here, $x_{1}$ is the first component of $\mathbf{x}$ and $\mbox{Re}(\lambda)$
is the real part of $\lambda$. We assert that $x_{1}$ is nonzero. This
conclusion can be proved by the method of contradiction. If $x_{1}$ is $0$,
starting with the first equation of (26), we have
$x_{2}\in\mbox{Ker}(\mathcal{B}_{1}^{T})$. However, $\mathcal{B}_{1}^{T}$ must
be of full column rank because $\mathcal{B}_{1}\mathcal{B}^{T}_{1}=I$. It
follows that $x_{2}$ is 0. Repeating this procedure as presented in (26), and
noting that $\mathcal{B}_{i}^{T}$ is of full column rank as
$\mathcal{B}_{i}\mathcal{B}^{T}_{i}=I$, we will have $x_{3}=...=x_{n}=0$. It
contradicts the fact that $\mathbf{x}$ is nonzero. Therefore, from (46), we
conclude that $\mbox{Re}(\lambda)>0$.

##  4 Additive Schur complement based preconditioners for twofold saddle point
problems

In this section, we discuss how to design preconditioners when the twofold
saddle point system is permuted into (2). Such a form of the linear system
arises naturally from domain decomposition methods [29] with two subdomains.
Without confusion, we denote

$\begin{array}[]{ll}A=\left[\begin{array}[]{cc}A_{1}&0\\\\[2.84526pt]
0&A_{3}\\\\[2.84526pt] \end{array}\right],\qquad
B^{T}=\left[\begin{array}[]{ccc}B^{T}_{1}\\\\[2.84526pt] C_{2}\\\\[2.84526pt]
\end{array}\right],\qquad
C=\left[\begin{array}[]{ccc}C_{1}&B^{T}_{2}\end{array}\right].\end{array}$
(47)

Here, $A_{3}$ can degenerate to $0$. If $A_{1}$ and $A_{3}$ are invertible,
then $A$ is invertible, and the Schur complement simply reads as

$S=A_{2}+CA^{-1}B^{T}.$

Furthermore, if $S$ is assumed to be invertible, then the following results
hold [12, 18].

######  Proposition 5.

If $A$ and $S$ are invertible and the preconditioner is

$\mathcal{Q}_{1}=\left[\begin{array}[]{cc}A&0\\\\[2.84526pt]
C&-S\\\\[2.84526pt] \end{array}\right],$

then

$\mathcal{Q}_{1}^{-1}\mathcal{A}=\left[\begin{array}[]{cc}I&A^{-1}B^{T}\\\\[2.84526pt]
0&I\\\\[2.84526pt] \end{array}\right],$

which satisfies the polynomial equation $(\lambda-1)^{2}=1$. Moreover, if
$A_{2}=0$ and the preconditioner is

$\mathcal{Q}_{D_{1}}=\left[\begin{array}[]{cc}A&0\\\\[2.84526pt]
0&S\\\\[2.84526pt] \end{array}\right],$

then

$\mathcal{T}_{Q_{D_{1}}}=\mathcal{Q}_{D_{1}}^{-1}\mathcal{A}=\left[\begin{array}[]{cc}I&A^{-1}B^{T}\\\\[2.84526pt]
S^{-1}C&0\\\\[2.84526pt] \end{array}\right],$

which satisfies the polynomial equation
$(\lambda(\lambda-1))^{2}=\lambda(\lambda-1)$. The non-zero eigenvalues of
$\mathcal{T}_{Q_{D_{1}}}$ are

$1,~{}~{}\frac{1\pm\sqrt{5}}{2}.$

The proof of this proposition can be found in [23] and [18].

Based on the formulation (47), if the preconditioner is

$\mathcal{Q}_{2}=\left[\begin{array}[]{cc}A&0\\\\[2.84526pt]
C&S\\\\[2.84526pt] \end{array}\right],$

then

$\mathcal{Q}_{2}^{-1}\mathcal{A}=\left[\begin{array}[]{cc}I&A^{-1}B^{T}\\\\[2.84526pt]
0&-I\\\\[2.84526pt] \end{array}\right],$

which satisfies the polynomial equation $(\lambda-1)(\lambda+1)=0$. Moreover,
if $A_{2}=0$ and the preconditioner is

$\mathcal{Q}_{D_{2}}=\left[\begin{array}[]{cc}A&0\\\\[2.84526pt]
0&-S\\\\[2.84526pt] \end{array}\right],$

then

$\mathcal{T}_{Q_{D_{2}}}=\mathcal{Q}_{D_{2}}^{-1}\mathcal{A}=\left[\begin{array}[]{cc}I&A^{-1}B^{T}\\\\[2.84526pt]
-S^{-1}C&0\\\\[2.84526pt] \end{array}\right],$

which satisfies the polynomial
equation$(\lambda(\lambda-1))^{2}=-\lambda(\lambda-1)$. The non-zero roots are

$1,~{}~{}\frac{1\pm\sqrt{3}i}{2}.$

We note that the polynomial equations for $\mathcal{Q}_{i}$ and
$\mathcal{Q}_{D_{i}}$ have a lower degree than $\mathcal{P}_{i}$ and
$\mathcal{Q}_{P_{i}}$ discussed in Section 2. From a purely algebraic point of
view, if the GMRES method is applied, $\mathcal{Q}_{1}$ is more favored than
$\mathcal{P}_{1}$, $\mathcal{Q}_{D_{2}}$ is more favored than
$\mathcal{P}_{D_{1}}$. However, on the other hand, it might be not easy to
directly derive a proper approximation for the additive Schur complement for
some specific problems [22, 29]. For the nested Schur complement based
preconditioners, they actually provide a recursive approach for solving block-
tridiagonal systems.

##  5 Additive Schur complement based preconditioners from the viewpoint of
domain decomposition

In this section, we generalize the discussion of Section 4 to the $n$-tuple
case. We consider a system

$\mathcal{A}_{n}=\left[\begin{array}[]{cccc}A_{1}&0&\ldots&B_{1}^{T}\\\
0&-A_{2}&\ddots&\vdots\\\ \vdots&\ddots&\ddots&B_{n-1}^{T}\\\\[4.30554pt]
C_{1}&\ldots&C_{n-1}&(-1)^{n-1}A_{n}\end{array}\right].$ (48)

Formally, the system matrix (48) arises naturally from the domain
decomposition method for an elliptic problem with $n-1$ subdomains. Like in
Section 4, one can define $A$, $C$, $B^{T}$ and $D$. The additive Schur
complement preconditioner and the corresponding theory in Section 4 can still
be applied. The resulting preconditioned system satisfies a polynomial of
degree $2$ (block triangular preconditioner) or $4$ (block-diagonal
preconditioner).

For some special forms of (48), for example, 3-by-3 block systems, it is
possible to reorder [28, 22] the interior and interface unknowns subdomain by
subdomain in a sequential manner and reduce them to the block tridiagonal
systems discussed in Section 3. If an additive Schur complement based
preconditioner is applied, then the degree of the polynomial equation that the
preconditioned system satisfied is 2 (for block triangular) or 4 (block-
diagonal). In contrast, the preconditioners based on the nested Schur
complement satisfy polynomials with the degrees being equal to $n$ (block
triangular) or $n!$ (block-diagonal). Therefore, an additive Schur complement
based preconditioner is more favored than a nested Schur complemented based
preconditioner. The Schur complement of (48), which corresponds to the
Steklov-Poincaré operator, is an additive type. Usually, it is hard to
directly invert the Schur complement. Dirichlet-Neumann, (balancing) Neumann-
Neumann, fractional Sobolev norm based preconditioners, vertex space
preconditioners, or other preconditioners can be designed (for example, for
diffusion-dominated problems). Interested readers are referred to [22, 29] for
more details. For some specific PDE problems, for example, Biot model [19, 24]
or models from fluid dynamics [6, 8, 14], it is also possible to derive an
analytic approximation of the additive Schur complement from a Fourier
analysis point of view [7, 8]. We also conducted numerical experiments for a
3-field reformulation of Biot model. In our experiments, the Schur complements
are firstly approximated by using a Fourier analysis derivation, then are
further approximated by using incomplete Cholesky factorizations. The detailed
numerical experiments are provided in Section 6.

If $n>3$, we realize that not all block tridiagonal systems can be directly
permuted into the domain decomposition form (48), although the number of
nonzero blocks of (48) and that of (10) are the same for the non-degenerate
cases. Therefore, it is also of importance to investigate the nested Schur
complement based preconditioners.

##  6 Numerical experiments based on a 3-field formulation of Biot model

A 3-field of Biot model reads as [19, 24]:

$\displaystyle-2\mu\mbox{div}\left(\bm{\epsilon}(\bm{u})\right)+\nabla\xi=\bm{f},$
(49)
$\displaystyle-\mbox{div}\bm{u}-\frac{1}{\lambda}\xi+\frac{\alpha}{\lambda}p=0,$
(50)
$\displaystyle\left(\left(c_{0}+\frac{\alpha^{2}}{\lambda}\right)p-\frac{\alpha}{\lambda}\xi\right)_{t}-K\mbox{div}\left(\nabla
p-\rho_{f}\bm{g}\right)=Q_{s}.$ (51)

Here, $\bm{u}$ is the displacement, $p$ is the fluid pressure, and
$\xi=-\alpha p+\lambda\mbox{div}\bm{u}$ is the total pressure with $\alpha$
being Biot constant, $K$ is the permeability, $\bm{f}$ is the body force,
$Q_{s}$ is the source term.

$\bm{\epsilon}({\bm{u}}):=\frac{1}{2}\left(\nabla\bm{u}+\nabla\bm{u}^{T}\right),$

and

$\lambda=\frac{E\nu}{(1+\nu)(1-2\nu)},\quad\mu=\frac{E}{2(1+\nu)}.$ (52)

The boundary conditions are

$\displaystyle\bm{u}={\bm{0}}\quad\mbox{on}~{}\Gamma_{d},$ (53)
$\displaystyle\sigma(\bm{u})\bm{n}-\alpha
p\bm{n}=\bm{h}\quad\mbox{on}~{}\Gamma_{t},$ (54) $\displaystyle
p=0\quad\mbox{on}~{}\Gamma_{p},$ (55) $\displaystyle K\left(\nabla
p-\rho_{f}\bm{g}\right)\cdot\bm{n}=g_{2}\quad\mbox{on}~{}\Gamma_{f}.$ (56)

Without loss of generality, the Dirichlet boundary conditions in (53)-(56) are
assumed to be homogeneous. The initial conditions are:

$\displaystyle\bm{u}(0)=\bm{u}_{0}\quad\mbox{and}\quad p(0)=p_{0}.$ (57)

We apply the Taylor-Hood element pair for the discretization of $\bm{u}$ and
$\xi$, and $P_{1}$ elements for $p$. The corresponding finite element spaces
are denoted as $\bm{V}^{h}$, $Q^{h}$ and $W^{h}$. We apply a backward Euler
scheme with a time stepsize $\Delta t$ for (51). After the time and spatial
discretizations, the resulted linear system is of the following form.

$Ax=\left[\begin{array}[]{ccc}A_{\bm{u}}&B^{T}_{\bm{u}\xi}&0\\\
B_{\bm{u}\xi}&-A_{\xi}&B^{T}_{\xi p}\\\ 0&B_{\xi p}&A_{p}\\\
\end{array}\right]\left[\begin{array}[]{c}\bar{\bm{u}}\\\ \bar{{\xi}}\\\
\bar{p}\\\ \end{array}\right]=b.$ (58)

Here, $\bar{\bm{u}}$, $\bar{{\xi}}$, and $\bar{p}$ are unknowns associated
with $\bm{u}$, $\xi$, and $p$ respectively;
$A_{\bm{u}},B_{{\bm{u}}\xi},A_{\xi},B_{\xi p}$, and $A_{p}$ are the matrices
resulted from the following bilinear forms.

$\displaystyle
a_{\bm{u}}(\bm{u},\bm{v})=2\mu\left(\bm{\epsilon}(\bm{u}),\bm{\epsilon}({\bm{v}})\right),\quad
b_{\bm{u}\xi}(\bm{u},\phi)=-\left(\phi,{\rm div}\bm{u}\right),$ $\displaystyle
a_{\xi}(\xi,\phi)=\frac{1}{\lambda}(\xi,\phi),\quad b_{\xi
p}(\xi,\psi)=\frac{\alpha}{\lambda}\left(\xi,\psi\right),$ $\displaystyle
a_{p}(p,\psi)=-\left(\left(c_{0}+\frac{\alpha^{2}}{\lambda}\right)p,\psi\right)-\Delta
tK\left(\nabla p,\nabla\psi\right).$

We comment here that $A_{p}$ is negative definite according to the weak form
it corresponds to. Nevertheless, the presentation in this section is
consistent with the notations in the previous sections.

Applying a Fourier analysis derivation [7, 8, 9]: firstly, we note that the
differential operator for solving (49)-(57) is

$\mathcal{A}=\left[\begin{array}[]{cccc}-2\mu\partial^{2}_{x}-\mu\partial^{2}_{y}&-\mu\partial_{y}\partial_{x}&\partial_{x}&0\\\
-\mu\partial_{y}\partial_{x}&-\mu\partial^{2}_{x}-2\mu\partial^{2}_{y}&\partial_{y}&0\\\
-\partial_{x}&-\partial_{y}&-1/\lambda&\alpha/\lambda\\\
0&0&\alpha/\lambda&-(c_{0}+\frac{\alpha^{2}}{\lambda})+K\Delta
t\Delta\end{array}\right]$

The Fourier mode of $\hat{\mathcal{A}}$ is

$\hat{\mathcal{A}}=\left[\begin{array}[]{cccc}2\mu k^{2}+\mu l^{2}&\mu
kl&ik&0\\\ \mu kl&2\mu l^{2}+\mu k^{2}&il&0\\\
-ik&-il&-\frac{1}{\lambda}&\alpha/\lambda\\\
0&0&\alpha/\lambda&-(c_{0}+\frac{\alpha^{2}}{\lambda})-K\Delta
t(k^{2}+l^{2})\end{array}\right].$

Thus, the Fourier mode for
$S_{\xi}=A_{\xi}+B_{\bm{u}\xi}A_{\bm{u}}^{-1}B_{\bm{u}\xi}^{T}$ is

$\begin{array}[]{rl}\displaystyle\hat{\mathbb{S}}_{\xi}&\displaystyle=\frac{1}{\lambda}+[-ik~{}-il]\left[\begin{array}[]{cc}2\mu
k^{2}+\mu l^{2}&\mu kl\\\ \mu kl&2\mu l^{2}+\mu k^{2}\\\
\end{array}\right]^{-1}\left[\begin{array}[]{c}ik\\\ il\\\
\end{array}\right]=\frac{1}{\lambda}+\frac{1}{2\mu}.\end{array}$

Then, the Fourier mode for $S_{p}=A_{p}+B_{\xi p}S_{\xi}^{-1}B_{\xi p}^{T}$ is

$\displaystyle\begin{array}[]{rl}\hat{\mathbb{S}}_{p}&=-\left(c_{0}+\frac{\alpha^{2}}{\lambda}\right)-K\Delta
t(k^{2}+l^{2})+\frac{\alpha}{\lambda}\left(\frac{1}{\lambda}+\frac{1}{2\mu}\right)^{-1}\frac{\alpha}{\lambda}\\\
&=-(c_{0}+\frac{\alpha^{2}}{\lambda})-K\Delta
t(k^{2}+l^{2})+\frac{2\mu\alpha^{2}}{\lambda(\lambda+2\mu)}\\\
&=-(c_{0}+\frac{\alpha^{2}}{2\mu+\lambda})-K\Delta t(k^{2}+l^{2}).\end{array}$

From the above Fourier analysis, it easy to see that $\hat{\mathbb{S}}_{\xi}$
and $\hat{\mathbb{S}}_{p}$ correspond to the following matrices.

$\mathbb{S}_{\xi}=\left(\frac{1}{\lambda}+\frac{1}{2\mu}\right)M_{\xi}~{}\quad\mbox{and}\quad~{}\mathbb{S}_{p}=A_{p}+\frac{2\mu\alpha^{2}}{\lambda(\lambda+2\mu)}M_{p}.$

Here, $M_{\xi}$ and $M_{p}$ are mass matrices which correspond to the identity
operators in $Q^{h}$ and $W^{h}$, respectively. Therefore, one can take the
above two as the approximations of $S_{\xi}$ and $S_{p}$, respectively. In
implementation, one only needs to generate the mass matrices in $Q^{h}$ and
$W^{h}$ and put the proper scalings as shown above. Correspond to $P_{1}$ to
$P_{4}$ (3) to (5), the resulted preconditioners are denoted as
$\bar{P}_{D_{1}}$ to $\bar{P}_{D_{4}}$ by replacing the exact Schur
complements by the above approximations. Similarly, correspond to $P_{D_{1}}$
to $P_{D_{4}}$ in (6) and (9), the resulted preconditioners are denoted as
$\bar{P}_{1}$ to $\bar{P}_{4}$.

Table 1: Comparisons of the different preconditioners. GMRES iteration counts. The drop tolerance for the incomplete Cholesky factorizations are set to be $1.0E-3$. | Block Diagonal Precconditioners | Block Triangular Precconditioners
---|---|---
$N$ | $\widehat{P}_{D_{1}}^{-1}$ |  $\widehat{P}_{D_{2}}^{-1}$ | $\widehat{P}_{D_{3}}^{-1}$ |  $\widehat{P}_{D_{4}}^{-1}$ | $\widehat{P}_{1}^{-1}$ |  $\widehat{P}_{2}^{-1}$ | $\widehat{P}_{3}^{-1}$ |  $\widehat{P}_{4}^{-1}$
$16\times 16$ | 55 | 56 | 50 | 67 | 27 | 49 | 49 | 36
$32\times 32$ | 108 | 108 | 96 | 127 | 55 | 97 | 97 | 72
$64\times 64$ | 229 | 229 | 204 | 273 | 118 | 206 | 206 | 157
$128\times 128$ | 470 | 470 | 449 | 593 | 261 | 434 | 434 | 345
$256\times 256$ | 1057 | 1067 | 934 | 1767 | 569 | 947 | 947 | 782
Table 2: Comparisons of the different preconditioners. GMRES iteration counts. The drop tolerance for the incomplete Cholesky factorizations are set to be $1.0E-3$. | Block Diagonal Precconditioners | Block Triangular Precconditioners
---|---|---
$N$ | $\widehat{P}_{D_{1}}^{-1}$ |  $\widehat{P}_{D_{2}}^{-1}$ | $\widehat{P}_{D_{3}}^{-1}$ |  $\widehat{P}_{D_{4}}^{-1}$ | $\widehat{P}_{1}^{-1}$ |  $\widehat{P}_{2}^{-1}$ | $\widehat{P}_{3}^{-1}$ |  $\widehat{P}_{4}^{-1}$
$16\times 16$ | 36 | 36 | 32 | 41 | 15 | 31 | 31 | 19
$32\times 32$ | 51 | 51 | 45 | 57 | 24 | 45 | 45 | 29
$64\times 64$ | 92 | 92 | 80 | 97 | 45 | 84 | 84 | 54
$128\times 128$ | 187 | 188 | 167 | 205 | 97 | 170 | 170 | 119
$256\times 256$ | 373 | 373 | 353 | 455 | 206 | 350 | 350 | 267

For practical implementation, for each of $A_{\bm{u}}$, $\mathbb{S}_{\xi}$ and
$-\mathbb{S}_{p}$ (because $\mathbb{S}_{p}$ is negative definite), we apply an
incomplete Cholesky factorization, i.e.,

$\widehat{A}=LL^{T}+R$

with a drop tolerance. Then $\widehat{A}^{-1}$ is approximated by
$L^{-T}L^{-1}$. We denote the corresponding inexact approximations of
$\bar{P}_{D_{1}}$ to $\bar{P}_{D_{4}}$ as $\widehat{P}_{D_{1}}$ to
$\widehat{P}_{D_{4}}$. Similarly, we denote the corresponding inexact
approximations of $\bar{P}_{1}$ to $\bar{P}_{4}$ as $\widehat{P}_{1}$ to
$\widehat{P}_{4}$. The implementation code, numerical results, and eigenvalues
calculations are available at https://github.com/cmchao2005/Preconditioners-
for-Biot-model.

In our numerical experiments, $\nu=0.499$, and $\lambda$ and $\mu$ are
computed by using (52). All the other physical parameters are set to be $1$.
The computational domain is a unit square. The Dirichlet type boundaries
$\Gamma_{d}$ and $\Gamma_{p}$ are the two vertical lines. Other parts of the
boundary are Neumann type. In Table 1, we report the GMRES iteration counts
for all preconditioners with the drop tolerance for the incomplete
factorizations being $10^{-3}$. In Table 2, we summarize the GMRES iteration
counts for all preconditioners with the drop tolerance for the incomplete
factorizations being $10^{-4}$. From the numerical experiments, we see that
the results based on $\widehat{P}_{D_{1}}^{-1}$ are almost the same as those
based on $\widehat{P}_{D_{2}}^{-1}$, the results based on
$\widehat{P}_{2}^{-1}$ are almost the same as those based on
$\widehat{P}_{3}^{-1}$. In comparison, $\widehat{P}_{D_{3}}$ outperforms among
the block diagonal preconditioners, and $\widehat{P}_{1}$ outperforms among
the block triangular preconditioners. Note that we have theoretically proved
that ${P}_{D_{3}}$ and $P_{1}$ are positive stable preconditioners. The
numerical experiments clearly illustrate the advantages of the positive stable
preconditioners. Although we do not rigorously prove that $\bar{P}_{D_{3}}$
and $\bar{P}_{1}$ are positive stable preconditioners in this draft, the
corresponding numerical verifications are provided in our code.

##  7 Concluding remarks

In this paper, both nested Schur complement and additive Schur complement
based preconditioners are constructed for the twofold and block tridiagonal
linear systems. The polynomial equations of the preconditioned matrices are
analyzed. It is shown that by properly selecting the sign in front of each
Schur complement, some preconditioners are positive stable. Numerical
experiments based on the Biot model are provided to show that positive stable
preconditioners outperform other preconditioners. More clearly, when inexact
elliptic approximations are incorporated in, the inexact versions of the
positive stable preconditioners are more efficient.

For block tridiagonal systems, by comparing the theoretical analysis for the
nested Schur complement based preconditioners and that for the additive type
preconditioners, our argument is that permutation is important and necessary
when designing preconditioners. These results are instructive for devising the
corresponding inexact versions of the preconditioners and iterative methods
[2].

Acknowledgments. The authors would like to express their sincere thanks to
Jarle Sogn and the reviewers for the communications and constructive
suggestions. The work of M. Cai is supported in part by NIH-BUILD grant
through UL1GM118973, NIH-RCMI grant through U54MD013376, and the National
Science Foundation awards (1700328, 1831950). The work of G. Ju is supported
in part by the National Key R & D Program of China(2017YFB1001604). The work
of J. Li is partially supported by the NSF of China No. 11971221 and the
Shenzhen Sci-Tech Fund No. RCJC20200714114556020, JCYJ20170818153840322 and
JCYJ20190809150413261, and Guangdong Provincial Key Laboratory of
Computational Science and Material Design No. 2019B030301001.

## References

  * [1] J. Adler, F. Gaspar, X. Hu, C. Rodrigo, and L. Zikatanov. Robust block preconditioners for Biot’s model. Domain decomposition methods in science and engineering XXIV, 3-16, Lect. Notes Comput. Sci. Eng., 125, Springer, Cham, 2018.
  * [2] F. Beik, and M. Benzi, Iterative methods for double saddle point systems. SIAM J. Matrix Anal. Appl., 39(2), (2018), pp. 902-921.
  * [3] F. Beik and M. Benzi, Block preconditioners for saddle point systems arising from liquid crystal directors modeling. Calcolo, 55(3) (2018), pp. 1-16.
  * [4] M. Benzi, G. Golub, and J. Liesen, Numerical solution of saddle point problems. Acta Numerica, 14, (2005), pp. 1-137.
  * [5] S. Barnettand D. Šiljak, Routh’s algorithm: a centennial survey. SIAM Review, 19(3), (1977), pp.472-489.
  * [6] M. Cai, M. Mu, and J. Xu, Preconditioning techniques for a mixed Stokes/Darcy model in porous media applications, J. Comput. Appl. Math., 233 (2009), pp. 346-355.
  * [7] M. Cai, A. Nonaka, J. Bell, B. Griffith, and A. Donev, Efficient variable-coefficient finite-volume Stokes solvers, Comm. Comput. Phys., 2014, 16(5): 1263-1297.
  * [8] M. Cai, Analysis of some projection method based preconditioners for models of incompressible flow, Appl. Numer. Math., 2015, 90: 77-90.
  * [9] M. Cai, G. Ju, and J. Li, Overlapping Schwarz preconditioners for Biot model, preprint, 2022.
  * [10] Z. Cao, Positive stable block triangular preconditioners for symmetric saddle point problems, Appl. Numer. Math., 57 (2007), pp. 899-910.
  * [11] B. Cockburn and K. Shi, Devising HDG methods for Stokes flow: an overview, Comput. Fluids., 98 (2014): pp. 221-229.
  * [12] B. Fischer, A. Ramage, D. Silvester, and A. Wathen, Minimum residual methods for augmented systems, BIT Numer. Math., 1998, 38(3): pp. 527-543.
  * [13] F. Gantmacher (J.L. Brenner translator) Applications of the Theory of Matrices, (1959), New York: Interscience.
  * [14] G. Gatica, R. Oyarzua, F. Sayas, A twofold saddle point approach for the coupling of fluid flow with nonlinear porous media flow, IMA J. Numer. Anal., 32 (2012), pp. 845-887.
  * [15] Q. Hong, J. Kraus, M. Lymbery, and F. Philo. Parameter-robust Uzawa-type iterative methods for double saddle point problems arising in Biot’s consolidation and multiple-network poroelasticity models. Math. Models Methods Appl. Sci. 30 (2020), no. 13, pp. 2523-2555.
  * [16] J. Howell and N. Walkington, Inf-sup conditions for twofold saddle point problems, Numer. Math., 118 (2011), pp. 663-693.
  * [17] N. Huang and C. Ma, Spectral analysis of the preconditioned system for the $3\times 3$ block saddle point problem. Numer. Algorithms, 2019, 81 (2): pp. 421-444.
  * [18] I. Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051.
  * [19] J. Lee, K. Mardal, and R. Winther, Parameter-robust discretization and preconditioning of Biot’s consolidation model, SIAM J. Sci. Comput., 39 (2017), pp. A1–A24.
  * [20] K. Mardal and R. Winther, Preconditioning discretizations of systems of partial differential equations. Numer. Lin. Alg. Appl., 18(1), (2011) pp.1-40.
  * [21] K. Mardal, B. Nielsen, and M. Nordaas, Robust preconditioners for PDE-constrained optimization with limited observations. BIT Numer. Math., (2017), 57(2), pp. 405-431.
  * [22] T. Mathew, Domain decomposition methods for the numerical solution of partial differential equations (Vol. 61). Springer Science & Business Media, 2008.
  * [23] M. Murphy, G. Golub, and A. Wathen, A note on preconditioning for indefinite linear systems, SIAM J. Sci. Comput., 21 (2000), pp. 1969-1972.
  * [24] R. Oyarzua and R. Ruiz-Baier, Locking-free finite element methods for poroelasticity, SIAM J. Numer. Anal., 54 (5) (2016), pp. 2951-2973.
  * [25] J. Pearson and A. Wathen, A new approximation of the Schur complement in preconditioners for PDE-constrained optimization, Numer. Linear Alg. Appl., 19 (2012), pp. 816–829.
  * [26] A. Ramage and J. Gartland, A preconditioned nullspace method for liquid crystal director modeling, SIAM J. Sci. Comput., 35 (2013), pp. B226–B247.
  * [27] J. Sogn and W. Zulehner, Schur complement preconditioners for multiple saddle point problems of block tridiagonal form with application to optimization problems, IMA J. Numer. Anal., 39 (2019), pp. 1328–1359.
  * [28] Y. Saad, Iterative methods for sparse linear systems, SIAM, 2003.
  * [29] A. Toselli, O. Widlund, Domain decomposition methods-algorithms and theory, Springer Science & Business Media, New York, 2006.
  * [30] C. Wang and P. Sun, An augmented Lagrangian Uzawa iterative method for solving double saddle-point systems with semidefinite (2,2) block and its application to DLM/FD method for elliptic interface problems, Commun. Comput. Phys. 30 (2021), no. 1, pp. 124-143.