\documentclass{mc}
%%%%% journal  info -- DO NOT CHANGE %%%%%%%%%
\setcounter{page}{177}
\renewcommand\thisnumber{x}
\renewcommand\thisyear {2021}
\renewcommand\thismonth{xxx}
\renewcommand\thisvolume{26}
\renewcommand\datereceived{September 22, 2020}
\renewcommand\dateaccepted{January 8, 2021}


% \usepackage{color}
% \usepackage{graphicx,epsfig}

\newcommand{\cred}[1]{{\color{blue}  #1}}
\def\h{\hspace{-0.0cm}}
\def\A{{\mathcal A}}
\def\I{{\mathcal I}}
\def\x{{\bf x}}
\def\b{{\bf b}}
\def\n{{\bf n}}
\def\t{\theta}
\def\vec{\textrm{vec}}
\def\ran{\textrm{rand}}
\def\sp{\textrm{span}}
\def\diag{\textrm{diag}}
\def\trd{\textrm{tridiag}}
\def\pent{\textrm{pentadiag}}
\def\tr{\textrm{tr}}
\def\cond{\textrm{cond}}

\usepackage[displaymath,mathlines,pagewise]{lineno}

\begin{document}
%% \begin{linenumbers}

\markboth{D.\,K.\,Salkuyeh, H.\,Aslani and Z.\,-Z.\,Liang}{An APSS preconditioner for the three-by-three block saddle point problems}

\title[How to use the {\sf mc.cls} class file]{An alternating positive semidefinite splitting preconditioner for the three-by-three block saddle point problems\thanks{The work of Davod Khojasteh Salkuyeh is partially supported by the University of Guilan and the work of Zhao-Zheng Liang is supported by the National Natural Science Foundation of China (No 11801242).}}

\author[Davod Khojasteh Salkuyeh, Hamed Aslani and Zhao-Zheng Liang]{Davod Khojasteh Salkuyeh\affil{1,2}\comma\corrauth, Hamed Aslani\affil{3} and Zhao-Zheng Liang\affil{4}}

\address{\affilnum{1}\ Faculty of Mathematical Sciences, University of Guilan, Rasht, Iran\\
\affilnum{2}\  Center of Excellence for Mathematical Modelling, Optimization and Combinational Computing (MMOCC), University of Guilan, Rasht, Iran\\
\affilnum{3}\ Faculty of Mathematical Sciences, University of Guilan, Rasht, Iran\\
\affilnum{4}\ School of Mathematics and Statistics, Lanzhou University, Lanzhou, P.\,R.\,China
}

\emails{{\tt khojasteh@guialn.ac.ir}\,\,(D.\,K.\,Salkuyeh), {\tt hamedaslani525@gmail.com}\,\,(H.\,Aslani), {\tt liangzz@lzu.edu.cn}\,\,(Z.\,-Z.\,Liang)}


\begin{abstract}
 Using the idea of the dimensional splitting method, we present an iteration method for solving three-by-three block saddle
 point problems which  appear in linear programming and finite element discretization of the Maxwell equation. We prove that  the method is unconditionally convergent. Then the induced preconditioner is used to accelerate the convergence of the GMRES method for solving the system. Numerical results are presented to compare the performance of the method with some existing ones.
\end{abstract}

\keywords{saddle point, block, dimensional, split, preconditioner, GMRES}

\ams{65F10, 65F50, 65F08}

\maketitle

\section{Introduction} \label{Sec1}
We consider the block three-by-three saddle point problems of the form
\begin{equation}\label{MSys}
{\cal K}\x\equiv	
\begin{pmatrix}
A & B^T & 0 \\
B &  0  & C^T \\
0 &  C  & 0
\end{pmatrix}
\begin{pmatrix}
x \\
y \\
z
\end{pmatrix}
=
\begin{pmatrix}
f \\
g \\
h
\end{pmatrix},
\end{equation}
where $A\in \Bbb{R}^{n\times n}$, $B\in \Bbb{R}^{m\times n}$ and $C\in \Bbb{R}^{l\times m}$. Here, $f\in \Bbb{R}^{n}$, $g\in \Bbb{R}^{m}$ and $h\in \Bbb{R}^{l}$.
In this case, the coefficient matrix of system \eqref{MSys} is of order ${\bf n}\times {\bf n}$, where ${\bf n} = n + m + l$. Systems of the form \eqref{MSys} appear in many applications such as finite element methods for solving time-dependent Maxwell equations with discontinuous coefficients in general three-dimensional Lipschitz polyhedral domains \cite{ChenSIAM}, solving quadratic programs \cite{HanSIAM} and the Picard iteration method for a class of mixed finite element scheme for stationary magnetohydrodynamics models \cite{Hu-Arxiv}.

System (\ref{MSys}) can be equivalently rewritten in the following form:
\begin{equation}\label{TSys}
\mathcal{A} {\bf x}\equiv
\begin{pmatrix}
A & B^T & 0 \\
-B &  0  & -C^T \\
0 &  C  & 0
\end{pmatrix}
\begin{pmatrix}
x \\
y \\
z
\end{pmatrix}
=
\begin{pmatrix}
f \\
-g \\
h
\end{pmatrix}={\bf b}.
\end{equation}
Unlike the coefficient matrix of system \eqref{MSys}, the matrix $\A$ is nonsymmetric; however,
since the matrix $\A+\A^T$ is symmetric positive semidefinite, we see that the matrix $\A$ has the following good
properties:\\
1. $\A$ is semipositive real, that is, $v^TAv \geq 0$, for all $v\in\Bbb{R}^{\n}$;\\
2. $\A$ is positive semistable, which means that $\Re(\lambda)\geq 0$ for all $\lambda\in \sigma(\A)$, where $\sigma(\A)$ denotes the spectrum of $\A$. \\
The above results are significant for Krylov subspace methods like GMRES (see \cite[p. 420]{BenziJCM} and   \cite[Theorem 6.30]{SaadBook}).

Recently, some iteration methods and preconditioning techniques have been presented for system \eqref{TSys}.
In \cite{HuangMa}, Huang and Ma analyzed the spectral properties of the preconditioned matrix  ${\cal P}_D^{-1}\A$, where the preconditioner ${\cal P}_D$ is of the form
\begin{equation}\label{PDPrecon}
{\cal P}_D=\begin{pmatrix}
A &  0  & 0 \\
0 &  S  & 0 \\
0 &  0   & CS^{-1}C^T
\end{pmatrix},
\end{equation}
with $S=BA^{-1}B^T$. They proved that when $A$ is symmetric positive definite (SPD) and $C$ and $B$ are of full rank, then all the matrices $A$, $S$ and $CS^{-1}C^T$ are SPD and the preconditioned matrix ${\cal P}_D^{-1}\A$ has six different eigenvalues. Therefore, the GMRES method for the system ${\cal P}_D^{-1}\A\x={\cal P}_D^{-1}\b$ will converge in at most six iterations.
In the implementation of the preconditioner ${\cal P}_D$ within a Krylov subspace iteration method, at each iteration, a system of linear equations of the form ${\cal P}_Dw=r$ should be solved. To do so, three systems with coefficient matrices $A$, $S$ and $CS^{-1}C^T$ should be solved.  Hence, these systems can be solved directly using the Cholesky factorization, or inexactly using the conjugate gradient (CG) method. Obviously, forming the matrices $S$ and $CS^{-1}C^T$ is too expensive. Moreover, these matrices would be dense even if the matrices $A$, $B$ and $C$ are sparse.
On the other hand, in each matrix-vector multiplication of solving the system with the matrix $S$ using the CG method a system with the coefficient matrix $A$ should be solved,  which is very costly. Furthermore, this would be even more complicated for the matrix $CS^{-1}C^T$.
Here, we should also note that if the matrix $A$ is ill-conditioned, then the matrix $S$ as well as $CS^{-1}C^T$ are often ill-conditioned, too. To overcome this problem, in \cite{HuangMa}, the authors proposed to use $\hat{A}$ and $\hat{S}$ as approximations of $A$ and $S$, respectively. However, they solved the obtained subsystems directly.  In \cite{Xie-Li-CAMWA}, Xie and Li presented the following three preconditioners:
\begin{eqnarray}\label{P123}
\begin{split}
{\cal P}_1&=\begin{pmatrix}
A &   0  & 0 \\
B &  -S  & C^T \\
0 &   0  & -CS^{-1}C^T
\end{pmatrix},\quad
{\cal P}_2=\begin{pmatrix}
A &   0  & 0 \\
B &  -S  & C^T \\
0 &   0  & CS^{-1}C^T
\end{pmatrix},\\[1ex]
{\cal P}_3&=\begin{pmatrix}
A &   B^T  & 0 \\
B &  -S    & 0 \\
0 &   0    & -CS^{-1}C^T
\end{pmatrix}.
\end{split}
\end{eqnarray}
It was proved that (see \cite{Xie-Li-CAMWA})
\[
{\cal T}_1=\sigma({\cal P}_1^{-1}{\cal K})=\{-1,1\},\quad {\cal T}_2=\sigma({\cal P}_2^{-1}{\cal K})=\{1\}, \quad {\cal T}_3=\sigma({\cal P}_3^{-1}{\cal K})=\{\pm\frac{1}{2},1\},
\]
and that under the assumptions that $A$ is SPD and $C$ and $B$ are of full rank,  all the preconditioned matrices  ${\cal T}_i\, (i=1,2,3)$ have a minimal polynomial of degree 3. Here, $\sigma(.)$ denotes the spectrum of the matrix. Therefore, the GMRES method for solving the corresponding preconditioned systems converges in at most three iterations. Obviously, the implementation of the preconditioners is quite similar to the preconditioner ${\cal P}_D$, and the one with ${\cal P}_3$ is even more complicated. In \cite{CaoSS}, Cao presented the shift-splitting (SS) iteration method for solving system \eqref{TSys}. Next in \cite{WangSS}, Wang and Zhang generalized the SS iteration method to solve   system \eqref{TSys}. It has been proved that the SS iteration method and its generalized version are unconditionally convergent (see \cite{CaoSS,WangSS}). Huang \cite{Huang-NumerAlgor} and Huang et al. \cite{Huang-NLWA} presented the Uzawa-type methods for solving system \eqref{MSys}.

In \cite{DS-Benzi}, Benzi and Guo proposed the dimensional splitting method and the corresponding dimensional split preconditioner for saddle point problems. Then in \cite{DS-Ma}, Ke and Ma corrected an error in \cite{DS-Benzi}. Next, a modified version of the dimensional splitting method was presented by Cao et al. in \cite{DS-Cao} for solving  generalized saddle point problems.
When the (2,2) block of the matrix $\mathcal{A}$ is symmetric positive definite, using the idea of the dimensional splitting method \cite{DS-Benzi}, Liang and Zhang \cite{Liang} proposed  alternating positive semidefinite splitting (APSS) for solving the problem. However, in our problem, the (2,2) block of the matrix $\mathcal{A}$ is zero. In this paper,  we apply the APSS iteration method for solving system \eqref{TSys} and investigate its convergence properties. We prove that the method is unconditionally convergent. Then we apply the induced preconditioner to accelerate the convergence of the GMRES method for solving the system.

The rest of the paper is organized as follows. In  Section \ref{Sec2}, we give conditions for the nonsingularity of the matrix $\A$. Section \ref{Sec3} is devoted to the  application of the APSS method to solving system \eqref{TSys}. Convergence of the method is presented in Section \ref{Sec4}. We present some implementation issues for the APSS preconditioner in Section \ref{Sec5}. Numerical results are presented in Section \ref{Sec6}. Some concluding remarks are presented in Section \ref{Sec7}.

\section{Existence of the solution}\label{Sec2}

We state the following theorem for the existence of a solution to system \eqref{TSys}.
\begin{theorem}\label{ExistThm}
	Assume that $C$ has full row rank and $A$ is symmetric positive semidefinite. If ${\cal N}(A) \cap {\cal N}(B)=\{0\}$,
	${\cal N}(B^T) \cap {\cal N}(C)=\{0\}$ and
	${\cal R}(B) \cap {\cal R}(C^T)=\{0\}$, then the matrix ${\cal A}$ defined in Eq. \eqref{TSys} is nonsingular.
	\end{theorem}
\begin{proof}
		Let $\x=(u;v;p)$ and $\A\x=0$. Therefore,
		\begin{eqnarray}
		Au+B^T v&\h=\h& 0, \label{Exist01}\\
		Bu+C^Tp &\h=\h& 0, \label{Exist02} \\
		Cv &\h=\h& 0.  \label{Exist03}
		\end{eqnarray}	
		Pre-multiplying Eq. \eqref{Exist01} by $u^T$  and by using Eq. \eqref{Exist02} and Eq. \eqref{Exist03} it follows
		\begin{eqnarray*}
			0 &\h=\h& u^TAu+u^TB^T v=u^TAu+(Bu)^T v\\
			&\h=\h& u^TAu-p^TCv\\ %\qquad  \text{(by using Eq. \eqref{Exist02})}\\
			&\h=\h& u^TAu. %\qquad \text{(by using Eq. \eqref{Exist03})}
		\end{eqnarray*}
		Now, since the matrix $A$ is symmetric positive semidefinite, we deduce that $Au=0$. Therefore, from Eq. \eqref{Exist01} we get $B^Tv=0$.
		This along with Eq. \eqref{Exist03} shows that $v\in {\cal N}(B^T) \cap {\cal N}(C)$. Thus, $v=0$. On the other hand, from Eq. \eqref{Exist02} we have  $Bu=-C^Tp \in {\cal R}(B) \cap {\cal R}(C^T)$. Hence, $Bu=0$ and $C^Tp=0$. So, $p=0$, since $C$ has full row rank. Finally, from $Bu=0$ and $Au=0$, we conclude that $u\in {\cal N}(A) \cap {\cal N}(B)$, which reveals that $u=0$.  Hereupon, $\x=(u;v;p)=0$, which completes the proof.
	\end{proof}


It is easy to prove  that if ${\cal R}(B)\subseteq {\cal N}(C)$, then ${\cal R}(B)\cap {\cal R}(C^T)=\{0\}$. So in Theorem \ref{ExistThm}, one may replace the condition ${\cal R}(B)\cap {\cal R}(C^T)=\{0\}$ by  ${\cal R}(B)\subseteq {\cal N}(C)$.   The following corollaries can be easily deduced  using Theorem \ref{ExistThm} and its proof.
\begin{corollary}[see {\cite{Huang-NumerAlgor,Huang-NLWA,Xie-Li-CAMWA}}]
	Assume that $A$ is symmetric positive definite  and the matrices  $B$  and $C$ are of  full row rank. Then the matrix ${\cal A}$ defined in Eq. \eqref{TSys} is nonsingular.
\end{corollary}

\begin{corollary}[see {\cite{Xie-Li-CAMWA}}]\label{Cor2}$\,$
\begin{enumerate}
	\item[$($i$)$]	Assume that $A$ is a symmetric positive semidefinite and the matrices  $B$  and $C$ are of  full row rank. Then the matrix ${\cal A}$ defined in Eq. \eqref{TSys} is nonsingular if and only if ${\cal N}(A) \cap {\cal N}(B)=\{0\}$.
\item[$($ii$)$]	Assume that $A$ is a symmetric positive definite matrix and the matrix $C$ has full row rank. Then the matrix ${\cal A}$ defined in Eq. \eqref{TSys} is nonsingular if and only if ${\cal N}(B^T) \cap {\cal N}(C)=\{0\}$.
\end{enumerate}
\end{corollary}


\section{The APSS method}\label{Sec3}

We split the matrix $\A$ into $\A = \A_1 + \A_2$, where
\begin{equation*}\label{MSys0}
\A_1=
\begin{pmatrix}
A & B^T & 0 \\
-B &  0  & 0 \\
0 &  0  & 0
\end{pmatrix},
\quad
\A_2=
\begin{pmatrix}
0 & 0 & 0 \\
0 &  0  & -C^T \\
0 &  C  & 0
\end{pmatrix}.
\end{equation*}
Both matrices $\A_1$ and $\A_2$ are positive semidefinite. So, using the splittings
\[
\A= (\alpha \I + {\A}_1)-(\alpha \I - \A_2) = (\alpha \I + \A_2)  - (\alpha  \I - \A_1),
\]
where $\alpha>0$ is a given real parameter and  $\I$ is the identity matrix of order $\n$, we state the APSS iteration
\begin{equation}\label{DS-Met}
\begin{cases}
(\alpha \I+\A_1){\x^{(k+\frac{1}{2})}}=(\alpha \I-\A_2){\x^{(k)}}+\b, \\
(\alpha \I+\A_2){\x^{(k+1)}}=(\alpha \I-\A_1){\x^{(k+\frac{1}{2})}}+\b, \\
\end{cases} k=0,1,2,\ldots,
\end{equation}
where $\x^{(0)}\in\Bbb{R}^{\n}$ is a given parameter.
It follows from Lemma 1 in \cite{SD-Benzi-SIMAX} that both matrices $\alpha \I+\A_1$ and $\alpha \I+\A_2$ are nonsingular.
Computing ${\bf x}^{(k+\frac{1}{2})}$ from the first equation in \eqref{DS-Met} and substituting  in the second one results in the stationary iteration
\begin{equation}\label{DS-Stat}
{\bf x}^{(k+1)}={\cal T}_{\alpha} {\bf x}^{(k)}+{\bf c},
\end{equation}
where
\[
{\cal T}_{\alpha}=\left( \alpha \I + \A_2 \right)^{-1}\left( \alpha \I - \A_1 \right)
\left( \alpha \I + \A_1 \right)^{-1}\left( \alpha \I -\A_2 \right),
\]
and
\[
{\bf c}=2 \alpha\left( \alpha \mathcal I + \A_2 \right)^{-1} \left( \alpha \mathcal I + \A_1 \right)^{-1} {\bf b}.
\]
Similarly to the Hermitian and skew-Hermitian splitting (HSS) iteration method \cite{HSS}, if we set
\[
{\cal P}_{\alpha}= \frac{1}{2\alpha} (\alpha \I + \A_1) (\alpha \I + \A_2), \qquad
{\cal Q}_{\alpha}= \frac{1}{2\alpha} (\alpha \I - \A_1) (\alpha \I - \A_2),
\]
then $\A={\cal P}_{\alpha}-{\cal Q}_{\alpha}$ and
\begin{equation}\label{RoPa}
{\cal T}_\alpha= {\cal P}_{\alpha}^{-1}{\cal Q}_{\alpha}=\I-{\cal P}_{\alpha}^{-1}\A.
\end{equation}

In each iteration of the APSS iteration method, two systems with the coefficient matrices $\alpha \I + \A_1$ and $\alpha \I + \A_2$ should be solved. The first step of the $k$th iteration of APSS can be written as
\begin{equation}\label{InSys1}
\begin{pmatrix}
\alpha I+A & B^T & 0 \\
-B & \alpha I  & 0 \\
0 &  0  & \alpha I
\end{pmatrix}
\begin{pmatrix}
x^{(k+\frac{1}{2})} \\
y^{(k+\frac{1}{2})} \\
z^{(k+\frac{1}{2})}
\end{pmatrix}
=
\begin{pmatrix}
\alpha x^{(k)}+ f \\
\alpha y^{(k)}+C^T z^{(k)} -g \\
-Cy^{(k)}+\alpha z^{(k)}+ h
\end{pmatrix}
\equiv
\begin{pmatrix}
f^{(k)} \\
g^{(k)} \\
h^{(k)}
\end{pmatrix}.
\end{equation}
From Eq. \eqref{InSys1}, we immediately get $z^{(k+\frac{1}{2})}=h^{(k)}/\alpha$ and the reduced system
\begin{equation}\label{InSys1-Reds}
\begin{pmatrix}
\alpha I+A & B^T  \\
-B & \alpha I  \\
\end{pmatrix}
\begin{pmatrix}
x^{(k+\frac{1}{2})} \\
y^{(k+\frac{1}{2})}
\end{pmatrix}
=
\begin{pmatrix}
f^{(k)} \\
g^{(k)}
\end{pmatrix},
\end{equation}
From the second equation in \eqref{InSys1-Reds}, we obtain
\[
y^{(k+\frac{1}{2})}=\frac{1}{\alpha}(g^{(k)}+Bx^{(k+\frac{1}{2})}),
\]
and substituting it into the first equation of  \eqref{InSys1-Reds} gives the system
\[
(\alpha I + A + \frac{1}{\alpha} B^TB) x^{(k+\frac{1}{2})}= f^{(k)}-\frac{1}{\alpha} B^T g^{(k)},
\]
for $x^{(k+\frac{1}{2})}$. Summarizing the above results gives the following algorithm for computing $\x^{(k+\frac{1}{2})}$:
\begin{algorithm}\label{Alg1} Computation of $\x^{(k+\frac{1}{2})}$:\\
	$1$. Set $z^{(k+\frac{1}{2})}=h^{(k)}/\alpha$.\\
	$2$. Solve $(\alpha I + A + \frac{1}{\alpha} B^TB) x^{(k+\frac{1}{2})}= f^{(k)}-\frac{1}{\alpha} B^T g^{(k)}$ for $x^{(k+\frac{1}{2})}$. \\
	$3$. Set $y^{(k+\frac{1}{2})}=\frac{1}{\alpha}(g^{(k)}+Bx^{(k+\frac{1}{2})})$.\\
	$4$. Set $\x^{(k+\frac{1}{2})}=(x^{(k+\frac{1}{2})}; y^{(k+\frac{1}{2})}; z^{(k+\frac{1}{2})})$.
\end{algorithm}

The second step of the $k$th iteration of APSS can be written as
\begin{eqnarray*}\label{InSys2}
\begin{split}
\begin{pmatrix}
\alpha I &    0      & 0 \\
0  & \alpha I  & -C^T \\
0 &  C        & \alpha I
\end{pmatrix}
\begin{pmatrix}
x^{(k+1} \\
y^{(k+1)} \\
z^{(k+1)}
\end{pmatrix}
&=
\begin{pmatrix}
(\alpha I-A) x^{(k+\frac{1}{2})}-B^T y^{(k+\frac{1}{2})}+ f \\
B x^{(k+\frac{1}{2})} +  \alpha y^{(k+\frac{1}{2})} -g \\
\alpha z^{(k+\frac{1}{2})}+ h
\end{pmatrix}\\
&\equiv
\begin{pmatrix}
f^{(k+\frac{1}{2})} \\
g^{(k+\frac{1}{2})} \\
h^{(k+\frac{1}{2})}
\end{pmatrix}.
\end{split}
\end{eqnarray*}
From the first equation of the above system we get $x^{(k+1}=f^{(k+\frac{1}{2})}/\alpha$. After eliminating  the variable
$x^{(k+1)}$ from the above system we find the following reduced system
\begin{equation}\label{InSys2-Red}
\begin{pmatrix}
\alpha I  & -C^T \\
C        & \alpha I
\end{pmatrix}
\begin{pmatrix}
y^{(k+1)} \\
z^{(k+1)}
\end{pmatrix}
=
\begin{pmatrix}
g^{(k+\frac{1}{2})} \\
h^{(k+\frac{1}{2})}
\end{pmatrix}.
\end{equation}
From the first equation in \eqref{InSys2-Red} we deduce that
\[
y^{(k+1)}=\frac{1}{\alpha}(g^{(k+\frac{1}{2})}+C^T z^{(k+1)} ),
\]
and substituting $y^{(k+1)}$ in the second equation of \eqref{InSys2-Red} results in the following system for $z^{(k+1)}$:
\[
(\alpha I+\frac{1}{\alpha} CC^T) z^{(k+1)} = h^{(k+\frac{1}{2})} - \frac{1}{\alpha} C g^{(k+\frac{1}{2})}.
\]
Summarizing the above results gives the following algorithm for computing $\x^{(k+1)}$.
\begin{algorithm}\label{Alg2} Computation of $\x^{(k+1)}$:\\
	$1$. Set $x^{(k+1)}=f^{(k+\frac{1}{2})}/\alpha$.\\
	$2$. Solve $(\alpha I+\frac{1}{\alpha} CC^T) z^{(k+1)} = h^{(k+\frac{1}{2})} - \frac{1}{\alpha} C g^{(k+\frac{1}{2})}$ for $z^{(k+1)}$. \\
	$3$. Set $y^{(k+1)}=\frac{1}{\alpha}(g^{(k+\frac{1}{2})}+C^T z^{(k+1)} )$. \\
	$4$. Set $\x^{(k+1)}=(x^{(k+1)}; y^{(k+1)}; z^{(k+1)})$.
\end{algorithm}


From algorithms \ref{Alg1} and \ref{Alg2} we see that two systems with the coefficient matrices $\alpha I + A + \frac{1}{\alpha} B^TB$ and $\alpha I+\frac{1}{\alpha} CC^T$ should be solved. Both of these matrices are symmetric positive definite. Therefore, the corresponding systems can be solved exactly using the Cholesky factorization, or inexactly using the conjugate gradient (CG) method.

\section{Convergence of the APSS iteration method}\label{Sec4}

It is well-known that the stationary iteration method \eqref{DS-Stat} converges to the solution of system \eqref{TSys} for any initial guess ${\bf x}^{(0)}$ if and only if $\rho({\cal T}_{\alpha})<1$. To prove the convergence, we first state the next lemma.
\begin{lemma}[Kellogg's lemma, see \cite{Marchuk}] \label{Kellogg}
If $A\in\Bbb{C}^{n\times n}$ is positive semidefinite, then
	\[
	\|(\alpha I+A)^{-1} (\alpha I-A)\|_2 \leq 1,
	\]  	
	for all $\alpha>0$. Moreover, if $A$ is positive definite, then
	\[
	\|(\alpha I+A)^{-1} (\alpha I-A)\|_2 < 1,
	\]
	for all $\alpha>0$.	
\end{lemma}

\begin{lemma}\label{Rolte1} Let $A$ be SPD and let matrices $B$ and $C$ be of full row rank. Then the spectral radius of the iteration matrix  ${\cal T}_{\alpha}$ in \eqref{DS-Stat} is less than or equal to unity, that is,	$\rho({\cal T}_{\alpha})\leq 1$, in which $\rho(.)$ denotes the spectral radius of the matrix.
\end{lemma}	
\begin{proof}
		Clearly, the matrix ${\cal T}_{\alpha}$ is similar to
		\[
		\hat{\cal T}_{\alpha}=\left( \alpha \I - \A_1 \right)\left( \alpha \I + \A_1 \right)^{-1}\left( \alpha \I -\A_2 \right) \left( \alpha \I + \A_2 \right)^{-1}.
		\]
		So,
		\begin{eqnarray*}
			\rho({\cal T}_{\alpha}) &=& \rho(\hat{\cal T}_{\alpha})
			\leq  \left\| \left( \alpha \I - \A_1 \right)\left( \alpha \I + \A_1 \right)^{-1}\right\|_2
			\left\|\left( \alpha \I -\A_2 \right) \left( \alpha \I + \A_2 \right)^{-1} \right\|_2.
		\end{eqnarray*}
		Since both matrices $\A_1$ and $\A_2$ are positive semidefinite, it follows from the Kellogg's lemma that
		\[
		\left\| \left( \alpha \I - \A_1 \right)\left( \alpha \I + \A_1 \right)^{-1}\right\|_2\leq1,\quad
		\left\|\left( \alpha \I -\A_2 \right) \left( \alpha \I + \A_2 \right)^{-1} \right\|_2\leq 1.
		\]
		Therefore, we conclude that $\rho({\cal T}_{\alpha}) = \rho(\hat{\cal T}_{\alpha})\leq 1$.
	\end{proof}	


It follows from Lemma \ref{Rolte1} that for the convergence of the APSS iteration method it is enough to prove that
$\rho({\cal T}_{\alpha})=1$ never happens. To do so, we state the following lemma.

\begin{lemma}\label{Lemma3}
	Let $A$ be SPD and matrices $B$ and $C$ be of full row rank. Then the following are equivalent:\\
\begin{enumerate}	
\item[$($i$)$]
The matrix
	\[
	{\cal G}_\alpha=\begin{pmatrix}
	A &    B^T+\frac{1}{\alpha^2} B^TC^TC      & 0 \\
	-B &    0      & -C^T \\
	0 &    C      & 0
	\end{pmatrix}
	\] 	
	does not have any purely imaginary eigenvalue.\\
\item[$($ii$)$]
The spectral radius of the iteration matrix  ${\cal T}_{\alpha}$ in \eqref{DS-Stat} is strictly less than 1, that is,	$\rho({\cal T}_{\alpha})< 1$.
\end{enumerate}
\end{lemma}	
\begin{proof}
		Similarly to the proof of Lemma 2 in \cite{DS-Benzi}, let $\lambda$ be an eigenvalue of ${\cal T}_\alpha$. Then, it follows from Eq. \eqref{RoPa} that $\lambda=1-\mu$,  where $\mu$ is an eigenvalue of the matrix ${\cal P}_{\alpha}^{-1}\A$.
		Let $(\mu,\x)$ be an eigenpair of ${\cal P}_{\alpha}^{-1}\A$.  Then, we have $\A \x=\mu {\cal P}_{\alpha} \x $, which is equivalent to
		\[
		\A \x = \frac{\mu}{2 \alpha} (\A_1 \A_2 +\alpha \A+\alpha^2 \I)\x,
		\]
		or
		\begin{equation}\label{AA1A2}
		\left( 1-\frac{1}{2} \mu \right) \A \x =\frac{\mu \alpha}{2} (\I+\frac{1}{\alpha^2} \A_1 \A_2) \x.
		\end{equation}
		Direct computation reveals that
		\[
		{\cal H}_{\alpha}:=\I+\frac{1}{\alpha^2} \A_1 \A_2= \begin{pmatrix}
		I &  0  &-\frac{1}{\alpha^2} B^TC^T \\
		0 &  I  &  0 \\
		0 &  0  &  I
		\end{pmatrix},
		\]
		which is obviously nonsingular. Since both matrices $\A$ and ${\cal H}_{\alpha}$ are nonsingular we deduce that
		$\mu \neq 0$ and $\mu\neq 2$.  Therefore, from \eqref{AA1A2} we have
		\[
		{\cal H}_{\alpha}^{-1} \A \x =\frac{\mu \alpha}{2-\mu} \x.
		\]
		This shows that
		\begin{equation}\label{tetmu}
		\theta:=\frac{\mu \alpha}{2-\mu}
		\end{equation}
		is an eigenvalue of
		\begin{equation}\label{Galpha}
		{\cal H}_{\alpha}^{-1} \A=\begin{pmatrix}
		A &    B^T+\frac{1}{\alpha^2} B^TC^TC      & 0 \\
		-B &    0      & -C^T \\
		0 &    C      & 0
		\end{pmatrix}={\cal G}_{\alpha}.
		\end{equation}
		From \eqref{tetmu} we see that
		\[
		\mu=\frac{2\theta}{\alpha+\theta},
		\]
		and as a result
		\[
		\lambda=1-\mu=1-\frac{2\theta}{\alpha+\theta}=\frac{\alpha-\theta}{\alpha+\theta}.
		\]
		Now from Lemma \ref{Rolte1} we have
		\[
		|\lambda|=\left|\frac{\alpha-\theta}{\alpha+\theta}\right|\leq 1,
		\]
		and $|\lambda|=1$ if and only if
		$|{\theta-\alpha}|=|{\theta+\alpha}|$, which is itself equivalent to
		\[
		(\Re(\theta)-\alpha)^2+\Im(\theta)^2= (\Re(\theta)+\alpha)^2+\Im(\theta)^2.
		\]
		The latter equation is equivalent to $\Re(\theta)=0$. Therefore, $\rho({\cal T}_{\alpha})=1$ if and only if
		${\cal G}_{\alpha}$ has at least one purely imaginary eigenvalue.
	\end{proof}
%=========================================================================
\begin{theorem}\label{Thm2}
	Let $A$ be SPD and matrices $B$ and let $C$ be of full row rank. Then the spectral radius of the iteration matrix  ${\cal T}_{\alpha}$ in \eqref{DS-Stat} is strictly less than 1, that is,	
	\[
	\rho({\cal T}_{\alpha})<1,\quad \forall \alpha>0,
	\]
	which means that the APSS iteration method unconditionally converges to the solution of \eqref{TSys}.
\end{theorem}
\begin{proof}
		According to Lemma \ref{Lemma3}, all we need is to prove that the matrix ${\cal G}_{\alpha}$ defined in
		\eqref{Galpha} has no purely imaginary eigenvalue. Let $(\t,\x)$ be an eigenpair of the matrix ${\cal G}_{\alpha}$
		with $\|\x\|_2=1$. Clearly, the matrix ${\cal G}_{\alpha}$ is nonsingular, therefore $\t\neq 0$. Letting $\x=(u;v;p)$, it follows from  ${\cal G}_{\alpha}\x=\t \x$ that
		\begin{eqnarray}
		Au+B^T(I+\frac{1}{\alpha^2}C^TC)v &\h=\h&\t u, \label{Tet1}\\
		-Bu-C^Tp &\h=\h&\t v, \label{Tet2} \\
		Cv &\h=\h& \t p.  \label{Tet3}
		\end{eqnarray}
		
		
		If $u=0$, then it follows from  $B^T(I+\frac{1}{\alpha^2}C^TC)v=0$ that $v=0$, since $B^T$ is of full column rank and the matrix $I+\frac{1}{\alpha^2}C^TC$ is SPD. Substituting $v=0$ in \eqref{Tet3} yields $p=0$. Therefore, $\x=(u;v;p)=0$, which is impossible since $\x$ is an eigenvector. Therefore, $u\neq 0$.
		
		If $v=0$, then from Eq. \eqref{Tet1} we obtain $Au=\t u$. Thus,
		\[
		\t=\frac{1}{\|u\|_2^2} u^*Au>0,
		\]
		since $u\neq 0$ and $A$ is SPD. This means that $\Re(\t)=\t>0$.
		
		If $p=0$, then from Eq. \eqref{Tet3} we get $Cv=0$. Therefore, from eqs. \eqref{Tet1} and \eqref{Tet2} we obtain the eigenvalue problem
		\begin{equation}\label{Thm-Red}
		\begin{pmatrix}
		A & B^T  \\
		-B & 0  \\
		\end{pmatrix}
		\begin{pmatrix}
		u \\
		v
		\end{pmatrix}
		=
		\t \begin{pmatrix}
		u \\
		v
		\end{pmatrix}.
		\end{equation}
		Eq. \eqref{Thm-Red} shows that $\t$ is an eigenvalue of the saddle matrix
		\begin{equation*}\label{InSys1-Red}
		{\cal S}:=\begin{pmatrix}
		A & B^T  \\
		-B & 0  \\
		\end{pmatrix}.
		\end{equation*}
		Therefore, since $A$ is SPD and $B$ is of full rank, it follows  from  \cite[Lemma 1.1 (iv)]{SD-Benzi-SIMAX}  that $\Re(\t)>0$.
		
		According to the above results, here we assume that $u\neq 0$, $v\neq 0$ and $p \neq 0$.  From ${\cal G}_{\alpha}\x=\t \x$ and $\|\x\|_2=1$, we get $\t=\x^* {\cal G}_{\alpha} \x$. Hence,
		\begin{eqnarray}
		\nonumber\Re(\t) &\h=\h&\frac{1}{2} \x^* ({\cal G}_\alpha+ {\cal G}_\alpha^T) \x
		=\frac{1}{2}\begin{pmatrix}
		u^* & v^* & p^*
		\end{pmatrix}
		\begin{pmatrix}
		2A &  \frac{1}{\alpha^2} B^TC^TC & 0\\
		\frac{1}{\alpha^2} C^TCB  &  0  &  0 \\
		0  &  0  &  0
		\end{pmatrix} \begin{pmatrix}
		u \\
		v \\
		p
		\end{pmatrix} \\
		&\h=\h& u^*Au +\frac{1}{2\alpha^2} \left(u^*B^TC^TCv+v^*C^TCBu\right). \label{Rt1}
		\end{eqnarray}
		On the other hand, from Eq. \eqref{Tet1} we deduce
		\begin{equation}\label{tusu1}
		u^*Au+u^*B^T(I+\frac{1}{\alpha^2} C^TC)v=\t u^*u,
		\end{equation}
		and, taking the conjugate of both sides, the latter equation gives
		\begin{equation}\label{tusu2}
		u^*Au+v^*(I+\frac{1}{\alpha^2} C^TC)Bu=\bar{\t}u^*u.
		\end{equation}
		From eqs. \eqref{tusu1} and \eqref{tusu2} we get
		\begin{equation}\label{tusu3}
		u^*Au+\frac{1}{2}(u^*B^Tv+v^*Bu)+\frac{1}{2\alpha^2} \left(u^*B^TC^TCv+v^*C^TCBu\right)=\Re(\t) u^*u.
		\end{equation}
		From eqs. \eqref{tusu3}  and \eqref{Rt1} we obtain
		\[
		\Re(\t)+\frac{1}{2}(u^*B^Tv+v^*Bu)=\Re(\t) u^*u.
		\]
		Now, by contradiction we assume that $\Re(\t)=0$. In this case, from the above equation we deduce that
		\begin{equation}\label{tusu3b}
		u^*B^Tv+v^*Bu=0.
		\end{equation}
		On the other hand, we have
		\begin{eqnarray*}
			\Re(\t) &\h=\h& u^*Au +\frac{1}{2\alpha^2} \left(u^*B^TC^TCv+v^*C^TCBu\right)\\
			&\h=\h& u^*Au +\frac{1}{2\alpha^2} \left(\t u^*B^TC^Tp+ \bar{\t} p^*CBu\right) \hspace{3.2cm} \text{(from Eq. \eqref{Tet3})}\\
			&\h=\h& u^*Au +\frac{1}{2\alpha^2} \left(\t u^*B^T(-Bu-\t v)+ \bar{\t} (-Bu-\t v)^*Bu\right) \hspace{0.7cm} \text{(from Eq. \eqref{Tet2})}\\
			&\h=\h& u^*Au -\frac{1}{2\alpha^2} \left( 2\Re(\t) \|Bu\|_2^2+\t^2 u^*B^Tv+ \bar{\t}^2 v^*Bu \right).
		\end{eqnarray*}
		Now, if $\Re(\t)=0$, then $\t=i\xi$, where $\xi\neq 0$. Therefore, from the above equation we see that
		\begin{eqnarray*}
			0 &\h=\h&  u^*Au + \frac{\xi^2}{2\alpha^2} \left(u^*B^Tv+ v^*Bu \right)\\
			&\h=\h&  u^*Au, \qquad \text{(from Eq. \eqref{tusu3b})},
		\end{eqnarray*}
		which is a contraction since $u\neq 0$ and $A$ is SPD. Therefore, the proof is completed.
	\end{proof}		


\section{APSS preconditioner and its implementation issues}\label{Sec5}

In the previous section we have proved that under the conditions of Theorem \ref{Thm2}  we have $\rho({\cal T}_{\alpha})<1$. Therefore, the eigenvalues of ${\cal T}_{\alpha}$ are included in a circle centered at the origin with radius 1. On the other hand, from \eqref{RoPa} we have
\[
{\cal P}_{\alpha}^{-1}\A=\I- {\cal T}_\alpha.
\]
Hence, the eigenvalues of ${\cal P}_{\alpha}^{-1}\A$ are contained in a circle centered at $(1,0)$ with radius 1. Therefore, ${\cal P}_{\alpha}$ can be used as a preconditioner for the system $\A \x=\b$. Since the pre-factor $1/(2\alpha)$ in the matrix  ${\cal P}_{\alpha}$ has no effect on the preconditioned system, we drop it. Hence, the matrix
\[
{\cal M}_{\alpha}= (\alpha \I + \A_1) (\alpha \I + \A_2)
\]
can be used to accelerate the convergence of the Krylov subspace methods like GMRES or its restarted version to solve the system  $\A \x=\b$.
In each iteration of the GMRES method a linear system of the equations of the form ${\cal M}_{\alpha}v=r$ should be solved.
This can be done in two steps. In the first step, the system
\[
(\alpha \I + \A_1)w=r
\]
is solved for $w$ and then the system
\[
(\alpha \I + \A_2)v=w
\]
for $v$. These systems can be solved using algorithms \ref{Alg1} and \ref{Alg2}.

\begin{remark}
	According to part (i) of Lemma \ref{Cor2}, if the matrix $A$ is symmetric positive semidefinite and the matrices  $B$  and $C$ are of  full row rank with ${\cal N}(A) \cap {\cal N}(B)=\{0\}$,  then the matrix ${\cal A}$ is nonsingular. In this case, the proposed iteration method is still well-defined; however, convergence of the method is not guaranteed. Nevertheless, from Lemma \ref{Rolte1} we have $\rho({\cal T}_{\alpha})\leq 1$ and the GMRES method can be applied to the corresponding preconditioned system. On the other hand, using the idea of \cite{SD-Benzi-SIMAX}, one may state a convergent iteration method using the extrapolation technique.
\end{remark}


\section{Numerical experiments}\label{Sec6}
We present some numerical results to show the efficiency of the proposed preconditioner. We present three examples
and in all examples we first apply a symmetric diagonal scaling for the matrix $\A$. To do so, we replace the coefficient matrix
${\A}$ by the matrix ${\cal D}^{-\frac{1}{2}} \A {\cal D}^{-\frac{1}{2}} $, where ${\cal D}=\diag (\| {\A}_{1}\|_2,\ldots,\|{\A}_{\bf  n}\|_2)$, in which
${\A}_{j}$ is the $j$th column of the matrix $\A$. The right-hand side vector of the system is set to $\b={\cal A}{\bf e}$, where ${\bf e}$ is a vector of all ones. We use the flexible version of the GMRES(50) \cite{SaadBook,FGMRES}, FGMRES(50), for solving the systems.
The iteration is started from a zero vector and terminated as soon as the residual 2-norm is reduced by a factor of $10^{-6}$. The maximum number of iterations is set to be 20\,000.

For the  APSS preconditioner we use the implementation method presented in Section \ref{Sec5}. The subsystems are solved using the conjugate gradient (CG) method. The CG method is started from a zero vector and the iteration is stopped as soon as the residual 2-norm is reduced by a factor of $10^3$. The maximum number of  CG iterations is set to be 200. We compare the numerical results of the APSS preconditioner with those of the  preconditioners ${\cal P}_D$  and ${\cal P}_2$ defined in eqs. \eqref{PDPrecon} and \eqref{P123}, respectively.
In the implementation of the preconditioners ${\cal P}_D$  and ${\cal P}_2$ in the FGMRES(50) method three systems with the coefficient matrices $A$, $S=BA^{-1}B$ and $M=CS^{-1}C^T$ should be solved, which are symmetric positive definite. The system with the coefficient matrix $A$ can be solved using the CG method; however, there is no  obvious implementation method to solve inexactly the systems with the coefficient matrices $S$ and $CS^{-1}C^T$. Hence, using the idea of \cite{Huang-NumerAlgor}, we set
\[
\hat{A}=LL^T,
\]
where $L$ is computed using the incomplete Cholesky factorization with dropping tolerance $10^{-8}$. To do so, we use the \textsc{Matlab} command:
\[
\verb"L=ichol(A, struct('type','ict','droptol',1e-8))".
\]
We also set
\[
\hat{S}=\diag(B\hat{A}^{-1}B^T).
\]
Now, instead of using the preconditioners ${\cal P}_D$ and ${\cal P}_2$,  we use
\[
\hat{{\cal P}}_D=\begin{pmatrix}
\hat{A} &  0  & 0 \\
0 &  \hat{S}  & 0 \\
0 &  0   & C\hat{S}^{-1}C^T
\end{pmatrix}, \quad
\hat{{\cal P}}_2=\begin{pmatrix}
\hat{A} &  0  & 0 \\
B &  -\hat{S}  & C^T \\
0 &  0   & C\hat{S}^{-1}C^T
\end{pmatrix}, \quad
\]
respectively. In the implementation of these preconditioners the matrix $C\hat{S}^{-1}C^T$ is formed explicitly  and the corresponding system is solved using the Cholesky factorization. All runs are implemented in \textsc{Matlab} R2017, equipped with a laptop with 1.80 GHz central processing unit (Intel(R) Core(TM) i7-4500), 6 GB memory and Windows 7 operating system.

Numerical results are presented in the tables. In the tables, ``Iters" stands for the number of iterations for the convergence. The elapsed CPU time (in seconds) is denoted by ``Time". The value of $R_k$ defined by
\[
R_k=\frac{\|\b -\A \x^{(k)}\|_2}{\|\b\|_2}
\]
is also reported, where  $\x^{(k)}$ is the computed solution at iteration  $k$. In the tables, $nnz$ stands for the number of nonzero entries of the matrix. A dagger ($\dag$) shows that the  method has not converged in the maximum number of iterations. Also,  $\S$ shows that
the preconditioner has not been computed in  500 seconds.
Finally, ``NA" (for not applicable) means that the coefficient matrix does not satisfy the assumptions of Theorem \ref{Thm2}.
We now present the examples.

\begin{example} \label{Ex1}
	We consider problem \eqref{TSys} with (see \cite{Huang-NumerAlgor,Xie-Li-CAMWA})
	\begin{align*}
	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}},\\
	B&=(I \otimes F,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}},
	\end{align*}
	where
	\begin{align*}
	T=\frac{1}{h^{2}}  \operatorname{tridiag}(-1,2,-1) \in \mathbb{R}^{p \times p}, \quad F=\frac{1}{h}  \operatorname{tridiag}(0,1,-1) \in \mathbb{R}^{p \times p},
	\end{align*}
	and $
	E=\operatorname{diag}\left(1, p+1, \ldots, p^{2}-p+1\right)
	$ in which $\otimes$ denotes the Kronecker product and $h={1}/{(p+1)}$ the discretization meshsize.
	We present the numerical results for $p=16,32,64,128,$ and $256$. We may use the idea of \cite{CaoCAM} to estimate the optimal value of the parameter $\alpha$. However, for the APSS preconditioner the value of $\alpha$ is set to be $0.005$ for all dimensions in this example. We have experimentally found that this choice often gives quite suitable results for this problem. It is necessary to mention that $\alpha=0.005$ is not the optimal value of the parameter. Numerical results are presented in Table \ref{Tbl1}.
	%======================================================================================================================     	
	\begin{table}[h!]
\centering
			\begin{tabular}{|p{2.7cm}|p{1.5cm}|p{1.5cm}|p{1.5cm}|p{1.5cm}|p{1.5cm}|}\hline
				$p,{\bf n},nnz$  & & No Prec.& $\mathcal{M}_{\alpha}$ & $\hat{\mathcal{P}}_{D}$ & $\hat{\mathcal{P}}_{2}$ \\	\hline
				$p=16$          &  Iters     &  425    & 15     & 29     & 22     \\
				${\bf n}=1024$  &  Time      &  4.02   & 0.04   & 0.04   & 0.04   \\
				$nnz=5408$      &  $R_k$     & 8.6e-7  & 3.1e-7 & 5.9e-7 & 8.8e-7 \\  \hline
				
				$p=32$          &  Iters     &  949    & 13     & 32     & 25     \\
				${\bf n}=4096$  &  Time      & 83.43   & 0.09   & 0.32   & 0.25   \\
				$nnz=22080$     &  $R_k$     & 9.9e-7  & 5.3e-7 & 6.9e-7 & 7.5e-7 \\  \hline
				
				$p=64$           &  Iters     & $\dag$  & 13     & 34     & 28     \\
				${\bf n}=16384$  &  Time      &   -     & 0.47   & 8.69   & 5.99   \\
				$nnz=89216$      &  $R_k$     &   -     & 9.1e-7 & 7.3e-7 & 5.7e-7 \\  \hline
				
				$p=128$          &  Iters     & $\dag$  & 22     & 68     & 69     \\
				${\bf n}=65536$  &  Time      &   -     & 3.10   & 174.16 & 175.89   \\
				$nnz=358656$     &  $R_k$     &   -     & 6.4e-7 & 4.8e-7 & 5.6e-7 \\  \hline
				
				$p=256$          &  Iters     & $\dag$  & 51     & $\S$   & $\S$     \\
				${\bf n}=262144$ &  Time      &   -     & 37.47  & -      &   -   \\
				$nnz=1438208$    &  $R_k$     &   -     & 6.6e-7 & -      &   -  \\  \hline
				
			\end{tabular}
		\caption{\it Numerical results for FGMRES(50) to solve Example \ref{Ex1}}
\label{Tbl1}
	\end{table}

\begin{figure}[h!]
\centering
			\includegraphics[height=4.5cm,width=5.5cm]{Fig1a.eps}\includegraphics[height=4.5cm,width=5.5cm]{Fig1b.eps}\\
			\includegraphics[height=4.5cm,width=5.5cm]{Fig1c.eps}\includegraphics[height=4.5cm,width=5.5cm]{Fig1d.eps}
		\caption{\it Eigenvalue distribution ${\cal A}$ ,  ${\mathcal M}_{\alpha}^{-1} \mathcal{A}$,  $\hat{\mathcal P}_{D}^{-1} \mathcal{A}$ and  $\hat{\mathcal P}_2^{-1} \mathcal{A}$  with $p=8$ for Example \ref{Ex1}} \label{Fig1} 	
\end{figure}

\newpage
\noindent
Numerical results show that the three preconditioners significantly reduce the number of iterations and the CPU time of FGMRES(50). However, the APSS preconditioner outperforms the others from both the number of iterations and the CPU time point of view. The main problem with the preconditioners
	$\hat{{\cal P}}_D$ and $\hat{{\cal P}}_2$ is their need for a large amount of CPU time to set up.

Fig. \ref{Fig1} displays the eigenvalues of matrices ${\cal A}$ ,  ${\mathcal M}_{\alpha}^{-1} \mathcal{A}$,  $\hat{\mathcal P}_{D}^{-1} \mathcal{A}$ and  $\hat{\mathcal P}_2^{-1} \mathcal{A}$  for $p=8$. In this case, the matrix $A$ is of dimension ${\bf n}=256$. It is seen that the eigenvalues of ${\mathcal M}_{\alpha}^{-1} \mathcal{A}$ are clustered   in a circle centered at $(1,0)$ with radius 1. Another observation which can be posed here is that the eigenvalues of ${\mathcal M}_{\alpha}^{-1} \mathcal{A}$ are more clustered around the point $(1,0)$ than the others.
\end{example}


\begin{example}\label{Ex2}
	We consider the three-by-three block saddle point problem \eqref{TSys}, where  (see \cite{Huang-NumerAlgor,Xie-Li-CAMWA})
	\begin{align*}
	A=\operatorname{diag}\left(2 W^{T} W+D_{1}, D_{2}, D_{3}\right) \in \mathbb{R}^{n \times n},
	\end{align*}
	is a block-diagonal matrix,
	\begin{align*}
	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{align*}
	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{align*}
	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.
	\\[1ex]
	d_{j}^{(3)}&=10^{-5}(j+\tilde{p})^{2}
	\text { for } 1 \leq j \leq 2 \tilde{p},
	\end{align*}
	and
	\begin{align*}
	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{align*}
	The numerical results are  listed in Table \ref{Tbl2} for different values of $p$ ($p=16,32,64,128,$ $256$).
	In the APSS preconditioner, the value of $\alpha$ is set to be $0.4$ for all dimensions. We have experimentally found  that this value often gives good results for this problem.
	
	As the numerical results show, the APSS preconditioner improves significantly the number of the iterations and the CPU time of the GMRES method. As we see, the preconditioner $\hat{\mathcal P}_{D}$ improves neither the number of iterations nor the CPU time of the GMRES for all sizes, except for $p=16$.  On the other hand, $\hat{\mathcal P}_2$ improves slightly the number of iterations; however, this is not the case for the CPU time. Therefore, the APSS preconditioner is superior to the others from both the number of iterations and the CPU time point of view.
	
Fig. \ref{Fig2} shows the eigenvalues of matrices ${\cal A}$ ,  ${\mathcal M}_{\alpha}^{-1} \mathcal{A}$,  $\hat{\mathcal P}_{D}^{-1} \mathcal{A}$ and  $\hat{\mathcal P}_2^{-1} \mathcal{A}$  for $p=8$. In this case, the matrix $A$ is of dimension ${\bf n}=256$. It is seen that the eigenvalues of ${\mathcal M}_{\alpha}^{-1} \mathcal{A}$ are clustered   in a circle centered at $(1,0)$ with radius 1. We also see that the eigenvalues of ${\mathcal M}_{\alpha}^{-1} \mathcal{A}$ are more clustered around the point $(1,0)$ than the others.
	
	\begin{table}[h!]
		\centering
			\begin{tabular}{|p{2.7cm}|p{1.5cm}|p{1.5cm}|p{1.5cm}|p{1.5cm}|p{1.5cm}|}\hline
				$p,{\bf n},nnz$  & & No Prec.& $\mathcal{M}_{\alpha}$ & $\mathcal{P}_{D}$ & $\mathcal{P}_{2}$ \\	\hline
				$p=16$          &  Iters     & 150     & 31     & 96      & 70     \\
				${\bf n}=2080$  &  Time      &  0.08   & 0.04   & 0.07    & 0.06   \\
				$nnz=9972$      &  $R_k$     & 9.4e-7  & 8.8e-7 & 1.0e-6  & 9.8e-7 \\  \hline
				
				$p=32$          &  Iters     &  176    & 32     & 222     & 158     \\
				${\bf n}=4096$  &  Time      & 0.36    & 0.14   & 0.58    & 0.47   \\
				$nnz=22080$     &  $R_k$     & 1.0e-6  & 8.0e-7 & 9.3e-7  & 9.9e-7 \\  \hline
				
				$p=64$           &  Iters     & 213    & 31     & 275     & 186     \\
				${\bf n}=16384$  &  Time      & 1.53   & 0.42   & 3.16    & 2.59   \\
				$nnz=89216$      &  $R_k$     & 9.8e-7 & 9.9e-7 & 1.0e-6  & 9.8e-7 \\  \hline
				
				$p=128$          &  Iters     & 204    & 30     & 263     & 178     \\
				${\bf n}=65536$  &  Time      & 7.59   & 1.67   & 26.22   & 21.94   \\
				$nnz=358656$     &  $R_k$     & 9.7e-7 & 8.4e-7 & 9.9e-7  & 9.8e-7 \\  \hline
				
				$p=256$          &  Iters     & 195    & 29     & 248     & 167    \\
				${\bf n}=262144$ &  Time      & 32.01  & 8.41   & 853.62  & 832.01 \\
				$nnz=1438208$    &  $R_k$     & 9.6e-7 & 7.5e-7 & 9.7e-7  & 9.6e-7    \\  \hline
			\end{tabular}
\caption{\it Numerical results of FGMRES(50) to solve Example \ref{Ex2}}
\label{Tbl2}
	\end{table}
	
\begin{figure}[h!]
		\centering
			\includegraphics[height=4.5cm,width=5.5cm]{Fig2a.eps}\includegraphics[height=4.5cm,width=5.5cm]{Fig2b.eps}\\
			\includegraphics[height=4.5cm,width=5.5cm]{Fig2c.eps}\includegraphics[height=4.5cm,width=5.5cm]{Fig2d.eps}
\caption{\it Eigenvalue distribution ${\cal A}$ ,  ${\mathcal M}_{\alpha}^{-1} \mathcal{A}$,  $\hat{\mathcal P}_{D}^{-1} \mathcal{A}$ and  $\hat{\mathcal P}_2^{-1} \mathcal{A}$  with $p=8$ for Example \ref{Ex2}}
\label{Fig2}	
\end{figure}
\end{example}

\begin{example}\label{Ex3}
	In this example, we consider a quadratic program \cite{HanSIAM,Huang-NumerAlgor}
	\begin{eqnarray}
	&\h\h& \min_{x\in\Bbb{R}^n,y\in\Bbb{R}^l} \frac{1}{2}x^TAx+r^Tx+q^Ty \label{CUTER}\\
	\nonumber &\h\h&~~~s.t.:~~~ Bx+C^Ty=b,
	\end{eqnarray}
	where $r\in\Bbb{R}^n$ and $q\in\Bbb{R}^l$. To solve the above problem we define the Lagrange function
	\[
	L(x, y, \lambda) =\frac{1}{2}x^TAx+r^Tx+q^Ty+\lambda^T(Bx+C^Ty-b),
	\]
	where the vector $\lambda\in\Bbb{R}^{m}$ is the Lagrange multiplier. Then the Karush-Kuhn-Tucker necessary conditions of \eqref{CUTER} are as follows (see \cite{Bertsekas}):
	\[
	\nabla_{x} L(x, y, \lambda)=0,\quad \nabla_{y} L(x, y, \lambda)=0 \quad \text{and} \quad \nabla_{\lambda} L(x, y, \lambda)=0.
	\]
	It is easy to see that these equations give a system of linear equations of the form \eqref{MSys}. In this example, we have chosen the matrices $A$, $B$ and $C$ from the CUTEr collection \cite{CUTErCol}. Numerical results along with the generic properties of the test matrices are presented in Table \ref{Tbl3}. For the test matrices MOSARQP1 and STCQP2, the matrix $C$ is not full rank. Hence,  the matrix  $C\hat{S}^{-1}C^T$ is symmetric positive semidefinite, which shows that the preconditioners ${\cal P}_D$
	and ${\cal P}_2$ are singular. Therefore, they can not be used for preconditioning the system. Similarly, for the test matrices UBH1, AUG2D, CONT201 and  CONT300, the matrix $A$ is symmetric positive semidefinte. Hence the matrix $S$ as well as the matrices ${\cal P}_D$ and ${\cal P}_2$  can not be formed. For all test matrices in this example, the parameter $\alpha$ is the optimal one (the one with a minimum number of iterations) which has been found experimentally.
	
	The numerical results presented in Table \ref{Tbl3} show that the APSS preconditioner is very effective in reducing both the CPU time and the number of iterations for the convergence. Moreover, for the matrices YAO and liswet12 we see that the APSS preconditioner outperforms the two other preconditioners. It is necessary to mention that, although in some cases the matrices $A$ and $C$ are not of full rank, the APSS preconditioner works well.

\newpage
\noindent
\begin{table}[h!]
		\centering
			\begin{tabular}{|p{2.7cm}|p{1.5cm}|p{1.5cm}|p{1.5cm}|p{1.5cm}|p{1.5cm}|}\hline
				Matrix, ${\bf n},nnz$  & & No Prec.& $\mathcal{M}_{\alpha}$ & $\mathcal{P}_{D}$ & $\mathcal{P}_{2}$ \\	\hline
				YAO             &  $\alpha$  & -       & 0.5    & -        & -     \\
				${\bf n}=6004$  &  Iters     & 85      & 24     & 41       & 33     \\
				$nnz=18006$     &  Time      &  0.17   & 0.07   & 0.08     & 0.07   \\
				&  $R_k$     & 7.7e-7  & 8.6e-7 & 7.6e-7   & 8.2e-7 \\  \hline
				
				
				MOSARQP1        &  $\alpha$  & -       & 0.05   & -       & -     \\
				${\bf n}=5700$  &  Iters     & 115     & 6      & NA      & NA     \\
				$nnz=14434$     &  Time      &  0.32   & 0.03   & -       & -   \\
				&  $R_k$     & 8.9e-7  & 5.2e-7 & -       & - \\  \hline
				
				STCQP2          &  $\alpha$  & -       & 0.25   & -       & -     \\
				${\bf n}=10246$ &  Iters     & 316     & 45     & NA      & NA     \\
				$nnz=83979$     &  Time      &  1.30   & 0.48   & -       & -   \\
				&  $R_k$     & 9.9e-7  & 7.1e-7 & -       & - \\  \hline
				
				
				liswet12        &  $\alpha$  & -       & 0.5    & -       & -     \\
				${\bf n}=30004$ &  Iters     & 77      & 21     & 38      & 29     \\
				$nnz=90006$     &  Time      &  0.82   & 0.22   & 1.28    & 1.04   \\
				&  $R_k$     & 9.9e-7  & 9.0e-7 & 9.1e-7  & 7.4e-7 \\  \hline
				
				UBH1            &  $\alpha$  & -       & 0.095  & -       & -     \\
				${\bf n}=48018$ &  Iters     & 15969   & 1582   & NA      & NA     \\
				$nnz=138021$    &  Time      & 261.32  & 66.78  & -       & -     \\
				&  $R_k$     & 1.0e-7  & 1.0e-6 & -       & -     \\  \hline
				
				AUG2D           &  $\alpha$  & -       & 0.5    & -       & -     \\
				${\bf n}=50400$ &  Iters     & 77      & 25     & NA      & NA     \\
				$nnz=140200$    &  Time      &  1.03   & 0.43   & -       & -   \\
				&  $R_k$     & 9.6e-7  & 8.3e-7 & -       & - \\  \hline
				
				CONT201         &  $\alpha$  & -       & 0.25    & -      & -     \\
				${\bf n}=120992$&  Iters     & 183     & 26      & NA     & NA     \\
				$nnz=489592$    &  Time      & 6.30    & 1.95    & -      & -   \\
				&  $R_k$     & 9.4e-7  & 9.5e-7  & -      & - \\  \hline
				
				CONT300         &  $\alpha$  & -       & 0.25    & -      & -     \\
				${\bf n}=271492$&  Iters     & 143     & 19      & NA     & NA     \\
				$nnz=1101892$   &  Time      & 1542    & 3.22    & -      & -   \\
				&  $R_k$     & 9.6e-7  & 7.7e-7  & -      & - \\  \hline
					
			\end{tabular}
		\caption{\it Numerical results for GMRES to solve Example \ref{Ex3}}
\label{Tbl3}
\end{table}
\end{example}



\section{Conclusions} \label{Sec7}

We have applied the APSS iteration method for solving three-by-three block saddle point problems. Convergence of the method has been proved. The induced preconditioner has been used to accelerate the convergence of the FGMRES method for solving the system.
Numerical results showed that the proposed preconditioner is very effective. Moreover, we have seen that it outperforms two recently proposed preconditioners.

\section*{Acknowledgement}

The authors would like to thank the referee for careful reading of the paper and giving useful comments and suggestions.



%%%% Bibliography  %%%%%%%%%%
\begin{thebibliography}{99}

\bibitem{HSS} {\sc Z.\,Z.\,Bai, G.\,H.\,Golub, M.\,K.\,Ng}, {\it Hermitian and skew-Hermitian splitting methods for
non-Hermitian positive definite linear systems}, SIAM J. Matrix Anal. Appl. \textbf{24}(2003), 603--626.

\bibitem{BenziJCM} {\sc M.\,Benzi}, {\it Preconditioning techniques for large linear
systems: A survey},   J. Comput. Phys.  \textbf{182}(2002) 418--477.


\bibitem{SD-Benzi-SIMAX} {\sc M.\,Benzi, G.\,H.\,Golub,} {\it A preconditioner for generalized saddle point problems}, SIAM J.
Matrix Anal. Appl. \textbf{26}(2004), 20--41.

\bibitem{Benzi-Acta} {\sc M.\,Benzi, G.\,H.\,Golub, J.\,Liesen,} {\it Numerical solution of saddle point problems}, Acta
Numer. \textbf{14}(2005), 1--137.

\bibitem{DS-Benzi} {\sc M.\,Benzi, X.\,-P.\,Guo,} {\it A dimensional split preconditioner for Stokes and linearized Navier-Stokes equations}, Appl. Numer. Math. \textbf{61}(2011), 66--76.

\bibitem{Bertsekas} {\sc D.\,P.\,Bertsekas,}  {\it Nonlinear Programming}, 2nd Ed., Athena Scientific, Athena, 1999.	

\bibitem{CaoCAM} {\sc Y.\,Cao}, {\it A block positive-semidefinite splitting preconditioner for
generalized saddle point linear systems}, J. Comput. Appl. Math. \textbf{374}(2020), 112787.


\bibitem{CaoSS} {\sc Y.\,Cao}, {\it Shift-splitting preconditioners for a class of block three-by-three saddle point problems}, Appl. Math. Lett. \textbf{96}(2019), 40--46.

\bibitem{DS-Cao} {\sc Y.\,Cao, L.\,-Q.\,Yao, M.\,-Q.\,Jiang}, {\it A modified dimensional split preconditioner for generalized
saddle point problems}, J. Comput. Appl. Math. \textbf{250}(2013), 70--82.

\bibitem{ChenSIAM} {\sc Z.\,Chen, Q.\,Du, J.\,Zou,} {\it Finite element methods with matching and nonmatching meshes for Maxwell equations with discontinuous coefficients}, SIAM J. Numer. Anal. \textbf{37}(2000), 1542--1570.

\bibitem{CUTErCol} {\sc N.\,I.\,M.\,Gould, D.\,Orban, P.\,L.\,Toint}, {\it CUTEr and SifDec, a constrained and unconstrained testing
environment}, revisited, ACM Trans. Math. Softw. \textbf{29}(2003), 373--394.

\bibitem{HanSIAM} {\sc D.\,Han, X.\,Yuan}, {\it Local linear convergence of the alternating direction method of multipliers for quadratic programs}, SIAM J. Numer. Anal. \textbf{51}(2013), 3446--3457.

\bibitem{Hu-Arxiv} {\sc K.\,B.\,Hu,  J.\,C.\,Xu},  {\it Structure-preserving finite element methods for stationary MHD models},
arXiv:1503.06160.

\bibitem{Huang-NumerAlgor} {\sc N.\,Huang}, {\it Variable parameter Uzawa method for solving a class of block three-by-three saddle point problems}, Numer. Algor.  {\textbf 85}(2020), 1233--1254.

\bibitem{HuangMa} {\sc N.\,Huang, C.\,-F.\,Ma}, {\it Spectral analysis of the preconditioned system for the $3 \times 3$ block saddle point problem}, Numer. Algor. \textbf{81}(2019), 421--444.

\bibitem{Huang-NLWA} {\sc N.\,Huang, Y.\,-H.\,Dai,  Q.\,Hu}, {\it Uzawa methods for a class of block three-by-three saddle-point problems}, Numer. Linear Algebra Appl. \textbf{26}(2019), e2265.

\bibitem{DS-Ma} {\sc Y.\,Ke, C.\,F.\,Ma}, {\it A note on ``A dimensional split preconditioner for Stokes and linearized
Navier-Stokes equations''}, Appl. Numer. Math. \textbf{108}(2016), 223--225.

\bibitem{Liang} {\sc Z.\,-Z.\,Liang, G.\,-F.\,Zhang}, {\it Alternating positive semidefinite splitting preconditioners
for double saddle point problems}, Calcolo \textbf{56}(2019), 26.


\bibitem{Marchuk} {\sc G.\,I.\,Marchuk}, {\it Methods of Numerical Mathematics}, Springer-Verlag, New York, 1984.	

\bibitem{SaadBook} {\sc Y.\,Saad}, {\it Iterative methods for sparse linear systems}, PWS Press, New York, 1995.

\bibitem{FGMRES} {\sc Y.\,Saad}, {\it A flexible inner-outer preconditioned GMRES algorithm}, SIAM J. Sci. Comput. \textbf{14}(1993), 461--469.

\bibitem{WangSS} {\sc L.\,Wang, K.\,Zhang}, {\it Generalized shift-splitting preconditioner for saddle point problems with block
three-by-three structure}, OALib. Journal \textbf{6}(2019), e5968.

\bibitem{Xie-Li-CAMWA} {\sc X.\,Xie, H.\,-B.\,Li}, {\it A note on preconditioning for the $3 \times 3$ block saddle point problem}, Comput. Math.  Appl.  \textbf{79}(2020), 3289--3296.
\end{thebibliography}
%% \end{linenumbers}
\end{document}

