Accelerated RHSS Iteration Method for Stabilized Saddle-Point Problems

Zhenghui Song, Pingping Zhang

School of Science, Chongqing University of Posts and Telecommunications, Chongqing, China.

DOI: 10.4236/jamp.2022.104069                

Abstract

For stabilized saddle-point problems, we apply the two iteration parameters idea for regularized Hermitian and skew-Hermitian splitting (RHSS) method and establish accelerated RHSS (ARHSS) iteration method. Theoretical analysis shows that the ARHSS method converges unconditionally to the unique solution of the saddle point problem. Finally, we use a numerical example to confirm the effectiveness of the method.

Keywords

Stabilized Saddle-Point Problems, Regularized Hermitian and Skew-Hermitian Splitting, Iteration Parameters, Convergence Property

1. Introduction

Consider the stabilized saddle-point problem

$Ax\equiv \left(\begin{array}{cc}B& E\\ -E^{\mathrm{<mo>\; *\; </mo>}}& C\end{array}\right)\left(\begin{array}{c}y\\ z\end{array}\right)=\left(\begin{array}{c}f\\ g\end{array}\right)\equiv b\mathrm{,}$ (1)

where $B\in {ℂ}^{p\times p}$ is a Hermitian positive definite matrix, $C\in {ℂ}^{q\times q}$ is a Hermitian positive semidefinite matrix, $E\in {ℂ}^{p\times q}$ is rectangular matrix, $E^{\mathrm{<mo>\; *\; </mo>}}\in {ℂ}^{q\times p}$ is the complex conjugate transpose of E and $f\in {ℂ}^p$ and $g\in {ℂ}^q$. In addition, the null spaces of the matrices C and E do not overlap, i.e., $null\left(C\right)\cap null\left(E\right)=\left\{0\right\}$, then in accordance with [1], we know that the stabilized saddle-point problem (1) admits a unique solution; see also [2] and its reference therein. For stabilized saddle-point problem, it has a broad background in scientific computing and engineering applications. For example, it frequently appears in Navier-Stokes equations, finite element methods, regularized and weighted least-squares problems. For more details, see [3] [4] [5] [6] [7].

When $C=0$, the stabilized saddle-point problem (1) degenerated into standard saddle-point problem. Benzi and Golub [8] proposed the Hermitian and skew-Hermitian splitting (HSS) iteration method to solve the standard saddle-point problem. This iteration method is a direct generalization of the HSS iteration method initially proposed in [9]. To improve the convergence rate of the HSS iteration method, Bai, Golub and Pan [10] proposed the preconditioned HSS (PHSS) iteration method, and then Bai and Golub [11] proposed the accelerated HSS (AHSS) iteration method by calculating the HSS iteration sequence with different parameters. Bai and Benzi [12] proposed the regularized Hermitian and skew-Hermitian splitting (RHSS) iteration method to solve the standard saddle-point problem. The biggest highlight is that the condition number of iterative systems can be improved by regularizing the matrix. In addition, Behzadi and Abdollahi [13] combined the AHSS iteration method and the RHSS iteration method to build the accelerated RHSS (ARHSS) iteration method for the standard saddle-point linear system, and experimental results show that the ARHSS iteration method is better than the HSS iteration method and the RHSS iteration method.

Recently, based on the regularized Hermitian and skew-Hermitian splitting of the stabilized saddle-point matrix A in (1):

$\begin{array}{c}A=\left(\begin{array}{cc}B& 0\\ 0& Q_1+\omega C\end{array}\right)+\left(\begin{array}{cc}0& E\\ -E^{*}& -Q_1+\left(1-\omega \right)C\end{array}\right)=H_{+}\left(\omega \right)+S_{-}\left(\omega \right)\\ =\left(\begin{array}{cc}0& E\\ -E^{*}& Q_1+\left(1+\omega \right)C\end{array}\right)+\left(\begin{array}{cc}B& 0\\ 0& -Q_1-\omega C\end{array}\right)=S_{+}\left(\omega \right)+H_{-}\left(\omega \right),\end{array}$ (2)

