---
author:
- |
    **Hamed Aslani$^\dag$, Davod Khojasteh Salkuyeh$^\ddag$[^1], Fatemeh Panjeh Ali Beik$^\S$\
    *[$^\dag$Faculty of Mathematical Sciences, University of Guilan, Rasht, Iran]{}*\
    *[$^\ddag$Faculty of Mathematical Sciences, and Center of Excellence for Mathematical Modelling,]{}*\
    *[Optimization and Combinational Computing (MMOCC), University of Guilan, Rasht, Iran]{}*\
    *[$^\S$Department of Mathematics, Vali-e-Asr University of Rafsanjan,P.O. Box 518, Rafsanjan, Iran]{}*\
    **
title: '**On the preconditioning of three-by-three block saddle point problems**'
---

\
[**Abstract.**]{} We establish a new iterative method for solving a class of large and sparse linear systems of equations with three-by-three block coefficient matrices having saddle point structure. Convergence properties of the proposed method are studied in details and its induced preconditioner is examined for accelerating the convergence speed of generalized minimal residual (GMRES) method. More precisely, we analyze the eigenvalue distribution of the preconditioned matrix. Numerical experiments are reported to demonstrate the effectiveness of the proposed preconditioner.\
[*Keywords*]{}: [iterative methods, sparse matrices, saddle point, convergence, preconditioning, Krylov methods. ]{}\
[*AMS Subject Classification*]{}: 65F10, 65F50, 65F08.\
\

Introduction
============

Consider the following three-by-three block system of linear equations, $$\label{eq1}
\mathcal{A} {\bf x} \equiv\left(\begin{array}{ccc}
{A} & {B^{T}} & {0} \\
{B} & {0} & {C^{T}} \\
{0} & {C} & {0}
\end{array}\right)\left(\begin{array}{l}
{x} \\
{y} \\
{z}
\end{array}\right)=\left(\begin{array}{l}
{f} \\
{g} \\
{h} 
\end{array}\right),$$ where $A\in \mathbb{R}^{n\times n}$, $B\in \mathbb{R}^{m\times n}$, $C\in \mathbb{R}^{l\times m}$, $f\in \mathbb{R}^n$, $g\in \mathbb{R}^m$ and $h\in \mathbb{R}^l$ are known, and ${\bf x}=\left(x; y; z\right)$ is an unknown vector to be determined. Here, the <span style="font-variant:small-caps;">Matlab</span> symbol $(x;y;z)$ is utilized to denote the vector $(x^{T},y^{T},z^{T})^{T}.$

In the sequel, we assume that $A$ is a symmetric positive definite matrix and the matrices $B$ and $C$ have full row rank. These assumptions guarantee the existence of a unique solution of ; see [@A2] for further details.

Evidently matrix $\cal A$ can be regarded as a $2\times 2$ block matrix using the following partitioning strategy, $$\label{part}
\mathcal{A} = \left( {\begin{array}{cc|c}
    A & {B^T } & {0}  \\
    B & 0 & C^T   \\
        \hline
    0 & C & 0  \\
    \end{array}} \right).$$ As seen, the above block matrix has a saddle point structure. Hence, we call Eq. by three-by-three block saddle point problem.

Linear system of the form arises from many practical scientific and engineering application backgrounds, e.g., the discrete finite element methods for solving time-dependent Maxwell equation with discontinuous coefficient [@A3; @A4; @A5; @A6], the least squares problems[@A7], the Karush-Kuhn-Tucker (KKT) conditions of a type of quadratic program [@A8] and so on. Since the matrices $A,$ $B$ and $C$ in are large and sparse, the solution of is suited by iterative methods. In practice, stationary iterative methods may converge too slowly or fail to converge. For this reason they are usually combined with acceleration schemes, like Krylov subspace methods [@A9]. Here, we focus on preconditioned Krylov subspace methods, especially, the preconditiond GMRES method.

As seen, the coefficient matrix $\mathcal{A}$ in Eq. can be considered in a two-by-two block form given by . The observation was used in the literature for constructing preconditioners to improve the convergence speed of Krylov subspace methods for solving , such as block triangular preconditioners [@A11; @A12; @A13], shift-splitting preconditioners [@A14] and parameterized preconditioners [@A15]. Recently, Huang and Ma [@A1] proposed the following block diagonal preconditioner, $$\label{eq992}
\mathcal{P}_{D}=\left(\begin{array}{ccc}
{A} & {0} & {0} \\
{0} & {S} & {0} \\
{0} & {0} & {C S^{-1} C^{T}}
\end{array}\right),$$ for solving in which $S=B A^{-1} B^{T}.$ They also derive all the eigenpairs of preconditioned matrix. Xie and Li [@A2] presented the following three preconditioners $${\cal P}_1=\begin{pmatrix}
A &   0  & 0 \\
B &  -S  & C^T \\
0 &   0  & CS^{-1}C^T 
\end{pmatrix},~
{\cal P}_2=\begin{pmatrix}
A &   0  & 0 \\
B &  -S  & C^T \\
0 &   0  & -CS^{-1}C^T 
\end{pmatrix},~
{\cal P}_3=\begin{pmatrix}
A &   B^T  & 0 \\
B &  -S    & 0 \\
0 &   0    & -CS^{-1}C^T 
\end{pmatrix},$$ and analyzed spectral properties of corresponding preconditioned matrices in the case $S=BA^{-1}B^T$. The reported numerical results in [@A2] show that the above preconditioners can significantly improve the convergence speed of GMRES method. It can be observed that the preconditioner ${\cal P}_1$ outperforms other preconditioners in terms of both required CPU time and number of iterations for the convergence. Here, we consider the following equivalent form of : $$\label{eq111}
{\cal B}{\bf x}\equiv \left(\begin{array}{ccc}
{A} & {B^{T}} & {0} \\
-{B} & {0} & -{C^{T}} \\
{0} & {C} & {0}
\end{array}\right)\left(\begin{array}{l}
{x} \\
{y} \\
{z}
\end{array}\right)=\left(\begin{array}{l}
{f} \\
{-g} \\
{h} 
\end{array}\right)=\bf{b}.$$ Although the coefficient matrix of the system is not symmetric, it has some desirable properties. For instance, the matrix [$\mathcal{B}$]{} is positive semidefinite, i.e., $\mathcal{B}+\mathcal{B}^{T}$ is symmetric positive semidefinite. This is a significant for the GMRES method. In fact, the restarted version of GMRES($m$) converges for all $m\geq 1$. Recently, some iterative schemes have been extended in the literature for solving . For instance, Cao [@CAOAML] presented the shift-splitting method. In [@Huang-NumerAlgor; @Huang-NLWA], the Uzawa-type methods were developed. In this work, we present a new type of iterative method for solving three-by-three block saddle point problem . Next, we extract a preconditioner from the presented iterative method and examine its performance for speeding up the convergence of GMRES.

The remainder of this paper organized as follows. Before ending this section, we present notations and basic preliminaries used in next sections. In section \[sec2\], we propose a new iterative method for solving and study its converges properties. In section \[sec3\], we extract a preconditioner from the proposed method and analyze the spectrum of preconditioned matrix. Brief discussions are given in section \[sec4\] about practical implementation of the preconditioner. In section \[sec5\], we report some numerical results and brief concluding remarks are included in section \[sec6\].

[Throughout this paper, the identity matrix is denoted by $I$. The symbol $x^{*}$ is used for the conjugate transpose of the vector $x.$ For any square matrix $A$ with real eigenvalues, the minimum and maximum eigenvalues of $A$ are indicated by $\lambda_{\min} (A)$ and $\lambda_{\max} (A)$, respectively. The notation $\rho(A)$ stands for the spectral radius of $A.$ The matrix $A\in \mathbb{R}^{n \times n}$ is called symmetric positive definite (SPD), if $A^{T}=A$ and $x^{T} A x>0$ for all nonzero $x\in \mathbb{R}^{n}.$ Similarly, the matrix $A$ is called symmetric positive semidefinite (SPSD), if $A^{T}=A$ and $x^{T}A x \geqslant 0$ for all $x\in \mathbb{R}^{n}.$ We write $A \succ 0$ $ (A \succeq 0 ),$ if $A$ is SPD (SPSD). For two given matrices $A$ and $B,$ $A\succ B$ $(A \succeq B)$ means that $A-B \succ 0$ $(A-B \succeq 0).$ The matrix $A\in \Bbb{R}^{n\times n}$ is said to be positive (semi-) definite, if $A+A^T$ symmetric positive (semi-) definite. For any matrix $W$, we shall write its null space as $null(W).$ The norm $\parallel . \parallel$ indicates the 2-norm. ]{}