where $Q_1$ is a Hermitian positive semidefinite matrix and $\omega$ is a non-negative constant. Bai [2] extended the RHSS iteration method for the standard saddle-point problem to the stabilized saddle-point problem as follows:

$\{\begin{array}{l}\left(\alpha I+H_{+}\left(\omega \right)\right)x^{\left(k+\frac{1}{2}\right)}=\left(\alpha I-S_{-}\left(\omega \right)\right)x^{\left(k\right)}+b,\\ \left(\alpha I+S_{+}\left(\omega \right)\right)x^{\left(k+1\right)}=\left(\alpha I-H_{-}\left(\omega \right)\right)x^{\left(k+\frac{1}{2}\right)}+b,\end{array}$ (3)

where I is the identity matrix and $\alpha$ is a positive constant. And Bai proved the RHSS iteration method converges unconditionally to the exact solution of the stabilized saddle-point problem (1). Numerical results on stabilized saddle-point problems in [2] show that the RHSS iteration method significantly outperforms the Hermitian and skew-Hermitian splitting iteration method in iteration counts and computing times.

In this paper, inspired by Behzadi and Abdollahi [13] and Bai [2], we know that the accelerated RHSS iteration method has advantages over the HSS iteration method and the RHSS iteration method on the standard saddle-point problem. So we extend the accelerated RHSS iteration method for the standard saddle-point problem to the stabilized saddle-point problem, and establish the ARHSS iteration method for the stabilized saddle-point problem. Then we confirm the convergence of the ARHSS iteration method. Finally, experimental results show that the ARHSS iteration method is more effective than the HSS iteration method and the RHSS iteration method for the stabilized saddle-point problem (1).

The remainder of this work is organized as follows. In Section 2, we propose the ARHSS iteration method for the stabilized saddle-point problem (1). In Section 3, we give the convergence property of the ARHSS iteration method. In Section 4, A numerical experiment is given to show the effectiveness of the method. In Section 5, we end the paper with some brief conclusions.

2. The ARHSS Iteration Method

In this section, we derive the ARHSS iteration method for solving (1). Let $Q\in {ℂ}^{q\times q}$ be a Hermitian positive semidefinite matrix and $\omega$ be a prescribed non-negative parameter. Then we can split the stabilized saddle-point matrix $A\in {ℂ}^{n\times n}$ in (1), and obtain the regularized Hermitian and skew-Hermitian (RHS) splitting:

$\begin{array}{c}A=\left(\begin{array}{cc}B& 0\\ 0& Q+\omega C\end{array}\right)+\left(\begin{array}{cc}0& E\\ -E^{*}& -Q+\left(1-\omega \right)C\end{array}\right)=H_{+}+S_{-}\\ =\left(\begin{array}{cc}0& E\\ -E^{*}& Q+\left(1+\omega \right)C\end{array}\right)+\left(\begin{array}{cc}B& 0\\ 0& -Q-\omega C\end{array}\right)=S_{+}+H_{-}.\end{array}$ (4)

By applying the regularized Hermitian and Skew-Hermitian (RHS) splitting to (4), we then obtain the iteration scheme

$\{\begin{array}{l}\left(\Lambda +H_{+}\right)x=\left(\Lambda -S_{-}\right)x+b\\ \left(\Lambda +S_{+}\right)x=\left(\Lambda -H_{-}\right)x+b\end{array}$

where

$\Lambda =\left(\begin{array}{cc}\alpha I& 0\\ 0& \beta I\end{array}\right)$

with $\alpha$ and $\beta$ positive constants.

By iterating alternatively between these two fixed-point systems as

$\left(\Lambda +H_{+}\right)x^{\left(k+1/2\right)}=\left(\Lambda -S_{-}\right)x^{\left(k\right)}+b\mathrm{,}$

and

$\left(\Lambda +S_{+}\right)x^{\left(k+1\right)}=\left(\Lambda -H_{-}\right)x^{\left(k+1/2\right)}+b\mathrm{,}$

or in their blockwise forms

$\begin{array}{l}\left(\begin{array}{cc}\alpha I+B& 0\\ 0& \beta I+Q+\omega C\end{array}\right)\left(\begin{array}{c}{y}^{\left(k+1/2\right)}\\ z^{\left(k+1/2\right)}\end{array}\right)\\ =\left(\begin{array}{cc}\alpha I& -E\\ E^{\mathrm{<mo>\; *\; </mo>}}& \beta I+Q-\left(1-\omega \right)C\end{array}\right)\left(\begin{array}{c}{y}^{\left(k\right)}\\ z^{\left(k\right)}\end{array}\right)+\left(\begin{array}{c}f\\ g\end{array}\right)\mathrm{,}\end{array}$ (5)

and

$\begin{array}{l}\left(\begin{array}{cc}\alpha I& E\\ -E^{\mathrm{<mo>\; *\; </mo>}}& \beta I+Q+\left(1+\omega \right)C\end{array}\right)\left(\begin{array}{c}{y}^{\left(k+1\right)}\\ z^{\left(k+1\right)}\end{array}\right)\\ =\left(\begin{array}{cc}\alpha I-B& 0\\ 0& \beta I+Q+\omega C\end{array}\right)\left(\begin{array}{c}{y}^{\left(k+1/2\right)}\\ z^{\left(k+1/2\right)}\end{array}\right)+\left(\begin{array}{c}f\\ g\end{array}\right)\mathrm{.}\end{array}$ (6)

So, we obtain the ARHSS iteration method for solving the stabilized saddle-point problem (1) as follows:

The ARHSS iteration method. Let $\alpha$ and $\beta$ be positive constants, $\omega$ be a non-negative constant and $Q\in {ℂ}^{p\times p}$ be a Hermitian positive semidefinite matrix. Give an initial guess $x^{\left(0\right)}={\left(y^{{\left(0\right)}^{\ast }}\mathrm{,}{z}^{{\left(0\right)}^{\ast }}\right)}^{\ast }\in {ℂ}^n$, for $k=0,1,2,\cdots$ until the iteration sequence $\left\{x^{\left(k\right)}\right\}=\left\{{\left(y^{{\left(k\right)}^{\ast }}\mathrm{,}{z}^{{\left(k\right)}^{\ast }}\right)}^{\ast }\right\}\subset {ℂ}^n$ converges, compute the next iterate $x^{\left(k+1\right)}={\left(y^{{\left(k+1\right)}^{\ast }}\mathrm{,}{z}^{{\left(k+1\right)}^{\ast }}\right)}^{\ast }\in {ℂ}^n$ according to following procedure:

(i) solve for $y^{\left(k+1/2\right)}\in {ℂ}^p$ from the linear subsystem

$\left(\alpha I+B\right)y^{\left(k+1/2\right)}=\alpha I{y}^{\left(k\right)}-E{z}^{\left(k\right)}+f;$

(ii) compute

$f^{\left(k+1/2\right)}=\left(\alpha I-B\right)y^{\left(k+1/2\right)}+f\mathrm{<mo>\; ;\; </mo>}$

and

$g^{\left(k+1/2\right)}=\left(\beta I+Q+\left(\omega -1\right)C\right)z^{\left(k\right)}+E^{\mathrm{<mo>\; *\; </mo>}}{y}^{\left(k\right)}+2g\mathrm{<mo>\; ;\; </mo>}$

(iii) solve for $z^{\left(k+1\right)}\in {ℂ}^q$ from the linear subsystem

$\left(\beta I+Q+\left(1+\omega \right)C+\frac{1}{\alpha }{E}^{\mathrm{<mo>\; *\; </mo>}}E\right)z^{\left(k+1\right)}=\frac{1}{\alpha }{E}^{\mathrm{<mo>\; *\; </mo>}}{f}^{\left(k+1/2\right)}+g^{\left(k+1/2\right)}\mathrm{<mo>\; ;\; </mo>}$

(iv) compute

$y^{\left(k+1\right)}=\frac{1}{\alpha }\left(-E{z}^{\left(k+1\right)}+f^{\left(k+1/2\right)}\right)\mathrm{.}$

3. The Convergence Property of the ARHSS Iteration Method

In this section we prove the unconditional convergence of the ARHSS iteration method. We know that calculating the spectral radius of the iterative matrix $L\left(\alpha \mathrm{,}\beta \mathrm{,}\omega \right)$ is the easiest way to prove the convergence of ARHSS iteration method. From (5) and (6), we easily know $L\left(\alpha \mathrm{,}\beta \mathrm{,}\omega \right)$ can be equivalently rewritten as

$L\left(\alpha \mathrm{,}\beta \mathrm{,}\omega \right)=M{\left(\alpha \mathrm{,}\beta \mathrm{,}\omega \right)}^{-1}N\left(\alpha \mathrm{,}\beta \mathrm{,}\omega \right)\mathrm{,}$

where

$M\left(\alpha \mathrm{,}\beta \mathrm{,}\omega \right)=\frac{1}{2}\left(\begin{array}{cc}\alpha I+B& \frac{1}{\alpha }\left(\alpha I+B\right)E\\ -E^{\mathrm{<mo>\; *\; </mo>}}& \beta I+Q+\left(\omega +1\right)C\end{array}\right)$

and

$N\left(\alpha \mathrm{,}\beta \mathrm{,}\omega \right)=\frac{1}{2}\left(\begin{array}{cc}\alpha I-B& -\frac{1}{\alpha }\left(\alpha I-B\right)E\\ E^{\mathrm{<mo>\; *\; </mo>}}& \beta I+Q+\left(\omega -1\right)C\end{array}\right)$

To prove the convergence of the ARHSS iteration method, we first introduce a lemma.

Lemma 3.1. ( [2], Theorem 3.1) For the stabilized saddle-point problem (1), the RHSS iteration method (3) converges unconditionally to the exact solution of the stabilized saddle-point problem (1).

In the following theorem, we prove the unconditional convergence of the ARHSS iteration method.

Theorem 3.2. For the stabilized saddle-point problem (1), assume that $B\in {ℂ}^{p\times p}$ is Hermitian positive definite, $C\in {ℂ}^{q\times q}$ is Hermitian positive semidefinite and $E\in {ℂ}^{p\times q}$ satisfies $null\left(C\right)\cap null\left(E\right)=0$. Let $\alpha ,\beta >0$ be prescribed positive constant, $\omega$ be a given non-negative constant and $Q\in {ℂ}^{q\times q}$ be a given Hermitian positive semidefinite matrix. Then it holds that $\rho \left(L\left(\alpha \mathrm{,}\beta \mathrm{,}\omega \right)\right)<1$ , i.e., the ARHSS iteration method converges unconditionally to the exact solution of the stabilized saddle-point problem (1).

Proof. There are two cases to consider.

Case I: If $\beta \ge \alpha$, according to the ARHS splitting, we obtain

$\begin{array}{c}A=\left(\Lambda +H_{+}\right)-\left(\Lambda -S_{-}\right)\\ =\left(\alpha I+\left(\begin{array}{cc}B& 0\\ 0& \tilde{Q}+\omega C\end{array}\right)\right)-\left(\alpha I-\left(\begin{array}{cc}0& E\\ -E^{\mathrm{<mo>\; *\; </mo>}}& -\tilde{Q}+\left(1-\omega \right)C\end{array}\right)\right)\end{array}$ (7)

and

$\begin{array}{c}A=\left(\Lambda +S_{+}\right)-\left(\Lambda -H_{-}\right)\\ =\left(\alpha I+\left(\begin{array}{cc}0& E\\ -E^{*}& \tilde{Q}+\left(1+\omega \right)C\end{array}\right)\right)-\left(\alpha I-\left(\begin{array}{cc}B& 0\\ 0& -\tilde{Q}-\omega C\end{array}\right)\right)\end{array}$ (8)

where $\tilde{Q}=\left(\beta -\alpha \right)I+Q$. Since $\tilde{Q}$ is Hermitian positive semidefinite matrix, relations (7) and (8) are exactly the RHS splitting of A. By Lemma 3.1, this split converges to the exact solution of the saddle point linear system (1).

Case II: If $\beta <\alpha$, according to the ARHS splitting, we obtain

$\begin{array}{c}A=\left(\Lambda +H_{+}\right)-\left(\Lambda -S_{-}\right)\\ =\left(\beta I+\left(\begin{array}{cc}B& 0\\ 0& \tilde{Q}+\omega C\end{array}\right)\right)-\left(\beta I-\left(\begin{array}{cc}0& E\\ -E^{\mathrm{<mo>\; *\; </mo>}}& -\tilde{Q}+\left(1-\omega \right)C\end{array}\right)\right)\end{array}$ (9)

and

$\begin{array}{c}A=\left(\Lambda +S_{+}\right)-\left(\Lambda -H_{-}\right)\\ \mathrm{<mo>\; =\; </mo>}\left(\beta I+\left(\begin{array}{cc}0& E\\ -E^{\mathrm{<mo>\; *\; </mo>}}& \tilde{Q}+\left(1+\omega \right)C\end{array}\right)\right)-\left(\beta I-\left(\begin{array}{cc}B& 0\\ 0& -\tilde{Q}-\omega C\end{array}\right)\right)\end{array}$ (10)

where $\tilde{Q}=\left(\alpha -\beta \right)I+Q$. Since $\tilde{Q}$ is Hermitian positive semidefinite matrix, the relations (9) and (10) are exactly the RHS splitting of A. By Lemma 3.1, this split converges to the exact solution of the saddle point linear system (1). The desired result holds by Case I and Case II. Therefore, we complete the proof.

4. Numerical Experiment

In this section, a numerical example is given to illustrate the effectiveness of the ARHSS iteration method for the stabilized saddle-point problem (1). We compare the number of iteration steps (denoted as IT) and computing times in seconds (denoted as CPU) of the ARHSS iteration method with the HSS iteration method and the RHSS iteration method. All our tests are started from the initial vector $x^{\left(0\right)}=0$. All iterations are terminated as the current residuals satisfy $\Vert b-A{x}^{\left(k\right)}\Vert \le {10}^{-6}\times \Vert b\Vert$ or the number of iterations steps exceeds 5000. In addition, all the computations were implemented in MATLAB [version 9.1.0.441655 (R2016b)] in double precision on a personal computer with 1.60 GHZ central processing unit [Intel(R) Core(TM)i5-5250U] and 4.00 GB memory.

Example 4.1. ( [2], Example 4.3) Consider the nonlinear image restoration problem

$\min {\Vert \tilde{f}-s\left(Ky\right)\Vert }_2^2+\beta {\Vert y\Vert }_2^2\mathrm{,}$ (11)

where $\tilde{f}$ and $y\in {ℂ}^p$ represent the observed and the original images,respectively, $s\mathrm{<mo>\; :\; </mo>}{ℂ}^p\to ℂ$ denotes a nonlinear point spread function and $K=\left({\left[K\right]}_{ij}\right)\in {ℂ}^{p\times p}$ is a blurring matrix. When this regularized nonlinear least-squares problem (11) is solved by the Gauss-Newton iteration method,at each step for the currently available approximant $y_c=\left({\left[y_c\right]}_1\mathrm{,}{\left[y_c\right]}_2\mathrm{,}\cdots \mathrm{,}{\left[y_c\right]}_p\right)$ ,we need to solve a linear system of the form

$\left(\beta I+K^{⊺}{D}^{2}K\right)y=K^{⊺}D\left[\tilde{f}-s\left(K{y}_c\right)+DK{y}_c\right]$

to obtain the next approximant,where $D=diag\left({\left[D\right]}_1\mathrm{,}{\left[D\right]}_2\mathrm{,}\cdots \right)\mathrm{,}{\left[D\right]}_p\in {ℂ}^{p\times p}$ is a diagonal matrix with its diagonal entries being given by

${\left[D\right]}_i={\frac{\partial s}{\partial \xi }|}_{\xi ={\displaystyle {\sum }_{j=1}^p{\left[K\right]}_{ij}{\left[y_c\right]}_j}}\mathrm{,}i=\mathrm{1,2,}\cdots \mathrm{,}p\mathrm{.}$

This linear system can be equivalently reformulated into the stabilized saddle-point problem (1),in which

$B=D^{-2},E=K,C=\beta I$

and

$f=D^{-1}\left[\tilde{f}-s\left(K{y}_c\right)\right]+K{y}_c,g=0.$

In actual computations we set

$\tilde{f}={\left[\frac{254}{p}\mathrm{<mo>\; :\; </mo>}\frac{254}{p}\mathrm{<mo>\; :\; </mo>254}\right]}^{⊺}$

and

$y_c={\left[\left[0.5+\frac{508}{p}\mathrm{<mo>\; :\; </mo>}\frac{508}{p}\mathrm{<mo>\; :\; </mo>254.5}\right]\mathrm{,}\left[\mathrm{254.5\; <mo>\; :\; </mo>}-\frac{508}{p}\mathrm{<mo>\; :\; </mo>0.5}+\frac{508}{p}\right]\right]}^{⊺}$

and take $\beta ={10}^{-3}$ ,

$s\left(\xi \right)=30\log \left(\xi \right),{\left[K\right]}_{ij}=\frac{1}{\sqrt{2\pi }\mu }\exp \left(\frac{-{\left|i-j\right|}^2}{2{\mu }^2}\right),i,j=1,2,\cdots ,p,$

with $\mu =2$. Here,the notation “:”is a standard MATLAB symbol used in indicating a row vector of the form

$\left[x_f:\delta x:x_l\right]=\left[x_f,x_f+\delta x,\cdots ,x_f+k\delta x,\cdots ,x_l\right],$

for which $x_f$ and $x_l$ are the first and the last elements, $\delta x$ is a prescribed increment such that $x_l-x_f$ is its multiple and k is a non-negative integer ranging from 0 to $\frac{x_l-x_f}{\delta x}$.

For this example we have $p=q$, so that the dimension of the stabilized saddle-point matrix $A\in {\mathbb{B}}^{p\times p}$ is $n=2p$. Note that the Toeplitz matrix E is highly ill-conditioned with rapidly decaying singular values so that it is almost rank-deficient, especially when the problem size p becomes sufficiently large. Besides, we take the regularization matrix to be

$Q=\left(\alpha \gamma -\omega \right)C+\gamma E^{⊺}E\mathrm{.}$

Taking Q into (5) and (6), we get

$g^{\left(k+1/2\right)}=\left(\beta I+\left(\alpha \gamma -1\right)C+\gamma E^{\mathrm{<mo>\; *\; </mo>}}E\right)z^{\left(k\right)}+E^{\mathrm{<mo>\; *\; </mo>}}{y}^{\left(k\right)}+2g\mathrm{<mo>\; ;\; </mo>}$

and

$\left(\beta I+\left(\alpha \gamma +1\right)C+\frac{1}{\alpha }{E}^{\mathrm{<mo>\; *\; </mo>}}E\right)z^{\left(k+1\right)}=\frac{1}{\alpha }{E}^{\mathrm{<mo>\; *\; </mo>}}{f}^{\left(k+1/2\right)}+g^{\left(k+1/2\right)}\mathrm{.}$

So, in the specific algorithms, the parameter $\omega$ does not need to be determined. In our experiments, let the parameter $\gamma$ be taken as the optimal value of [ [2], Example 4.3], and the optimal values of the parameters $\alpha$ and $\beta$ can be determined by experience and trial runs. For example, first, suppose the ranges of $\alpha$ and $\beta$ are $\left(\mathrm{0,1}\right)$, and then divide it into 100 equal parts. Second, fix parameter $\alpha$ to a certain value in $\left(\mathrm{0,1}\right)$, and then find the corresponding optimal parameter $\beta$ by trial runs. Third, using the second step method, fix parameter $\beta$ to a certain value in $\left(\mathrm{0,1}\right)$, and then find the corresponding optimal parameter $\alpha$ by trial runs. Finally, we choose the best parameters $\alpha$ and $\beta$ through experience and trial runs.

In Table 1, we list the best parameters $\alpha \mathrm{,}\beta$ and regularization parameter $\gamma$ for the three iteration methods including HSS, RHSS and ARHSS with different mesh number p. Through these optimal parameters, we list the experimental results of the three methods in Table 2. The experimental results show that compared with HSS and RHSS, the ARHSS iteration method takes the least CPU and the smallest IT to converge to the solution of (11).

Table 1. Optimal parameter and regularization parameter.

Table 2. IT and CPU.

5. Conclusion

For stabilized saddle-point problem, we extend the accelerated RHSS iteration method for the standard saddle-point problem to the stabilized saddle-point problem. The convergence property of the ARHSS iteration method is carefully studied, and the numerical result illustrates that the ARHSS is more effective for the stabilized saddle-point problem. In this paper, optimal iteration parameters can only be determined experimentally. So, how to establish practical formulas for computing the optimal iteration parameters is worth further study.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Bai, Z.-Z. (2019) Regularized HSS iteration Methods for Stabilized Saddle-Point Problems. IMA Journal of Numerical Analysis, 39, 1888-1923. https://doi.org/10.1093/imanum/dry046
[2]	Wathen, A.J. and Silvester, D.J. (1993) Fast Iterative Solution of Stabilished Stokes Systems I: Using Simple Diagonal Preconditioners. SIAM Journal on Numerical Analysis, 30, 630-649. https://doi.org/10.1137/0730031
[3]	Silvester, D.J. and Kechkar, N. (1990) Stabilised Bilinear-Constant Velocity-Pressure Finite Elements for the Conjugate Gradient Solution of the Stokes Problems. Computer Methods in Applied Mechanics and Engineering, 79, 71-86. https://doi.org/10.1016/0045-7825(90)90095-4
[4]	Elman, H.C., Silvester, D.J. and Wathen, A.J. (2002) Performance and Analysis of Saddle Point Preconditioners for the Discrete Steady-State Navier-Stokes Equations. Numerische Mathematik, 90, 665-688. https://doi.org/10.1007/s002110100300
[5]	Benzi, M., Golub, G.H. and Liesen, J. (2005) Numerical Solution of Saddle Point Problems. Acta Numerica, 14, 1-137. https://doi.org/10.1017/S0962492904000212
[6]	Bai, Z.-Z. (2006) Structured Preconditioners for Nonsingular Matrices of Block Two-by-Two Structures. Mathematics of Computation, 75, 791-815. https://doi.org/10.1090/S0025-5718-05-01801-6
[7]	Wathen, A.J. (2015) Preconditioning. Acta Numerica, 124, 329-376. https://doi.org/10.1017/S0962492915000021
[8]	Benzi, M. and Golub, G.H. (2004) A Preconditioner for Generalized Saddle Point Problems. SIAM Journal on Matrix Analysis and Applications, 26, 20-41. https://doi.org/10.1137/S0895479802417106
[9]	Bai, Z.-Z., Golub, G.H. and Ng, M.K. (2003) Hermitian and Skew-Hermitian Splitting Methods for Non-Hermitian Positive Definite Linear Systems. SIAM Journal on Matrix Analysis and Applications, 24, 603-626. https://doi.org/10.1137/S0895479801395458
[10]	Bai, Z.-Z., Golub, G.H. and Pan, J.-Y. (2004) Preconditioned Hermitian and Skew-Hermitian Splitting Methods for Non-Hermitian Positive Semidefinite Linear Systems. Numerische Mathematik, 98, 1-32. https://doi.org/10.1007/s00211-004-0521-1
[11]	Bai, Z.-Z. and Golub, G.H. (2007) Accelerated Hermitian and Skew-Hermitian Splitting Iteration Methods for Saddle-Point Problems. IMA Journal of Numerical Analysis, 27, 1-23. https://doi.org/10.1093/imanum/drl017
[12]	Bai, Z.-Z. and Benzi, M. (2017) Regularized HSS Iteration Methods for Saddle-Point Linear Systems. BIT Numerical Mathematics, 57, 287-311. https://doi.org/10.1007/s10543-016-0636-7
[13]	Behzadi, R. and Abdollahi, F. (2018) Accelerated RHSS Iteration Methods for Saddle-Point Linear Systems. UPB Scientific Bulletin, Series A, 80, 153-162.