The proposed iteration scheme {#sec2}
=============================

Let us first consider the following splitting for the coefficient matrix $\mathcal{B}$ in : $$\label{Split}
\mathcal{B}=\mathcal{P}-\mathcal{R},$$ where $$\mathcal{P}=\left(\begin{array}{ccc}
{A} & {B^{T}} & {0} \\
{0} & {S} & {-C^{T}} \\
{0} & {C} & {0}
\end{array}\right), \quad
\mathcal{R}=\left(\begin{array}{ccc}
{0} & {0} & {0} \\
{B} & {S} & {0} \\
{0} & {0} & {0}
\end{array}\right),$$ in which $S$ is a given symmetric positive definite matrix. It is not difficult to verify that the matrix ${\cal P}$ is nonsingular. The iteration scheme associated with splitting is given by $$\label{qe3}
{\bf x}^{(k+1)}=\mathcal{G} {\bf x}^{(k)}+{\bf c}, \quad k=0,1,2, \ldots,$$ where ${\bf x}^{(0)}$ is an initial guess, $\mathcal{G}=\mathcal{P}^{-1} \mathcal{R} $ is the iteration matrix and ${\bf c}=\mathcal{P}^{-1} {\bf b}$.

Now, we present sufficient conditions under which the iterative scheme is convergent. To this end, we first need to recall the following theorem.

\[thorn\] [@Horn1 Theorem 7.7.3] Let $A$ and $B$ be two $n\times n$ real symmetric matrices such that $A$ is positive definite and $B$ is positive semidefinite. Then $A \succeq  B$ if and only if $\rho(A^{-1}B) \le 1,$ and $A\succ B$ if and only if $\rho(A^{-1}B) < 1$.\

\[th1\] Let $A\succ 0$, $S\succ 0$ and $B$ and $C$ be full row rank matrices. If $2S \succ B A^{-1} B^{T}$ then the iterative method converges to the unique solution of for any initial guess.

Let $\lambda $ be an arbitrary eigenvalue of $\mathcal{G}=\mathcal{P}^{-1} \mathcal{R}$ with the corresponding eigenvector $w=\left(x; y; z\right)$. Consequently, we have $\mathcal{R} w=\lambda \mathcal{P}w$ which is equivalent to say that

(A x+B\^[T]{} y) =0, \[eq55\]\
( S y - C\^[T]{} z) = B x+ S y, \[eq56\]\
C y =0 .\[eq57\]

Without loss of generality, we may assume that $\lambda \neq 0$. Obviously $y \neq 0$, otherwise in view of the positive definiteness of $A$ and the assumption that $C$ has full row rank we conclude that $x$ and $z$ are both zero vectors which is in contradiction with the fact that $(x;y;z)$ is an eigenvector. From Eqs. and we can deduce that $$x=-A^{-1} B^{T} y,\quad y^{*} C^{T}= 0.$$ Multiplying both sides of Eq. on the left by $y^{*}$ and substituting the preceding equalities, we derive $$\lambda =1 - \frac{y^{*} B A^{-1} B^{T} y}{y^{*} S y}.$$ This shows that the eigenvalues of ${\cal G}$ are all real. By Theorem \[thorn\], it is immediate to conclude that $\lambda_{\max}(S^{-1}BA^{-1}B^T)=\rho(S^{-1}BA^{-1}B^T)<2$ if and only if $2S \succ B A^{-1} B^{T}$. This fact together with Courant-Fisher inequality [@A9] can deduce that $$0 < \frac{y^{*} B A^{-1} B^{T} y}{y^{*} S y} \le \lambda_{\max}(S^{-1}BA^{-1}B^T) <2.$$ Therefore, we have $$|1 - \frac{y^{*} B A^{-1} B^{T} y}{y^{*} S y}|<1,$$ which completes the proof.

We finish this section with a remark providing alternative sufficient conditions for convergence of iterative method which are stronger than $2S \succ B A^{-1} B^{T}$, however, it might be easier to checked the following sufficient conditions in some cases. To do so, we first remind the following two lemmas. The first one is a consequence of Weyl’s Theorem, see [@Horn1 Theorem 4.3.1].

\[lem2.8\] Suppose that $A$ and $B$ are two Hermitian matrices. Then, $$\begin{aligned}
        \lambda_{\max}(A+B) & \leq & \lambda_{\max}(A) + \lambda_{\max}(B),\\
        \lambda_{\min}(A+B) & \geq & \lambda_{\min}(A) + \lambda_{\min}(B).
    \end{aligned}$$

\[lem2\] [@Zhang] Suppose that $A$ is a Hermitian negative definite matrix and $B$ is Hermitian positive semidefinite. Then the eigenvalues of $AB$ are real and satisfy $${\lambda _{\min}}(A){\lambda _{\min}}(B) \le {\lambda _{\max}}(AB)
    \le {\lambda _{\max}}(A){\lambda _{\min}}(B),$$ $${\lambda _{\min}}(A){\lambda _{\max}}(B) \le {\lambda _{\min}}(AB)
    \le {\lambda _{\max}}(A){\lambda _{\max}}(B).$$

Notice that $2S \succ B A^{-1} B^{T}$ is equivalent to say that all eigenvalues of $2S - B A^{-1} B^{T}$ are positive, i.e., $\lambda_{\min}(2S - B A^{-1} B^{T}) >0$. From Lemma \[lem2.8\], it can be seen that $$\begin{aligned}
    \lambda_{\max }\left (B A^{-1} B^{T}\right)<2 \lambda_{\min}(S) , \label{qe4}
    \end{aligned}$$ implies that $2S \succ B A^{-1} B^{T}$. Using Lemma \[lem2\], one can deduce that the condition is satisfied as soon as $$\|B\|^2 < 2\lambda_{\min }(A)\lambda_{\min }(S),$$ which follows from the fact that $$\lambda_{\max }\left (B A^{-1} B^{T}\right)=    \lambda_{\max }\left ( A^{-1} B^{T}B\right) \le \lambda_{\max }\left ( A^{-1} \right)\lambda_{\max }\left (B^{T}B\right)=\frac{\|B\|^2}{\lambda_{\min}(A)}.$$

The induced preconditioner and its spectral analysis {#sec3}
====================================================

From the splitting we have $${\cal P}^{-1}{\cal B}=I-{\cal P}^{-1}{\cal R}=I-{\cal G}.$$ Therefore, under the conditions of Theorem \[th1\] the eigenvalues of ${\cal P}^{-1}{\cal B}$ are contained in the interval $(0,2]$. Thus, $$\label{eq99}
     \mathcal{P}=\left(\begin{array}{ccc}
     {A} & {B^{T}} & {0} \\
     {0} & {S} & {-C^{T}} \\
     {0} & {C} & {0}
     \end{array}\right),$$ can be used as a preconditioner to accelerate the convergence of Krylov subspace methods like GMRES for solving the system .

In the succeeding theorem, we investigate the spectral properties of ${\cal P}^{-1}{\cal B}$ in more details.

\[th5\] Let $A$ be symmetric positive definite and $B$ and $C$ be of full row rank. Then all the eigenvalues of $\mathcal{P}^{-1} \mathcal{B}$ are real and nonzero. Furthermore, $\lambda =1 $ is an eigenvalue of algebraic multiplicity at least $n+l$ and its corresponding eigenvectors are of the form $(x;-S^{-1} B x;z)$ where $x\in\Bbb{R}^{n}$ and $ z \in \Bbb{C}^l$ such that $x,z$ are not simultaneously zero.\
The remaining eigenvalues of $\mathcal{P}^{-1} \mathcal{B}$ are of the form $$\lambda=\frac{y^{*} B A^{-1} B^{T} y}{y^{*} S y},$$ and the corresponding eigenvectors are of the form $(-A^{-1} B^{T}y;y;z)$ for all $0\neq y \in null(C)$ and arbitrary $z$.

Let $\lambda$ be an arbitrary eigenvalue of $\mathcal{P} ^{-1} \mathcal{B}$ with the corresponding eigenvector $(x;y;z)$, i.e.,

A x+B\^[T]{} y=(A x + B\^[T]{} y), \[eq22\]\
-B x - C\^[T]{} z=( S y - C\^[T]{} z), \[eq23\]\
C y =Cy .\[eq24\]

Let $x=0$. If $\lambda \neq 1,$ then by we have $B^{T}y=0,$ which shows taht $y=0.$ This along with leads to $C^{T}z = 0.$ Since $C$ is a full row rank matrix, then $z=0$. Consequently, we have $(x; y; z) = (0; 0; 0)$ which contradicts with the fact that $(x;y;z)$ is an eigenvector. If $\lambda = 1$, then by and the positive definiteness of $S$ we derive that $y=0.$ In addition, the corresponding eigenvectors are $(0;0;z),$ with $z \neq 0.$ In fact, $\lambda=1$ is an eigenvalue of $\mathcal{P}^{-1} \mathcal{B}$ with multiplicity $l$ corresponding eigenvector $(0;0;z)$ with $0\neq z \in \Bbb{C}^l$.

In the following, we consider the case that $x \neq 0.$ If $y=0$, then Eqs. and are reduced to $$\begin{aligned}
            A x=\lambda A x \quad \text{and} \quad -Bx -C^{T} z=-\lambda C^{T}z, \label{eq25}
            \end{aligned}$$ respectively. The first relation shows that $ \lambda = 1.$ By substituting it into the second equality of , we have $Bx=0.$ Therefore, the corresponding eigenvectors are of the form $(x;0;z)$ with $0 \ne x \in null(B)$ and $z\in\Bbb{R}^{l}$. Notice that, in general, we can observe that $\lambda=1$ and $(x;0;z)$ is an eigenpair of $\mathcal{P}^{-1} \mathcal{B}$ where $x \in null(B)$ and $x,z$ are not simultaneously zero.

In summary, using and in view of the positive definiteness of $S$ is the, we can conclude that $\lambda = 1 $ and $(x;-S^{-1} B x;z)$ is an eigenpair of $\mathcal{P}^{-1} \mathcal{B}$.

It is immediate to see that if $x$ and $y$ are both zero vectors then $\lambda=1$ and $z$ must be a nonzero vector. In rest of the proof, we assume that $x \neq 0$ and $y \neq 0.$ If $\lambda \neq 1$, then from Eqs. and , we observe that $x=-A^{-1} B^{T}y$ and $Cy=0$, respectively. Pre-multiplying both sides of from left by $y^{*}$ and substituting deduced $x$ and $z$ into , we get $$\lambda=-\frac{y^{*} Bx}{y^{*} S y}=\frac{y^{*} B A^{-1} B^{T} y}{y^{*} S y}.$$ Hence, the corresponding eigenvectors are of the form $(-A^{-1} B^{T}y;y;z)$ for all $0\neq y \in null(C)\subseteq \Bbb{R}^{m}$ and arbitrary $z$.

Let $S$ be an arbitrary symmetric positive definite matrix. From Theorem \[th5\] we see that the non-unit eigenvalues of the preconditioned matrix $\mathcal{P}^{-1} \mathcal{B}$ satisfies $$0<\frac{\lambda_{\min}(B A^{-1} B^{T})}{\lambda_{\max}(S)}\leq \lambda=\frac{y^{*} B A^{-1} B^{T} y}{y^{*} S y} \leq \frac{\lambda_{\max}(B A^{-1} B^{T})}{\lambda_{\min}(S)}.$$

\[minpol\] Under the assumptions of Theorem \[th5\], if $S=B A^{-1} B^{T},$ then all the eigenvalues of preconditiond matrix ${\cal H}={\cal P}^{-1}{\cal B}$ are equal to $1$ and its minimal polynomial is of degree 2.

Consider the matrix $\mathcal{P}_{D}$ defined in Eq. with $S=B A^{-1} B^{T}$. Obviously, $\mathcal{P}_{D}$ is symmetric positive definite, therefore there is a symmetric positive definite matrix $\mathcal{P}_{D}^\frac{1}{2}$ such that $\mathcal{P}_{D}=\mathcal{P}_{D}^\frac{1}{2}\mathcal{P}_{D}^\frac{1}{2}$. Similar to the proof of Theorem 3.1 in [@A2], we see that the matrix ${\cal H}$ is similar to the matrix $$\begin{aligned}
\hat{{\cal H}}&:=& \mathcal{P}_{D}^\frac{1}{2} {\cal H} \mathcal{P}_{D}^{-\frac{1}{2}}\\
&=&    \mathcal{P}_{D}^\frac{1}{2} {\cal P}^{-1}{\cal B} \mathcal{P}_{D}^{-\frac{1}{2}}  \\
&=&  \left(\mathcal{P}_{D}^{-\frac{1}{2}}{\cal P} \mathcal{P}_{D}^{-\frac{1}{2}}\right)^{-1}  \left( \mathcal{P}_{D}^{-\frac{1}{2}} {\cal B} \mathcal{P}_{D}^{-\frac{1}{2}} \right) \\
&=& \left(\begin{array}{ccc}
{I} & {M^{T}} & {0} \\
{0} & {I} & {-N^{T}} \\
{0} & {N} & {0}
\end{array}\right)^{-1}
\left(\begin{array}{ccc}
{I} & {M^{T}} & {0} \\
{-M} & {0} & {-N^{T}} \\
{0} & {N} & {0}
\end{array}\right),\end{aligned}$$ where $M=S^{-\frac{1}{2}} B A^{-\frac{1}{2}}$ and $N=(CS^{-1}C^T)^{-\frac{1}{2}} C S^{-\frac{1}{2}}$. It straightforward to verify that $MM^T=I$, $NN^T=I$ and $$\hat{{\cal H}}=I+
\left(\begin{array}{ccc}
M^T(I-N^TN)M & {M^{T}(I-N^TN)} & {0} \\
{(N^TN-I)M} & {N^TN-I} & 0 \\
{NM} & {N} & {0}
\end{array}\right).$$ Direct computation reveals that $(\hat{{\cal H}}-I)^2=0$. This shows that the minimal polynomial of $\hat{{\cal H}}$, as well as ${\cal H}$ is 2.

\[remG\] Theorem \[minpol\] shows that the complete version of the GMRES method for solving the system ${\cal P}^{-1} {\cal B}{\bf x}={\cal P}^{-1}\bf{b}$ will converge in two iterations in exact arithmetic.

Implementation of the preconditioner {#sec4}
====================================

In the implementation of the preconditioner ${\cal P}$ in a Krylov subspace method like GMRES, in each iteration, a vector of the form $v=\mathcal{P}^{-1}w$ should be computed. To this end, all we need is to solve $\mathcal{P} v=w$ for $v$. If we set $v=(v_{1};v_{2};v_{3})$ and $w=(w_{1};w_{2};w_{3})$ in which $v_{1}, w_{1} \in \mathbb{R}^{n},$ $v_{2}, w_{2} \in \mathbb{R}^{m}$ and $v_{3}, w_{3} \in \mathbb{R}^{l}$, then we need to solve the system $$\left(\begin{array}{ccc}
{A} & {B^{T}} & {0} \\
{0} & {S} & {-C^{T}} \\
{0} & {C} & {0}
\end{array}\right)\left(\begin{array}{l}
{v_{1}} \\
{v_{2}} \\
{v_{3}}
\end{array}\right)=\left(\begin{array}{l}
{w_{1}} \\
{w_{2}} \\
{w_{3}} 
\end{array}\right).$$ The following algorithm is given for solving the above linear system of equations.

Algorithm 1: Computation of $(v_1;v_2;v_2)=\mathcal{P} ^{-1}(w_1;w_2;w_3)$.

\
1: Set $t_{1}=w_{3}-C S^{-1} w_{2};$\
2: Solve $(C S^{-1} C^{T})v_{3}=t_{1}$ using the Cholesky factorization of $C S^{-1} C^{T}$;\
3: Set $t_{2}=w_{2}+C^{T} v_{3};$\
4: Solve $S v_{2}= t_{2}$ ;\
5: Set $t_{3}=w_{1}-B^{T} v_{2};$\
6: Solve $A v_{1}=t_{3}$ by the Cholesky factorization of $A.$

We end this section by pointing out to the choice of SPD matrix $S$. As seen, Remark \[remG\] shows that $S=BA^{-1}B^T$ leads to an ideal case. However, by this choice, the resulting algorithm can be costly in general cases. Basically, a preconditioner is called “optimal", if the number of preconditioned iterations is independent of the size of the problem and the amount of work per iteration scales linearly with the size of the problem. Notice that for our test problems, total work (and, approximately, the corresponding CPU-time) should grow by a factor of 4 each time the value of $p$ becomes twice.

In view of Remark \[remG\] and the above discussions, we are inspired to work with $S=\diag(BA^{-1}B^T)$ which is a quite inexpensive approximation of $BA^{-1}B^T$. Our numerical results shows that the number of preconditioned iterations (associated with $\mathcal{P}$) remains fixed for all problem sizes. The proposed preconditioner, while not quite optimal, scale well with increasing size for the first test example and it is almost perfectly scalable for the second example. We also used $S=\diag(BA^{-1}B^T)$ while working with preconditioners $\mathcal{P}_D$ and $\mathcal{P}_1$ because it leads to better time scalability than the case $S=BA^{-1}B^T$; see [@A2 Tables 1 and 2].

Numerical experiments {#sec5}
=====================

In this section, we numerically solve the three-by-three saddle point problem to examine the performance of proposed preconditioner in Section \[sec3\]. In order to compare the performance of our preconditioner with the recently proposed ones in the literature, test problems are taken from [[@A1; @A2]]{}. We use the right preconditioning in conjunction with GMRES. All runs were started from the initial zero vector and terminated once the current iteration (${\bf x}^{(k)}$) satisfies $$\frac{ \| {\bf b}-\mathcal{B} {\bf x} ^{(k)} \|}{\| {\bf b}\|}< 10^{-7},$$ or the maximum number of iterations exceeds 12000. In all tests, right-hand side vector ${\bf b}$ is set ${\bf b}=\mathcal{B}e$ where $e\in\Bbb{R}^{n+m+l}$ is vector of all ones. Numerical results are presented in the tables in which “IT" and “CPU" denote the number of iterations and elapsed CPU times in second, respectively. To show the accuracy of the methods we report the values $$Err=\frac{\|x^{(k)}- x^{*} \|} {\| x^{*} \|},$$ in the tables where $x^{*}$ stands for the exact solution of . All runs were performed in MATLAB R2019a on a PC with intel (R) Core(TM) i5-8265U CPU @ 1.60 GHz 8.GB.

[[@A1; @A2]]{} \[ex1\] Consider the saddle point problem with $$\begin{aligned}
        A=\left(\begin{array}{cc}
        {I \otimes T+T \otimes I} & {0} \\
        {0} & {I \otimes T+T \otimes I}
        \end{array}\right) \in \mathbb{R}^{2 p^{2} \times 2 p^{2}},
        \end{aligned}$$ $
        B=(I \otimes F \quad F \otimes I) \in \mathbb{R}^{p^{2} \times 2 p^{2}}$ and $C=E \otimes F \in \mathbb{R}^{p^{2} \times p^{2}}
        $ where $$\begin{aligned}
        T=\frac{1}{h^{2}} \cdot \operatorname{tridiag}(-1,2,-1) \in \mathbb{R}^{p \times p}, \quad F=\frac{1}{h} \cdot \operatorname{tridiag}(0,1,-1) \in \mathbb{R}^{p \times p},
        \end{aligned}$$ and $
        E=\operatorname{diag}\left(1, p+1, \cdots, p^{2}-p+1\right)
      $ in which $\otimes$ denotes Kronecker product and $h=\frac{1}{p+1}$ the discretization meshsize.\

           $p$                              
  -- ----- ---------- ---------- ---------- ----------
           16         32         48         64
     IT    870        3100       6759       11181
     CPU   0.74       85.95      884.22     3409.13
     Err   3.02e-08   8.86e-08   2.30e-05   9.62e-08
     IT    35         38         39         40
     CPU   0.07       0.28       0.93       3.35
     Err   1.06e-6    4.00e-06   5.07e-06   1.68e-05
     IT    26         29         31         32
     CPU   0.02       0.23       0.88       3.30
     Err   1.97e-07   9.11e-07   3.50e-07   6.07e-07
     IT    2          2          2          2
     CPU   0.02       0.21       0.83       3.19
     Err   2.33e-13   1.40e-12   4.72e-12   1.18e-11

  : Numerical results for GMRES to solve Example \[ex1\]

\[tab1\]

In this example, we take $S=\diag(B A^{-1}B^{T}).$ Table \[tab1\] shows the iteration counts and the CPU time of the GMRES with the preconditioners $\mathcal{P}_{D},$ ${\mathcal{P}_1}$ and $\mathcal{P}.$ To see the effectiveness of preconditioners, we have also reported the numerical results of the GMRES method without preconditioning. Numerical results illustrate that the preconditioners can drastically reduce the number of iterations and CPU time that the GMRES method requires for the convergence. As seen, ${\cal P}$ is superior to the other examined preconditioners. An interesting observation which can be posed here is that the GMRES with the preconditioner ${\cal P}$ for all the test problems converges in two iterations and the computed solution has the best accuracy among the preconditioners. Fig. \[fg1\] plots the eigenvalues of matrices ${\cal B}$ , $\mathcal P_{D}^{-1} \mathcal{B},$ ${\mathcal P}_1^{-1} \mathcal{B}$ and $\mathcal P^{-1} \mathcal{B}$ for $p=16.$ It is seen that the eigenvalues of $\mathcal P^{-1} \mathcal{B}$ are clustered more than the others.

![[Eigenvalue distributions of $\mathcal{B},$ $\mathcal{P}_{D}^{-1} \mathcal{B},$ ${\mathcal{P}}^{-1} \mathcal{B}$ and ${\mathcal{P}}^{-1} \mathcal{B}$ (from the left to right) with $p=16$ for Example \[ex1\]]{}. \[fg1\] ](edpm52n "fig:"){height="3.5cm" width="3.5cm"} ![[Eigenvalue distributions of $\mathcal{B},$ $\mathcal{P}_{D}^{-1} \mathcal{B},$ ${\mathcal{P}}^{-1} \mathcal{B}$ and ${\mathcal{P}}^{-1} \mathcal{B}$ (from the left to right) with $p=16$ for Example \[ex1\]]{}. \[fg1\] ](edpm52d "fig:"){height="3.5cm" width="3.5cm"} ![[Eigenvalue distributions of $\mathcal{B},$ $\mathcal{P}_{D}^{-1} \mathcal{B},$ ${\mathcal{P}}^{-1} \mathcal{B}$ and ${\mathcal{P}}^{-1} \mathcal{B}$ (from the left to right) with $p=16$ for Example \[ex1\]]{}. \[fg1\] ](edpm52l "fig:"){height="3.5cm" width="3.5cm"} ![[Eigenvalue distributions of $\mathcal{B},$ $\mathcal{P}_{D}^{-1} \mathcal{B},$ ${\mathcal{P}}^{-1} \mathcal{B}$ and ${\mathcal{P}}^{-1} \mathcal{B}$ (from the left to right) with $p=16$ for Example \[ex1\]]{}. \[fg1\] ](edpm52b "fig:"){height="3.5cm" width="3.5cm"}

[[@A1; @A2]]{}\[xe1\] Consider the three-by-three block saddle point problem for which $$\begin{aligned}
        A=\operatorname{diag}\left(2 W^{T} W+D_{1}, D_{2}, D_{3}\right) \in \mathbb{R}^{n \times n},
        \end{aligned}$$ is a block-diagonal matrix, $$\begin{aligned}
        B=\left[E,-I_{2 \widetilde{p}}, I_{2 \widetilde{p}}\right] \in \mathbb{R}^{m \times n} \text { and } \quad C=E^{T} \in \mathbb{R}^{l \times m},
        \end{aligned}$$ are both full row-rank matrices where $\tilde{p}=p^{2},$ $\hat{p}=p(p+1);$ $W=(w_{i j}) \in \mathbb{R}^{\hat{p} \times \hat{p}}$ with $w_{i j}=e^{-2\left((i / 3)^{2}+(j / 3)^{2}\right)}$; $D_{1}=I_{\widehat{p}}$ is an identity matrix; $D_{i}=\operatorname{diag}(d_{j}^{(i)}) \in \mathbb{R}^{2\tilde{p} \times 2\tilde{p}},$ $i=2,3,$ are diagonal matrices, with $$\begin{aligned}
        d_{j}^{(2)}=\left\{\begin{array}{ll}
        {1,} & {\text { for } \quad 1 \leq j \leq \tilde{p}}, \\
        {10^{-5}(j-\tilde{p})^{2}},  & \text { for }   {\tilde{p}+1 \leq j \leq 2 \tilde{p}},
        \end{array}\right.
        \end{aligned}$$ $$\begin{aligned}
        d_{j}^{(3)}=10^{-5}(j+\tilde{p})^{2}
        \text { for } 1 \leq j \leq 2 \tilde{p},
        \end{aligned}$$ and $$\begin{aligned}
        E=\left(\begin{array}{cc}
        {\widehat{E} \otimes I_{p}} \\
        {I_{p} \otimes \widehat{E}}
        \end{array}\right), \quad \widehat{E}=\left(\begin{array}{ccccc}
        {2} & {-1} & {} & {} & {}\\
        {} & {2} & {-1} & {} & {} \\
        {} & {} & {\ddots} & {\ddots} & {}\\
        {} & {} & {} & {2} & {-1}
        \end{array}\right) \in \mathbb{R}^{p \times(p+1)}.
        \end{aligned}$$ Similar to Example \[ex1\], we take $S=diag(B A^{-1}B^{T}).$ The numerical results are listed in Table \[tab2\] for different values of $p.$ As seen, the proposed preconditioner surpasses the others in terms of the iteration counts, CPU time and the accuracy of computed solution.

Fig. \[fig2\] displays the eigenvalue distribution of the original coefficient matrix, $\mathcal{P}_{D}^{-1} \mathcal{B},$ ${\cal P}_1^{-1} \mathcal{B}$ and $\mathcal{P}^{-1} \mathcal{B}$ for $p=16.$ As observed, eigenvalues of $\mathcal{P}^{-1} \mathcal{B}$ are more clustered around the point $(1,0)$ than the others.

           $p$                              
  -- ----- ---------- ---------- ---------- ----------
           16         32         48         64
     IT    261        557        1180       1815
     CPU   0.24       4.52       41.95      145.66
     Err   9.96e-08   9.90-08    9.99e-08   9.99e-08
     IT    72         111        98         80
     CPU   0.08       0.37       0.44       0.51
     Err   3.37e-07   7.16e-06   6.47e-05   1.88e-04
     IT    38         61         54         45
     CPU   0.02       0.12       0.17       0.24
     Err   3.46e-07   7.30e-08   6.15e-5    1.82e-04
     IT    15         23         20         17
     CPU   0.02       0.07       0.09       0.14
     Err   5.42e-08   7.30e-08   9.84e-08   7.34e-08

  : Numerical results for GMRES to solve Example \[xe1\]

\[tab2\]

![[Eigenvalue distributions of $\mathcal{B},$ ${\mathcal{P}_{D}}^{-1} \mathcal{B}, $ ${\mathcal{P}}_1^{-1} \mathcal{B}$ and ${\mathcal{P}}^{-1} \mathcal{B}$ with $p=16$ for Example \[xe1\]]{}. \[fig2\]](edpm53n "fig:"){height="3.5cm" width="3.5cm"} ![[Eigenvalue distributions of $\mathcal{B},$ ${\mathcal{P}_{D}}^{-1} \mathcal{B}, $ ${\mathcal{P}}_1^{-1} \mathcal{B}$ and ${\mathcal{P}}^{-1} \mathcal{B}$ with $p=16$ for Example \[xe1\]]{}. \[fig2\]](edpm53d "fig:"){height="3.5cm" width="3.5cm"} ![[Eigenvalue distributions of $\mathcal{B},$ ${\mathcal{P}_{D}}^{-1} \mathcal{B}, $ ${\mathcal{P}}_1^{-1} \mathcal{B}$ and ${\mathcal{P}}^{-1} \mathcal{B}$ with $p=16$ for Example \[xe1\]]{}. \[fig2\]](edpm53l "fig:"){height="3.5cm" width="3.5cm"} ![[Eigenvalue distributions of $\mathcal{B},$ ${\mathcal{P}_{D}}^{-1} \mathcal{B}, $ ${\mathcal{P}}_1^{-1} \mathcal{B}$ and ${\mathcal{P}}^{-1} \mathcal{B}$ with $p=16$ for Example \[xe1\]]{}. \[fig2\]](e3dpm53b "fig:"){height="3.5cm" width="3.5cm"}

Conclusions {#sec6}
===========

A new stationary iterative method was constructed for solving a class of three-by-three block saddle point problems. We analyzed the convergence properties of the elaborated stationary method. We further examined the performance of induced preconditioner from the proposed method. More precisely, the eigenvalue distribution of the preconditioned matrix was studied. Our numerical tests illustrated that the proposed preconditioner is more effective than the other tested preconditioners in the literature.

[99]{}

F. Assous, P. Degond, E. Heintze, P.A. Raviart, J. Segre, On a finite-element method for solving the three-dimensional Maxwell equations, J. Comput. Phys. 109 (1993) 222–237.

M. Benzi, M.K. Ng, Q. Niu, Z. Wang, A relaxed dimensional factorization preconditioner for the incompressible Navier-Stokes equations, J. Comput. Phys. 230 (2011) 6185–6202.

Y. Cao, Shift-splitting preconditioners for a class of block three-by-three saddle point problems, Appl. Math. Lett. 96 (2019) 40–46.

Z.-H. Cao, Positive stable block triangular preconditioners for symmetric saddle point problems, Appl. Numer. Math. 57 (2007) 899–910.

Z.-M. Chen, Q. Du, J. Zou, Finite element methods with matching and nonmatching meshes for Maxwell equations with discontinuous coefficients, SIAM J. Numer Anal. 37 (2000) 1542–1570.

C.-R. Chen, C.-F. Ma, A generalized shift-splitting preconditioner for singular saddle point problems, Appl. Math. Comput. 269 (2015) 947–955.

Z.-M. Chen, Q. Du, J. Zou, Finite element methods with matching and nonmatching meshes for Maxwell equations with discontinuous coefficients, SIAM J. Numer. Anal. 37 (2000) 1542–1570.

P. Ciarlet, J. Zou, Finite element convergence for the Darwin model to Maxwell’s equations, RAIRO Math. Modelling Numer. Anal. 31 (1997) 213–249.

H.C. Elman, D.J. Silvester, A.J. Wathen, Performance and analysis of saddle point preconditioners for the discrete steady-state Navier-Stokes equations, Numer. Math. 90 (2002) 665–688.

D.R. Han, X.M. Yuan, Local linear convergence of the alternating direction method of multipliers for quadratic programs, SIAM J. Numer. Anal. 51 (2013) 3446–3457.

R.A. Horn, C.R. Johnson, Matrix Analysis, Cambridge University Press, Cambridge, UK, 1985.

N. Huang, C.-F. Ma, Spectral analysis of the preconditioned system for the 3 $\times$ 3 block saddle point problem, Numer. Algorithms, 81 (2019) 421–444.

N. Huang, Variable parameter Uzawa method for solving a class of block three-by-three saddle point problems, Numerical Algorithms, 2020, <https://doi.org/10.1007/s11075-019-00863-y>.

N. Huang, Y.‐H. Dai, Q. Hu, Uzawa methods for a class of block three‐by‐three saddle‐point problems, Numer Linear Algebra Appl. 26 (2019) e2265.

Y.-F. Ke, C.-F. Ma, The parameterized preconditioner for the generalized saddle point problems from the incompressible Navier–Stokes equations, J. Comput. Appl. Math. 37 (2018) 3385–3398.

Y. Saad, Iterative Methods for Sparse Linear Systems, Second Edition, Society for Industrial and Applied Mathematics, Philadelphia, 2003.

X. Xie, H.-B. Li, A note on preconditioning for the $3 \times 3$ block saddle point problem, Comput. Math. Appl. 2020, <https://doi.org/10.1016/j.camwa.2020.01.022>.

J.-Y. Yuan, Numerical methods for generalized least squares problems, J. Comput. Appl. Math. 66 (1996) 571–584.

F. Zhang, Q. Zhang, Eigenvalue inequalities for matrix product, IEEE Trans. Automat. Control, 51 (2006), 1506–1509.

[^1]: Corresponding author. Emails: hamedaslani525@gmail.com (H. Aslani), khojasteh@guilan.ac.ir (D.K. Salkuyeh), f.beik@vru.ac.ir (F.P.A. Beik)
