# GMRES algorithms over 35 years

Qinmeng Zou School of Science, Beijing University of Posts and
Telecommunications, Beijing 100876, China<EMAIL_ADDRESS>
###### Abstract

This paper is about GMRES algorithms for the solution of nonsingular linear
systems. We first consider basic algorithms and study their convergence. We
then focus on acceleration strategies and parallel algorithms that are useful
for solving challenging systems. We also briefly discuss other problems, such
as systems with multiple right-hand sides, shifted systems, and singular
systems.

Keywords. linear systems; Gram-Schmidt; mixed precision; parallel computing;
preconditioning.

##  1 Introduction

Generalized minimal residual (GMRES) algorithms are widely used for solving
nonsymmetric linear systems arising from partial differential equations. They
possess some optimality property and behave reasonably well in practice. The
first and most well-known algorithm is that of Saad and Schultz, which was
initially introduced in a technical report in 1983, and formally published in
1986 [223]. Since then, numerous variants appeared, as well as studies of
their convergence and accuracy. In this paper we concentrate on GMRES
algorithms. We sketch the main developments and focus on the algorithmic
innovations in the past 35 years.

Consider the nonsingular linear system

$Ax=b,$ (1)

where $A\in\mathbb{R}^{N\times N}$ and $b\in\mathbb{R}^{N}$. Given an initial
approximation $x_{0}$, a projection method constructs a sequence of
approximate solutions $x_{n}$, $n=1,2,\dots$, such that

$x_{n}\in x_{0}+\mathcal{S}_{n},\quad r_{n}\perp\mathcal{C}_{n},$

where $r_{n}=b-Ax_{n}$ is the residual vector and $\mathcal{S}_{n}$ and
$\mathcal{C}_{n}$ are $n$-dimensional subspaces, called respectively search
space and constraint space. As $A$ is nonsingular, $x_{n}$ is uniquely
defined. In general, we build a sequence of nested search spaces

$\mathcal{S}_{1}\subset\mathcal{S}_{2}\subset\mathcal{S}_{3}\subset\dots$ (2)

to ensure finite termination; see, e.g., [167]. A minimal residual method
obtains the optimal approximation by minimizing the residual norm over all
candidate vectors in the search space

$\left\lVert r_{n}\right\rVert=\left\lVert b-Ax_{n}\right\rVert=\min_{x\in
x_{0}+\mathcal{S}_{n}}\left\lVert b-Ax\right\rVert,$ (3)

where $\left\lVert\cdot\right\rVert$ denotes the 2-norm or the corresponding
induced matrix norm. This class of methods can be interpreted as a projection
process with $\mathcal{C}_{n}=A\mathcal{S}_{n}$. In other words, they find an
optimal correction in $\mathcal{S}_{n}$, such that $r_{0}-r_{n}$ is the
orthogonal projection of $r_{0}$ onto $A\mathcal{S}_{n}$; see [81, 219, 80]
for further analysis. The Krylov subspace is the most broadly employed search
space and is defined by

$\mathcal{K}_{n}=\mathcal{K}_{n}(A,r_{0})=\operatorname{span}\\{r_{0},Ar_{0},\dots,A^{n-1}r_{0}\\}.$
(4)

These subspaces are nested in the sense of (2) and can be built up gradually
using only matrix-vector multiplication by $A$. The residual $r_{n}$ lies in
the next Krylov subspace $\mathcal{K}_{n+1}$. If $d$ denotes the grade of
$r_{0}$ with respect to $A$, then $A\mathcal{K}_{d}\subseteq\mathcal{K}_{d}$
holds. Therefore, the sequence of Krylov subspaces will eventually become
invariant. A projection method with $\mathcal{S}_{n}=\mathcal{K}_{n}$ is
called a Krylov subspace method. Practically a Krylov subspace method is
terminated as soon as the approximate solution is good enough. Ways to
formulate typical Krylov subspace methods and derivations of their properties
have been extensively discussed in the literature; see, e.g., [155, 105, 121,
220, 266, 243, 167, 186] and references therein. A comparison of Krylov
subspace methods can be found in a recent work [112]; see also [167, 222] for
a historical perspective.

Classical GMRES algorithms realize the so-called minimal residual Krylov
subspace method, or loosely called GMRES method, which is a projection process
satisfying (3) and taking (4) as search space. In this case, the projection
process can be rewritten as follows:

$\displaystyle x_{n}\in x_{0}+\mathcal{K}_{n},$ (5a) $\displaystyle r_{n}\perp
A\mathcal{K}_{n},$ (5b)

or equivalently,

$\left\lVert r_{n}\right\rVert=\left\lVert b-Ax_{n}\right\rVert=\min_{x\in
x_{0}+\mathcal{K}_{n}}\left\lVert b-Ax\right\rVert.$ (6)

We shall consider both sequential and parallel GMRES algorithms. Some variants
such as restarted or hybrid algorithms exhibit more complex behavior and may
no longer possess the finite termination property. Other crucial aspects that
can dramatically affect the performance of Krylov subspace algorithms include
the finite precision effect and communication costs; see [54] for a general
insight on the cost of the iterative computation and its role in Krylov
subspace methods. In what follows these topics shall be discussed to some
extent. We do not talk here about general-purpose preconditioning techniques
like incomplete factorization and algebraic multigrid, since these topics
deserve individual coverage by themselves, for which we refer the reader to
[27, 281]; see also [112] for experiments. For the similar reason, problem-
specific techniques will also not be covered, though they can be very
successful when prior knowledge of the source problem is available; see, e.g.,
[201, 109, 208]. We only mention certain techniques that are closely related
to GMRES iterations. Throughout most of the paper we restrict to the
nonsingular linear system (1). In Section 5, however, we briefly mention some
work where this restriction is lifted.

The paper is organized as follows. We begin, in Section 2, by formulating
basic GMRES algorithms, giving some equivalent formulations, and presenting
theoretical results in both exact and finite arithmetic. In Section 3 we
summarize various strategies for improving classical algorithms. In Section 4
we discuss parallel techniques that can reduce or hide communication overhead.
Sections 5 is devoted to a brief literature review of related topics. Finally,
concluding remarks are draw in Section 6.

##  2 Basic algorithms and convergence

The techniques introduced in this section form the basis for subsequent
discussions. In particular, the algorithm developed by Saad and Schultz [223]
is still in common use today. Methodologies have been developed to partially
characterize their convergence behavior. Note that the literature on Krylov
subspace methods uses the term “GMRES” almost as a synonym for the specific
algorithm developed in [223]. For the sake of clarity, here the latter is
referred to as “MGS-GMRES”, as the modified Gram-Schmidt process is applied in
[223] to orthogonalize successive basis vectors of (4). Later, the distinction
between the two terms shall often be blurred because in the literature GMRES
is almost exclusively implemented by means of a Gram-Schmidt-like process.

###  2.1 MGS-GMRES

Saad and Schultz proposed the first GMRES algorithm in [223], which uses the
Arnoldi process [9] to generate an orthonormal basis $\\{v_{1},\dots,v_{n}\\}$
of the Krylov subspace (4). The key relation used in the Arnoldi process is

$h_{j+1,j}v_{j+1}=Av_{j}-\sum_{i=1}^{j}h_{i,j}v_{i},$ (7)

where $v_{1}=r_{0}/\left\lVert r_{0}\right\rVert$ and $h_{i,j}$ are selected
such that $v_{j+1}$ is normalized and orthogonal to all the previous basis
vectors, that is,

$\displaystyle h_{i,j}=(Av_{j},v_{i}),\quad i=1,\dots,j,$ (8) $\displaystyle
h_{j+1,j}=\lVert Av_{j}-\sum_{i=1}^{j}h_{i,j}v_{i}\rVert,$ (9)

where $(\cdot,\cdot)$ denotes the Euclidean inner product. The above rules
correspond directly to the classical Gram-Schmidt (CGS) orthogonalization
process, leading to the so-called CGS-Arnoldi algorithm. In practice the
computation of $v_{j+1}$ by (7)–(9) undergoes a severe cancellation. Special
precautions, such as the reorthogonalization [61, 138] and modified Gram-
Schmidt (MGS) processes, have been recommended to ensure numerical stability.
The former is costly and generally used in the context of parallel algorithms
[141, 256], eigenvalue problems [221], and mixed precision [172], while MGS
has the same operation count as CGS, but with better numerical properties; see
[162] for a thorough presentation of Gram-Schmidt orthogonalization. Arnoldi
with MGS orthogonalization (MGS-Arnoldi) [215] can be viewed as the standard
version of Arnoldi iteration.

Algorithm 1 MGS-Arnoldi

1:$v_{1}=r_{0}/\left\lVert r_{0}\right\rVert$

2:for $j=1,2,\dots,n$ do

3: $w_{j}=Av_{j}$

4: for $i=1,2,\dots,j$ do

5: $h_{i,j}=(w_{j},v_{i}),\quad w_{j}=w_{j}-h_{i,j}v_{i}$

6: end for

7: $h_{j+1,j}=\left\lVert w_{j}\right\rVert$, $v_{j+1}=w_{j}/h_{j+1,j}$

8:end for

In Algorithm 1, we can see that lines 3–6 are mathematically equivalent to
(7)–(8). Let $\hat{P}_{j}$ denote the MGS projector. Then, according to
Algorithm 1, one finds
$\hat{P}_{j}=(I-v_{j}v_{j}^{\intercal})\dots(I-v_{1}v_{1}^{\intercal})$.
Recalling the grade $d$ of $r_{0}$ defined in the preceding section, we
observe that $h_{j+1,j}=0$ (line 7) if and only if $j=d$, that is, the maximal
dimension of the Krylov subspace is attained. In what follows we shall assume
$n<d$.

Let us write

$V_{n}=[v_{1},\dots,v_{n}]\in\mathbb{R}^{N\times n},\quad
H_{n}=[h_{i,j}]\in\mathbb{R}^{n\times n},\quad\mkern
3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern
1.0mu_{n}=\left[\begin{array}[]{c}H_{n}\\\
h_{n+1,n}e_{n}^{\intercal}\end{array}\right]\in\mathbb{R}^{(n+1)\times n},$

where $e_{i}$ denotes the $i$th column of the identity matrix of appropriate
order and $h_{i,j}=0$ for $i>j+1$. Therefore, $\mkern
3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern 1.0mu_{n}$ and $H_{n}$ are
upper Hessenberg matrices. It follows that

$\displaystyle AV_{n}$ $\displaystyle=V_{n+1}\mkern
3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern 1.0mu_{n}$ (10a)
$\displaystyle=V_{n}H_{n}+h_{n+1,n}v_{n+1}e_{n}^{\intercal},$ (10b)
$\displaystyle V_{n}^{\mathsf{T}}AV_{n}$ $\displaystyle=H_{n},$ (10c)

which can be interpreted as an orthogonal reduction of $A$ to upper Hessenberg
form. Let $y$ denote the coordinate vector of the correction $x_{n}-x_{0}$ in
the basis $V_{n}$. Then it follows from the minimal residual process (5)–(6)
that

$\displaystyle y$ $\displaystyle=\operatorname*{arg\,min}_{\hat{y}}\left\lVert
r_{0}-AV_{n}\hat{y}\right\rVert=\operatorname*{arg\,min}_{\hat{y}}\left\lVert\beta
e_{1}-\mkern 3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern
1.0mu_{n}\hat{y}\right\rVert,$ (11) $\displaystyle x_{n}$
$\displaystyle=x_{0}+V_{n}y,$

where $\beta=\left\lVert r_{0}\right\rVert$. Saad and Schultz [223] suggested
using Givens rotations that reduce the upper Hessenberg matrix $\mkern
3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern 1.0mu_{n}$ recursively to the
upper triangular matrix.

Algorithm 2 MGS-GMRES

1:$\beta=\left\lVert r_{0}\right\rVert,\quad v_{1}=r_{0}/\beta,\quad g=\beta
e_{1}$

2:for $j=1,2,\dots,n$ do

3: $w_{j}=Av_{j}$

4: for $i=1,2,\dots,j$ do

5: $h_{i,j}=(w_{j},v_{i}),\quad w_{j}=w_{j}-h_{i,j}v_{i}$

6: end for

7: $h_{j+1,j}=\left\lVert w_{j}\right\rVert,\quad v_{j+1}=w_{j}/h_{j+1,j}$

8: for $i=1,2,\dots,j-1$ do

9: $\gamma=c_{i}h_{i,j}+s_{i}h_{i+1,j},\quad
h_{i+1,j}=-s_{i}h_{i,j}+c_{i}h_{i+1,j},\quad h_{i,j}=\gamma$

10: end for

11: $\delta=\sqrt{h_{j,j}^{2}+h_{j+1,j}^{2}},\quad
s_{j}=h_{j+1,j}/\delta,\quad c_{j}=h_{j,j}/\delta$

12: $h_{j,j}=c_{j}h_{j,j}+s_{j}h_{j+1,j},\quad h_{j+1,j}=0$

13: $g_{j+1}=-s_{j}g_{j},\quad g_{j}=c_{j}g_{j},\quad\rho=\left\lvert
g_{j+1}\right\rvert$, if $\rho$ small enough, then $n=j$, go to 15

14:end for

15:Solve upper triangular linear system $H_{n}y=g_{1:n}$

16:$x_{n}=x_{0}+V_{n}y$

We show MGS-GMRES as Algorithm 2, where $\rho$ denotes the residual norm,
$g_{i}$ denotes the $i$th entry of vector $g$ and $g_{i:j}$ refers to the
subvector with elements indexed by $i$ through $j$. Lines 3–7 are the same as
in Algorithm 1, followed by Givens rotations

$\left(\begin{array}[]{c}h_{i,j}\\\
h_{i+1,j}\end{array}\right)\leftarrow\left(\begin{array}[]{cc}c_{i}&s_{i}\\\
-s_{i}&c_{i}\end{array}\right)\left(\begin{array}[]{c}h_{i,j}\\\
h_{i+1,j}\end{array}\right)$

to maintain a QR factorization of the upper Hessenberg matrix $\mkern
3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern 1.0mu_{j}$, which involves all
the previous rotations (lines 8–10) and a new rotation (lines 11–13) to
annihilate the subdiagonal entry. In particular, $s_{j}$ and $c_{j}$
correspond to the sine and cosine of the $j$th rotation angle such that
$h_{j+1,j}=0$; we refer the reader to [31] for a more robust implementation of
Givens rotations. An important by-product is that the residual norm
$\left\lvert g_{j+1}\right\rvert$ can be readily obtained at the end of each
step (line 13). Note that the CGS-based GMRES (CGS-GMRES) algorithm can be
easily derived by moving the vector update operations in line 5 out of the
for-loop. The CGS projector can be written as
$\tilde{P}_{j}=I-V_{j}V_{j}^{\mathsf{T}}$.

In [223], it was proved that MGS-GMRES never breaks down. Besides, in a
remarkable work by Paige et al. [204], it was proved that MGS-GMRES is
backward stable. By modifying the constraint (5b) as
$r_{n}\perp\mathcal{K}_{n}$, we can get the so-called FOM method [216]. A
detailed account of relations between GMRES and FOM can be found in [38, 59];
see also [220]. FOM has received far less attention partly because the
corresponding algorithms can break down. However, it was shown in [38] that
for the cases FOM performs poorly, nor will GMRES offers much help. In this
situation, the latter may still be preferred due to its minimal residual
property. GMRES has also been compared with CMRH and Lanczos-based methods;
see [227, 229, 140].

Algorithm 2 uses long recurrences to orthogonalize $Av_{j}$ against all the
previous basis vectors, the cost of which grows quadratically with the number
of steps. Saad and Schultz [223] considered restarting and truncation as
remedies. The truncation strategy keeps only a small number of previous basis
vectors in the orthogonalization phase, and computes the approximate solution
in a progressive manner; see [224]. In the restarted version, a parameter is
introduced to limit the size of the Krylov subspaces. The resulting scheme
outlined in Algorithm 3 is called restarted GMRES, denoted by GMRES($m$).
Practically we use the restarted scheme on the top level, and a specific GMRES
algorithm at the lower level; see [104] for an efficient implementation with
variable orthogonalization schemes.

Algorithm 3 GMRES($m$)

1:for $k=0,1,\dots$ do

2: Compute an $m$-dimensional basis $Z_{m}$ and the coordinate vector $y$
using a GMRES algorithm

3: $x_{k+1}=x_{k}+Z_{m}y$

4: Compute the residual norm $\rho$, if $\rho$ small enough, then stop

5:end for

Unfortunately, the desired properties such as finite termination and residual
minimization are lost. Both strategies may result in poor convergence or even
stagnation; see, e.g., [99, 251, 65, 17, 115, 249]. Moreover, using a large
$m$ may not necessarily lead to fast convergence (see [81, 90]). We shall
discuss other potential acceleration strategies in Section 3.

###  2.2 HH-GMRES

The minimal residual Krylov subspace method can be completely described by (5)
from a mathematical point of view. From a computational point of view,
however, mathematically equivalent algorithms may have quite different
numerical behavior. Since at that time the backward stability of MGS-GMRES
[204] was still unknown, motivated by numerical stability concerns, Walker
[277] in 1988 proposed an algorithm that uses Householder reflections to
orthogonalize the basis vectors; see also [278].

Algorithm 4 HH-GMRES

1:$w_{1}=r_{0}+\operatorname{sign}(r_{0,1})\left\lVert r_{0}\right\rVert
e_{1}$, $P_{1}=I-2w_{1}w_{1}^{\intercal}/(w_{1}^{\intercal}w_{1})$,
$\beta=-\operatorname{sign}(r_{0,1})\left\lVert r_{0}\right\rVert$,
$v_{1}=P_{1}e_{1}$

2:for $j=1,2,\dots,n$ do

3: $u=P_{j}P_{j-1}\dots P_{1}Av_{j}$,
$w_{j+1}=[0,\dots,0,u_{j+1},\dots,u_{N}]^{\mathsf{T}}$

4: $w_{j+1}=w_{j+1}+\operatorname{sign}(u_{j+1})\left\lVert
w_{j+1}\right\rVert e_{j+1}$,
$P_{j+1}=I-2w_{j+1}w_{j+1}^{\intercal}/(w_{j+1}^{\intercal}w_{j+1})$

5: $h_{j}=P_{j+1}u$, $v_{j+1}=P_{1}P_{2}\dots P_{j+1}e_{j+1}$

6:end for

7:Define $\mkern 3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern 1.0mu_{n}$ as
the first $n+1$ rows of the matrix $[h_{1},\dots,h_{n}]$

8:Compute $y$ such that $\left\lVert\beta e_{1}-\mkern
3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern 1.0mu_{n}y\right\rVert$ is
minimized (see lines 8–15 of Algorithm 2)

9:$x_{n}=x_{0}+V_{n}y$

Algorithm 4 depicts the Householder-based GMRES (HH-GMRES) algorithm. Here,
$\operatorname{sign}(\alpha)$ returns $1$ or $-1$ depending on the sign of
$\alpha$ and $r_{i,j}$ denotes the $j$th element of $r_{i}$. The signs (lines
1 and 4) are chosen to reduce the risk of subtractive cancellation. Line 4
corresponds directly to Householder reflections; the details can be found in
[118, 220]. HH-GMRES gradually generates a QR factorization

$[r_{0},Av_{1},\dots,Av_{n}]=P_{1}P_{2}\dots P_{n+1}[\beta
e_{1},h_{1},\dots,h_{n}].$

By discarding the zero rows of the upper triangular matrix and the first
column of both sides, we can obtain the Arnoldi relation (10). The Householder
vector $w_{j+1}$ contains $j$ zeros in the first rows (line 3) such that the
Householder transformation can leave the previous vectors
$h_{1},\dots,h_{j-1}$ unchanged. For the solution of line 8 and the estimate
of the residual norm at each iteration step, the same techniques used in
Algorithm 2 can be applied. Note that we only need to store the $w_{i}$’s
instead of the $P_{i}$’s. We also prefer not to explicitly store the $v_{i}$’s
but use the Qin Jiushao’s scheme (Horner’s scheme) to evaluate line 9 in a
recursive manner; see, e.g., [220].

Similar to MGS-GMRES, the cost of HH-GMRES increases dramatically as the
number of steps increases. Therefore, the restarted scheme in Algorithm 3 is
commonly employed as an outer loop. HH-GMRES was proved earlier to be backward
stable (see [73]), but can be roughly two times more expensive than MGS-GMRES.

###  2.3 Simpler GMRES

In 1994, Walker and Zhou [280] proposed another variant of the GMRES method,
called simpler GMRES (SGMRES), which does not require the QR factorization of
an upper Hessenberg matrix. Although rarely used in practice due to stability
issues, this algorithm can enhance our understanding of minimal residual
iterations.

We present a generalized version as given in [150] that contains the essential
features of SGMRES. Consider a nonorthogonal basis
$Z_{n}=\\{z_{1},\dots,z_{n}\\}$ of $\mathcal{K}_{n}$ and an orthonormal basis
$V_{n}=\\{v_{1}\dots,v_{n}\\}$ of $A\mathcal{K}_{n}$. We construct a matrix
relation

$AZ_{n}=V_{n}T_{n}$ (12)

such that $T_{n}=[t_{i,j}]\in\mathbb{R}^{n\times n}$ is an upper triangular
matrix. According to the constraint (5b), one finds
$r_{n}=(I-V_{n}V_{n}^{\mathsf{T}})r_{0}$, yielding the following update:

$\alpha_{n}=(r_{n-1},v_{n}),\quad r_{n}=r_{n-1}-\alpha_{n}v_{n}.$ (13)

In addition, since $x_{n}=x_{0}+Z_{n}y$ with some n-dimensional vector $y$, it
follows that $r_{n}=r_{0}-V_{n}T_{n}y$. Combining (5b) and (13), one obtains
that

$T_{n}y=V_{n}^{\mathsf{T}}r_{0}=[\alpha_{1},\dots,\alpha_{n}]^{\mathsf{T}}.$

Algorithm 5 Generalized Simpler GMRES

1:for $j=1,2,\dots,n$ do

2: Choose $z_{j}$ and compute $w_{j}=Az_{j}$

3: Orthonormalize $w_{j}$ against $v_{1},\dots,v_{j-1}$ to obtain
$t_{1,j},\dots,t_{j,j}$ and $v_{j}$

4: $\alpha_{j}=(r_{j-1},v_{j}),\quad r_{j}=r_{j-1}-\alpha_{j}v_{j}$

5: Compute the residual norm $\rho$, if $\rho$ small enough, then $n=j$, go to
7

6:end for

7:Solve $T_{n}y=[\alpha_{1},\dots,\alpha_{n}]^{\mathsf{T}}$

8:$x_{n}=x_{0}+Z_{n}y$

The generalized simpler GMRES framework (see [150]) is listed in Algorithm 5.
Choosing $Z_{n}=[r_{0}/\left\lVert r_{0}\right\rVert,V_{n-1}]$ results in
SGMRES, which was formalized by virtue of both MGS- and Householder-based
processes in [280]. Note that in the Householder variant there is no need to
store the $v_{i}$’s and thus line 4 should be modified accordingly; see [280]
for details. The choice $Z_{n}=[r_{0}/\left\lVert
r_{0}\right\rVert,\dots,r_{n-1}/\left\lVert r_{n-1}\right\rVert]$ was
introduced in [150], leading to a slightly different algorithm called
residual-based simpler GMRES (RB-SGMRES). Then, an adaptive variant combining
the two choices was proposed in [149], following the rule

$z_{n}=\begin{cases}\frac{r_{n-1}}{\left\lVert
r_{n-1}\right\rVert},&\left\lVert r_{n-1}\right\rVert\leq\omega\left\lVert
r_{n-2}\right\rVert,\\\ v_{n-1},&\text{otherwise},\end{cases}$ (14)

with $\omega\in[0,1]$. Notice that $\omega=0$ and $\omega=1$ correspond to
SGMRES and RB-SGMRES, respectively.

Theoretical aspects are collected in [280, 163, 164, 56, 150, 149], which
provides bounds for the condition number of $Z_{n}$ and allows for a better
understanding of the unstable behavior of Algorithm 5; see also [186]. It’s
interesting to observe that, along with the results in [128], SGMRES can be
employed to describe GMRES convergence; see [163].

###  2.4 Equivalent formulations

We briefly review some iterative schemes that are in some sense mathematically
equivalent to GMRES. Paige and Saunders [205] in 1975 proposed the MINRES
algorithm satisfying (6) for solving symmetric indefinite systems, which
motivated the development of MGS-GMRES. In 1982, Elman [88, 83] proposed the
GCR algorithm which amounted to generalizing the conjugate residual algorithm
(see, e.g., [220]) to the nonsymmetric case.

Algorithm 6 GCR

1:$q_{0}=r_{0}$

2:for $j=0,1,\dots,n$ do

3: $\alpha_{j}=(r_{j},Aq_{j})/(Aq_{j},Aq_{j})$

4: $x_{j+1}=x_{j}+\alpha_{j}q_{j},\quad r_{j+1}=x_{j}-\alpha_{j}Aq_{j}$

5: Compute the residual norm $\rho$, if $\rho$ small enough, then stop

6: $\beta_{i,j}=-(Ar_{j+1},Aq_{i})/(Aq_{i},Aq_{i}),\quad
q_{j+1}=r_{j+1}+\sum_{i=0}^{j}\beta_{i,j}q_{i}$

7:end for

In Algorithm 6, the term $Aq_{j}$ can be updated by means of the formula in
line 6 and the value of $Ar_{j}$ without an additional matrix-vector
multiplication. GCR generates a set of $A$-orthogonal basis vectors $q_{j}$
spanning the Krylov subspace (4). These vectors can be alternatively updated
by

$\beta_{i,j}=-(A^{2}q_{j},Aq_{i})/(Aq_{i},Aq_{i}),\quad
q_{j+1}=Aq_{j}+\sum_{i=0}^{j}\beta_{i,j}q_{i},$ (15)

leading to the so-called ORTHODIR algorithm, which was proposed by Jea and
Young [291, 147]. A truncated version of GCR was given earlier in [272], and
both GCR and ORTHODIR can be formulated in a cyclic manner like Algorithm 3 to
reduce the associated work and storage. Nonetheless, GCR can fail if the
symmetric part of $A$ is indefinite and the implementation of (15) can cause
stability problems, while it is known that MGS-GMRES is a cheaper and more
robust alternative; see, e.g., [223]. Similar reasoning holds for another
early variant proposed by Axelsson [10]. It was mentioned in [280] that there
is a close relationship between SGMRES and ORTHODIR; see also [150], in which
it was also shown that RB-SGMRES is closely related to GCR.

An early algorithm that realized (6) was already described by Khabaza [156] in
1963, which used $r_{0},Ar_{0},\dots,A^{n-1}r_{0}$ as the basis vectors and
thus suffered from numerical instabilities. In [82], Eirola and Nevanlinna
formalized a quasi-Newton method, later called EN method, that approximates
$A^{-1}$ by rank-one updates. It was shown in [82, 275] that a variant of EN
is mathematically equivalent to GMRES; see also [268] for additional comments.
More recent examples of equivalence relations involving quasi-Newton methods
can be found in [135, 136]. In 1988, on the other hand, Sidi gave a proof in
[233] that an acceleration technique called reduced rank extrapolation is
mathematically equivalent to GMRES; see [245] for a survey on vector
extrapolation methods. Later, it was established by Walker and Ni [279] that
another technique called Anderson acceleration is equivalent to GMRES in some
sense; we refer the reader to [279, 37] for more details on Anderson
acceleration and its relation to GMRES.

###  2.5 Convergence behavior

An important observation for studying GMRES convergence is that any vector in
the Krylov subspace (4) can be written in polynomial form, which implies

$r_{n}=p_{n}(A)r_{0},\quad p_{n}\in\mathcal{P}_{n},$

where $\mathcal{P}_{n}$ is the set of all polynomials $p$ of degree at most
$n$ such that $p(0)=1$. Then, (6) can be expressed as

$\left\lVert r_{n}\right\rVert=\left\lVert
p_{n}(A)r_{0}\right\rVert=\min_{p\in\mathcal{P}_{n}}\left\lVert
p(A)r_{0}\right\rVert.$ (16)

Let $\lambda_{i}(\cdot)$ denote the eigenvalues of a matrix and
$\lambda_{i}=\lambda_{i}(A)$. Assume that $A$ is diagonalizable so that
$A=X\Lambda X^{-1}$ with
$\Lambda=\operatorname{diag}(\lambda_{1},\dots,\lambda_{N})$. The residual
norm of GMRES satisfies

$\frac{\left\lVert r_{n}\right\rVert}{\left\lVert
r_{0}\right\rVert}\leq\kappa(X)\min_{p\in\mathcal{P}_{n}}\max_{i}\left\lvert
p(\lambda_{i})\right\rvert,$ (17)

where $\kappa(X)=\left\lVert X\right\rVert\left\lVert X^{-1}\right\rVert$
denotes the condition number of $X$. Note that $A$ being normal implies
$\kappa(X)=1$, in which case the bound in (17) is sharp (see [124, 152]). This
bound was presented in 1982 by Elman [88] in the context of the so-called GCR
algorithm; see also [83, 223]. It was also proved in [88] that if
$M=(A+A^{\mathsf{T}})/2$ is positive definite, then we have

$\frac{\left\lVert r_{n}\right\rVert}{\left\lVert
r_{0}\right\rVert}\leq\left(1-\frac{{\lambda_{\min}(M)}^{2}}{\lambda_{\max}(A^{\mathsf{T}}A)}\right)^{\frac{n}{2}},$

where $\lambda_{\min}(\cdot)$ and $\lambda_{\max}(\cdot)$ are respectively the
smallest and largest eigenvalues of a matrix in the absolute sense. This bound
was improved two decades later in [25, 24]. Now, it follows from (16) that

$\frac{\left\lVert r_{n}\right\rVert}{\left\lVert
r_{0}\right\rVert}\leq\psi_{n}(A)\leq\varphi_{n}(A),\quad\psi_{n}(A)=\max_{\left\lVert
v\right\rVert=1}\min_{p\in\mathcal{P}_{n}}\left\lVert
p(A)v\right\rVert,\quad\varphi_{n}(A)=\min_{p\in\mathcal{P}_{n}}\left\lVert
p(A)\right\rVert.$ (18)

In general, $\psi_{n}(A)$ is called worst-case GMRES value, as its purpose is
to study the worst-case behavior of the GMRES method. It is known that
$\psi_{n}(A)$ is attainable by the residual norm, that is, for each $n$ there
exists an initial vector $v$ such that $\left\lVert
r_{n}\right\rVert=\psi_{n}(A)$ (see, e.g., [258]). $\varphi_{n}(A)$ is called
ideal GMRES value, which was introduced in [129] to totally exclude the
influence of the initial vector. It was proved in [129, 170] that
$\varphi_{n}(A)$ has a unique minimizer, while as discussed in [100, 169], the
polynomial and unit-norm vector associated with $\psi_{n}(A)$ may not be
uniquely determined. Moreover, if $A$ is normal, then
$\psi_{n}(A)=\varphi_{n}(A)$ holds (see [124, 152]); otherwise, there exist
examples for which $\psi_{n}(A)<\varphi_{n}(A)$ (see [99, 260, 100]). We refer
the reader to [168, 169, 258, 100] and their references for further details on
(18).

The field of values plays a beneficial role in convergence analysis, which is
defined as

$\mathcal{F}(A)=\\{(Av,v):\left\lVert v\right\rVert=1,v\in\mathbb{C}^{N}\\}.$

Assume that the origin is outside $\mathcal{F}(A)$ and let

$\mu_{F}(A)=\min_{\lambda\in\mathcal{F}(A)}\left\lvert\lambda\right\rvert.$

Then, the following bound holds:

$\frac{\left\lVert r_{n}\right\rVert}{\left\lVert
r_{0}\right\rVert}\leq(1-\mu_{F}(A)\mu_{F}(A^{-1}))^{\frac{n}{2}};$ (19)

see [253, 80] for proofs based on the worst-case GMRES approximation; see also
[171] where the bound in (19) was proved to be hold for the ideal GMRES
approximation. Besides, an early account of the field of the values applied to
convergence analysis can be found in [79]. We refer the reader to the most
recent work [28] and the references therein for the use of the field of values
for convergence analysis. The polynomial numerical hull introduced in [200]
can be seen as a generalization of the field of values, which is defined as

$\mathcal{H}_{n}(A)=\\{\lambda\in\mathbb{C}:\left\lVert
p(A)\right\rVert\geq\left\lvert p(\lambda)\right\rvert,\deg(p)\leq n\\},$

from which a lower bound of the ideal GMRES approximation can be easily
derived. However, the determination of $\mathcal{H}_{n}(A)$ is very difficult
in its full generality. We refer to [122, 98, 123, 258] for related works.
Another group of results is obtained via the $\epsilon$-pseudospectrum [261].
Let $I$ denote the identity matrix of appropriate size. For a real
$\epsilon>0$, the $\epsilon$-pseudospectrum of $A$ is the set

$\Lambda_{\epsilon}=\\{\lambda\in\mathbb{C}:\left\lVert\lambda
I-A\right\rVert^{-1}\geq\epsilon^{-1}\\},$

from which an upper bound of $\varphi_{n}(A)$ can be obtained; see, e.g.,
[198, 262] for more details.

Linear contraction bounds for the residual norms can be somewhat misleading
since a highly nonlinear convergence behavior can not be adequately described
by a linear contraction. The book [167] written by Liesen and Strakoš is very
useful for getting more insight in the cost of computations (see also [54]),
which highlights the inherent nonlinear nature of Krylov subspace methods. In
[89], an experimental comparison of GMRES convergence bounds based on
eigenpairs, field of values and $\epsilon$-pseudospectrum for non-normal
matrices can be found, but none of these can be consistently better than the
others. In order to get more descriptive results, Titley-Peloquin et al. [259]
in 2014 derived several GMRES bounds that involve the initial residual vector;
see also [209] for an analysis with respect to nonstandard inner products.
More recently, Sacchi and Simoncini [226] gave a more descriptive convergence
analysis specifically for the case of localized ill-conditioning.

It was confirmed in a series of papers by Arioli, Greenbaum, Pták and Strakoš
[128, 126, 8] from 1994 to 1998 that convergence of GMRES for non-normal
matrices could not be determined solely by the distribution of eigenvalues. In
[128], Greenbaum and Strakoš studied the matrices $B$ for which the spaces
$A\mathcal{K}_{n}(A,r_{0})$ and $B\mathcal{K}_{n}(B,r_{0})$ are the same. As a
result, the convergence history generated for the $(B,r_{0})$ pair by GMRES is
the same as that generated for the $(A,r_{0})$ pair. The matrices $B$ that
satisfies this property are called GMRES-equivalent matrices (see, e.g., [163,
78]). Among other results, Greenbaum and Strakoš concisely state in [128] on
GMRES convergence that

> Any behavior that can be seen with the method can be seen with the method
> applied to a matrix having any nonzero eigenvalues.

This has been extended by Greenbaum et al. [126] along a different line, who
conclude that

> Any nonincreasing convergence curve can be obtained with GMRES applied to a
> matrix having any desired eigenvalues.

Finally, a complete parametrization of all $(A,b)$ pairs with prescribed
eigenvalues for which GMRES generates prescribed residual norms has been given
in [8]; see also [183] for further discussion. GMRES-equivalent matrices have
been employed by some subsequent studies [163, 78], while investigations of
“any behavior is possible”, like that in [128, 126, 8], have been continued in
the past ten years, concerning Ritz values [75, 76], harmonic Ritz values
[74], restarted GMRES [271, 77], and block GMRES [158]. Note that Ritz and
harmonic Ritz values are respectively the roots of FOM and GMRES residual
polynomials; see [119, 185]. We also mention [183] for a generalization to
other Krylov subspace methods.

Eigenvalue clusters associated with a parameterized model were employed in
[47] to explain GMRES convergence. Exact but complicated expressions for
residual norms were derived in [144, 165, 228]; see also [186]. Inspired by
the work [265] on the conjugate gradient method, van der Vorst and Vuik [267]
in 1993 investigated the superlinear convergence behavior of GMRES by means of
Ritz values. Simoncini and Szyld [242] addressed the same problem using
spectral projectors. Moret [189] discussed superlinear convergence of GMRES in
a complex separable Hilbert space; see also [33] for a generalization of
Moret’s result. Since any nonincreasing convergence curve is possible for
GMRES [126], on the other hand, stagnation may occur even for the original
GMRES method without restarting. This was first investigated by Brown [38] in
1991, and was then discussed in [293, 169, 244, 237, 182, 184]. Although most
work on convergence analysis considered full GMRES, the restarted version is
more practical, on which we refer the reader to, e.g., [151, 294, 295, 16,
270, 271, 77] for more details.

The above discussion assumes exact arithmetic. In finite precision arithmetic,
it was proved by Drkošová et al. [73] in 1995 that HH-GMRES is backward
stable. In 2006, it was proved by Paige et al. [204] that MGS-GMRES is also
backward stable; previous observations can be found in [127, 206]. Due to the
ill-conditioning of the matrix $Z_{n}$ in (12), SGMRES has serious stability
problems (see, e.g, [164, 56]), which can be partially remedied by using RB-
SGMRES [150] or the adaptive variant (14) [149]; see also the experiments in
[180].

##  3 Acceleration strategies

The focus in this section is on the acceleration strategies to tackle the
challenge that a basic GMRES algorithm might face. First, it is almost certain
that the use of preconditioners will be beneficial for Krylov subspace methods
when solving ill-conditioned problems. There are also numerous publications
concerning deflation and augmentation, some of which are highly related to
preconditioning. Then we briefly introduce weighted inner products and inexact
matrix-vector products, and end up this section by discussing mixed-precision
techniques.

###  3.1 Preconditioning

Preconditioning of the linear system is used to improve the convergence
behavior of Krylov subspace methods. One needs to find a nonsingular matrix
$M$ close to $A$ in some sense, commonly called preconditioning matrix or
preconditioner, such that the new system involving $M$ should be inexpensive
to solve, that is,

$\displaystyle M^{-1}Ax=M^{-1}b,$ (20) $\displaystyle AM^{-1}u=b,\quad Mx=u.$
(21)

The cases (20) and (21) correspond to left preconditioning and right
preconditioning, respectively. Another form called split preconditioning can
be employed when nearly symmetric systems are encountered; see, e.g., [220].
For GMRES, right preconditioning is often preferred because the residuals in
this case are identical to the true residuals, only the solution update
formula being transformed to $x_{n}=x_{0}+M^{-1}V_{n}y$, while in another case
left preconditioning leads to preconditioned residuals
$z_{n}=M^{-1}(b-Ax_{n})$. In addition, from an algorithmic point of view, one
only needs to transform line 3 in Algorithm 2 to $w_{j}=M^{-1}Av_{j}$ for the
case (20) and $w_{j}=AM^{-1}v_{j}$ for (21). High-quality preconditioners are
often constructed by incorporating problem-specific information; see, e.g.,
[201, 95, 96, 109, 208]. On the other hand, general-purpose preconditioners
have also been broadly exploited, such as incomplete factorization,
approximate inverse, and algebraic multigrid; see [27, 281] and references
therein.

The flexible MGS-GMRES variant proposed by Saad [217] in 1993 allows the
preconditioner to vary at each step, which shall be denoted by FGMRES as
widely used in the literature.

Algorithm 7 FGMRES

1:$\beta=\left\lVert r_{0}\right\rVert,\quad v_{1}=r_{0}/\beta$

2:for $j=1,2,\dots,n$ do

3: $z_{j}=M_{j}^{-1}v_{j},\quad w_{j}=Az_{j}$

4: for $i=1,2,\dots,j$ do

5: $h_{i,j}=(w_{j},v_{i}),\quad w_{j}=w_{j}-h_{i,j}v_{i}$

6: end for

7: $h_{j+1,j}=\left\lVert w_{j}\right\rVert,\quad v_{j+1}=w_{j}/h_{j+1,j}$

8:end for

9:Compute $y$ such that $\left\lVert\beta e_{1}-\mkern
3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern 1.0mu_{n}y\right\rVert$ is
minimized

10:$x_{n}=x_{0}+Z_{n}y$

As seen in Algorithm 7, $z_{j}=M_{j}^{-1}v_{j}$ should be explicitly stored in
$Z_{n}=[z_{1},\dots,z_{n}]$, and then reused in line 10. Therefore, the
Arnoldi-like relation $AZ_{n}=V_{n+1}\mkern
3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern 1.0mu_{n}$ holds. Unlike
standard right preconditioning, which constructs an orthonormal basis of the
Krylov subspace

$\operatorname{span}\\{r_{0},AM^{-1}r_{0},\dots,(AM^{-1})^{n-1}r_{0}\\},$

the subspace of solution estimates in FGMRES is no longer a standard Krylov
subspace. In addition, $H_{n}$ may be singular even if $A$ is nonsingular, and
thus the assumption of nonsingularity of $H_{n}$ must be made when deducing
exact solution from $h_{j+1,j}=0$; see [217, 220] for more discussion. A very
detailed analysis of FGMRES, in comparison with right preconditioned MGS-
GMRES, can be found in [7]. More recently, Greif et al. [130] developed the
multi-preconditioned GMRES scheme in 2017, which shares some similarities with
the block version of FGMRES [43, 42], but targets (1) instead of systems with
multiple right-hand sides. In practice, a restarted procedure should be used
in replacement of the simple version in Algorithm 7; see [103] for details on
the implementation of FGMRES.

Another interesting but less used variable preconditioning technique proposed
by van der Vorst and Vuik [268] in 1994 consists of the GCR algorithm as the
outer iterations and a GMRES algorithm as the inner iterations. The inner
algorithm constructs a preconditioner for the outer procedure. This new scheme
is called GMRESR, which was motivated by the EN algorithm [82]. Note that in
GMRESR the residuals are preconditioned, while in FGMRES, as shown in
Algorithm 7, the search directions are preconditioned. We refer to [274] for a
comparison of FGMRES and GMRESR. Extensions of GMRESR have been developed in
[64, 65] by requiring that the inner basis maintains orthogonality to the
outer basis, possibly combining the truncation strategy to reduce storage
costs.

The GMRES polynomial itself, defined as the polynomial in (16), can be used to
derive preconditioners for Krylov subspace algorithms. Liu et al. [174] in
2015 developed a polynomial preconditioned GMRES variant using the GMRES
polynomial. An improved approach was proposed in [176] by Loe and Morgan. The
idea can be summarized as follows:

  1. 1.

Run $m$ steps of MGS-GMRES;

  2. 2.

Compute the harmonic Ritz values (see (24)) which are the roots of the GMRES
polynomial;

  3. 3.

Order the roots and construct the polynomial;

  4. 4.

Apply a GMRES algorithm to the polynomial preconditioned system.

We refer the reader to [174, 178, 177, 176] for more details and experiments.
Polynomial preconditioners can be combined with other preconditioners and were
highlighted as effective tools for parallel computing. Early contributions of
hybrid or polynomial preconditioned GMRES algorithms can be found in [199,
254, 152, 269]. The stability problems arising in high degree preconditioners
were addressed in [91, 176, 290].

In order to overcome scaling difficulties for large number of processor cores,
two general iterative schemes were introduced in 2014 by McInnes et al. [181].
As an example, they presented a two-level hierarchical FGMRES algorithm. The
outer level applies FGMRES to the global problem, while the inner level
executes inexact GMRES within subgroups of cores. Therefore, the inner solver
serves as a variable preconditioner for the outer solver and entails less
communication overhead. The two-level scheme can be extended to multiple
levels with multiple solvers, leading to the hierarchical Krylov scheme.
Meanwhile, a nested version was also presented in [181]. The inner solvers in
the nested scheme apply inexact Krylov algorithms to the global system instead
of local systems.

###  3.2 Deflation and augmentation

The convergence of GMRES($m$) is often hampered by eigenvalues of small
magnitude and loss of information at the end of each cycle. Deflation
techniques can eliminate certain small eigenvalues or move certain small
eigenvalues away from the origin. Note that the word “deflation” has multiple
meanings in the literature which may lead to confusion. Augmentation amounts
to considering an $m$-dimensional subspace
$\mathcal{K}_{m_{1}}+\operatorname{span}\\{u_{1},\dots,u_{m_{2}}\\}$ with
$m=m_{1}+m_{2}$, adding information from previous cycles to the search space.

Let us write

$Z_{m}=[z_{1},\dots,z_{m}]=[v_{1},\dots,v_{m_{1}},u_{1},\dots,u_{m_{2}}],$
(22)

where the first $m_{1}$ vectors are the orthonormal Arnoldi vectors. In 1995,
Morgan [191] presented an augmented GMRES algorithm, later called GMRES-E in
[192], where $u_{1},\dots,u_{m_{2}}$ are chosen as the harmonic Ritz vectors
corresponding to the harmonic Ritz values of smallest magnitude; see also [55,
218] for more analysis of this approach. Given the relation

$AZ_{m}=V_{m+1}\mkern 3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern
1.0mu_{m},$ (23)

where $V_{m+1}$ forms a set of orthonormal vectors, it is known that solving
the harmonic eigenvalue problem

$Z_{m}^{\mathsf{T}}A^{\mathsf{T}}Z_{m}y_{i}=\frac{1}{\theta_{i}}Z_{m}^{\mathsf{T}}A^{\mathsf{T}}AZ_{m}y_{i}$

is equivalent to solving

$(H_{m}+h_{m+1,m}^{2}H_{m}^{-\mathsf{T}}e_{m}e_{m}^{\intercal})y_{i}=\theta_{i}y_{i},$
(24)

with $u_{i}=Z_{m}y_{i}$. Morgan [192] in 2000 proposed a mathematically
equivalent but more efficient algorithm, commonly called GMRES-IR, which uses
the implicit restarting technique with unwanted harmonic Ritz values as shifts
to generate augmented subspaces (see [160]). A particularly important
observation is that the augmented subspace can be formulated as a Krylov
subspace with a carefully chosen starting vector $\hat{r}_{0}$, that is,

$\operatorname{span}\\{\hat{r}_{0},A\hat{r}_{0},\dots,A^{m-1}\hat{r}_{0}\\}=\operatorname{span}\\{r_{0},Ar_{0},\dots,A^{m_{1}-1}r_{0},u_{1},\dots,u_{m_{2}}\\}.$

We mention that Le Calvez and Molina [159] developed independently an
implicitly restarted GMRES algorithm, which, however, resorts to a different
eigenvalue problem instead of (24). Further, Morgan [193] revisited GMRES-IR
in 2002 and found that there was room for improvement. The so-called GMRES
with deflated restarting (GMRES-DR) algorithm [193], based on the thick
restarting approach (see [283]), can be simpler and more stable than GMRES-IR;
see also [214] for a slightly modified variant. Besides, there are many
situations where variable preconditioners are worth to be considered, and thus
a flexible variant of GMRES-DR (FGMRES-DR) was suggested in [115]. We point
here to the more recent work of Liu et al. [174] who applied polynomial
preconditioning to GMRES-DR. In addition, the effectiveness of GMRES-DR was
studied by Morgan et al. [195]. The main idea is that the approximate
eigenvectors generated by GMRES-DR can be seen as pseudoeigenvectors. It was
shown in [195] by examples that GMRES-DR can also work for highly nonnormal
matrices.

Some deflation strategies are highly related to preconditioning. Kharchenko
and Yeremin [157] in 1995 suggested the use of low-rank transformations by
means of left and right Ritz vectors to construct a right preconditioner. As a
result, extremal eigenvalues would be translated into a vicinity of one. In
1996, Erhel et al. [94] employed a different approach to estimate the
invariant subspace and constructed a right preconditioner that can move the
small eigenvalues to a multiple large eigenvalue. A similar idea was later
suggested by Baglama et al. [12]. They described two algorithms in the context
of left preconditioning, which exploit implicit restarting to improve
flexibility. Moreover, an improved algorithm built on the previous work [191,
94, 12] was introduced in [41]; see also [197] for an adaptive strategy for
determining the restart frequency.

In 1996, de Sturler [64] proposed an inner-outer algorithm with augmented
basis, called GCRO, which was extended to the so-called GCROT algorithm in
[65] by incorporating truncation strategies. A beautiful account of Morgan’s
work, deflation by preconditioning, and GCROT can be found in the survey [81].
It was observed by Baker et al. [17] in 2005 that GMRES($m$) has alternating
behavior, namely, the residual vector often alternates direction at the end of
each cycle, leading to slow convergence and even stalling. They gave some
remedies to overcome this problem by combining the ideas of augmented
subspaces [191] and reserved error approximations [64], resulting in the
LGMRES($m_{1},m_{2}$) algorithm.

Algorithm 8 LGMRES($m_{1},m_{2}$)

1:$m=m_{1}+m_{2},\quad\beta=\left\lVert r_{0}\right\rVert,\quad
v_{1}=r_{0}/\beta$

2:for $k=0,1,\dots$ do

3: for $j=1,2,\dots,m$ do

4: if $j\leq m_{1}$ or $j-m_{1}>k$ then

5: $w_{j}=Av_{j}$

6: else

7: $w_{j}=Au_{k-(j-m_{1}-1)}$

8: end if

9: for $i=1,2,\dots,j$ do

10: $h_{i,j}=(w_{j},v_{i}),\quad w_{j}=w_{j}-h_{i,j}v_{i}$

11: end for

12: $h_{j+1,j}=\left\lVert w_{j}\right\rVert,\quad v_{j+1}=w_{j}/h_{j+1,j}$

13: end for

14: Compute $y$ such that $\left\lVert\beta e_{1}-\mkern
3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern 1.0mu_{m}y\right\rVert$ is
minimized

15: if $k-m_{2}<0$ then

16: $Z_{m}=[v_{1},\dots,v_{m_{1}+m_{2}-k},u_{k},\dots,u_{1}]$

17: else

18: $Z_{m}=[v_{1},\dots,v_{m_{1}},u_{k},\dots,u_{k-m_{2}+1}]$

19: end if

20: $u_{k+1}=Z_{m}y,\quad x_{k+1}=x_{k}+u_{k+1},\quad r_{k+1}=b-Ax_{k+1}$

21: $\beta=\left\lVert r_{k+1}\right\rVert$, if $\beta$ small enough, then
stop; else $v_{1}=r_{k+1}/\beta$

22:end for

Algorithm 8 satisfies (23) and generalizes GMRES-E to a GCRO-like vector
selection strategy. Notice that $u_{k+1}$ can be considered as an
approximation to the exact error, which is the key idea involved in both GCRO
and LGMRES; see [64, 17]. We refer the reader to [17] for further details and
numerical experiments; see also [142] for a similar approach. On the other
hand, combining the main features of GMRES-DR and GCRO, in 2006, an algorithm
called GCRO-DR was derived in [207] for sequences of linear systems. When
solving a single linear system, GCRO-DR is mathematically equivalent to GMRES-
DR.

It was Gutknecht [132] who attempted to present deflated and augmented Krylov
subspace methods in a common framework. Also, Gutknecht investigated the
breakdown conditions of deflated GMRES; we also refer the interested reader to
the companion paper [110]. These two papers might provide good starting point
for further investigations of deflation and augmentation.

###  3.3 Weighted inner products

In 1998, it was shown by Essai [97] that a suitably chosen inner product could
improve the performance of restarted MGS-GMRES. Let
$D=\operatorname{diag}(\delta_{1},\dots,\delta_{N})$ with $\delta_{i}>0$. The
$D$-inner product is defined as $(u,v)_{D}=v^{\intercal}Du$ for all
$u,v\in\mathbb{R}^{N}$, and then the associated $D$-norm is $\left\lVert
u\right\rVert_{D}=\sqrt{u^{\intercal}Du}$. From an algorithmic viewpoint, for
example in Algorithm 1, one only needs to replace lines 5 and 7 with the new
inner product and norm. The resulting algorithm shall be denoted by WGMRES. In
[97], the diagonal matrix $D$ was chosen as $\delta_{i}=\sqrt{N}\left\lvert
r_{0,i}\right\rvert/\left\lVert r_{0}\right\rVert$ and updated at each
restart.

Preconditioned variants were studied in [48]. Later, a modified version
motivated by the idea of [11] was proposed in [225], in which $\delta_{i}$ are
chosen randomly and there is no need to use Givens rotations. In [202], WGMRES
was combined with the augmentation technique, and then the new variant was
compared with LGMRES as illustrated in Algorithm 8. On the other hand, Embree
et al. [92] in 2017 applied weighting to GMRES-DR and observed that the mixed
scheme could be better than using deflation alone or weighting alone. In 2014,
Güttel and Pestana [134] proposed a new variant that applied the
orthogonalization process to the transformed matrix
$D^{\frac{1}{2}}AD^{-\frac{1}{2}}$ and starting vector $D^{\frac{1}{2}}r_{0}$
with respect to the Euclidean inner product, which is mathematically
equivalent to WGMRES. For more insight in nonstandard inner products, see
[209, 134, 92] and references therein.

###  3.4 Inexact matrix-vector products

An inexact process might be beneficial if the computation of the matrix-vector
product $Av_{j}$ is time-consuming. From a mathematical point of view, an
inexact matrix-vector product can be modeled by $(A+E_{j})v_{j}$ where $E_{j}$
is a perturbation matrix at $j$th iteration. From (10), one obtains that

$AV_{n}+[E_{1}v_{1},\dots,E_{n}v_{n}]=V_{n+1}\mkern
3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern 1.0mu_{n}.$

The inexact GMRES algorithm was first proposed in a technical report by Bouras
and Frayssé [34] in 2000 (later published in [35] in 2005), and then analyzed
by Simoncini and Szyld [241, 242], van den Eshof and Sleijpen [264], Giraud et
al. [114], and Sifuentes et al. [236]. A relaxation strategy on the inner
accuracy of GMRES can be found in [35]. It was shown in [241] that the early
matrix-vector products must be computed with high accuracy, but can be relaxed
as the iteration progresses, which explains the term “relaxation”. An
interesting review of inexact GMRES was given in [243]; see [235] for
performance evaluation. In 2013, Dolgov [72] employed relaxation techniques to
improve the performance of tensor GMRES; see [20, 72] for more details on
tensor formats.

###  3.5 Mixed precision

On modern computer architectures, single precision arithmetic (32 bits) is
cheaper than double precision arithmetic (64 bits) in the sense that
communication and computation costs grow with the size of the floating point
format [1]. In the 2008 revision of the IEEE standard, half precision (16
bits) was defined as a storage format. Mixed-precision numerical algorithms
can date back to more than five decades ago [188]. Recent advances in
hardware-level support (e.g., NVIDIA V100 GPU) for low-precision arithmetic
and the increasing impact of communication costs have promoted the reuse of
mixed-precision techniques; see [1] for a recent survey.

It is very natural to extend Algorithm 2 to the case of mixed precision. Here,
we only consider the two-precision algorithm, which exploits single precision
everywhere, except that

  1. 1.

The original system $Ax=b$ must be stored in double precision;

  2. 2.

Although $V_{n}y$ could be computed in single precision, the update of $x_{n}$
should be done in double precision;

  3. 3.

For the restarted version, the computation of residual vector should also be
done in double precision.

A recent paper by Gratton et al. [120] showed that using low precision in MGS-
GMRES could achieve the same convergence rate and final accuracy as the full
precision variant. They provided a theoretical analysis for the non-restarted
algorithm by means of inexact inner products; see also [4]. Then, Lindquist et
al. [172] focused on the experimental aspects of restarted mixed-precision
GMRES algorithms, including MGS-GMRES and a variant based on CGS with one
reorthogonalization (CGS2); see [117, 116] for error analysis of CGS2.

Algorithm 9 CGS2-Arnoldi

1:$v_{1}=r_{0}/\left\lVert r_{0}\right\rVert$

2:for $j=1,2,\dots,n$ do

3: $w_{j}=Av_{j}$

4: $\hat{h}=V_{j}^{\mathsf{T}}w_{j},\quad w_{j}=w_{j}-V_{j}\hat{h}$

5: $\tilde{h}=V_{j}^{\mathsf{T}}w_{j},\quad w_{j}=w_{j}-V_{j}\tilde{h}$

6: $[h_{1,j},\dots,h_{j,j}]^{\mathsf{T}}=\hat{h}+\tilde{h}$

7: $h_{j+1,j}=\left\lVert w_{j}\right\rVert$, $v_{j+1}=w_{j}/h_{j+1,j}$

8:end for

In Algorithm 9, we illustrate the CGS2-Arnoldi procedure. Motivated by the
idea in [120], a theoretical result for CGS2-GMRES was provided in [173],
along with more experimental results.

Another interesting class of mixed-precision algorithms involves iterative
refinement where low-precision GMRES is employed as inner solver. In 1992,
Turner and Walker [263] showed that GMRES($m$) can be regarded as an iterative
refinement process. In 2009, Arioli and Duff [6] proved that the FGMRES
algorithm preconditioned by low-precision factorization is backward stable.

Algorithm 10 GMRES-IR

1:Compute the LU factorization $A=LU$ in single precision

2:Solve $Mx_{0}=b$ in single precision with $M=LU$

3:$r_{0}=b-Ax_{0}$ in double precision

4:for $k=0,1,\dots$ do

5: Solve $M^{-1}Ad_{k}=M^{-1}r_{k}$ by a GMRES algorithm in single precision
with $M=LU$

6: $x_{k+1}=x_{k}+d_{k}$ in double precision

7: $r_{k+1}=b-Ax_{k+1}$ in double precision

8: Check convergence

9:end for

We sketch GMRES-based iterative refinement (GMRES-IR) in Algorithm 10; see,
e.g., [51, 175] for more details. Note that this is just a simple example
illustrated using two precisions, for which three or more precisions could
also be exploited by relaxing some of the operations. Carson and Higham [50]
in 2017 gave a new forward analysis of iterative refinement and suggested the
use of GMRES-IR. Then, in [51], they introduced and analyzed an iterative
refinement scheme that can possibly exploit three different precisions. It was
shown that GMRES-IR could provide the required relative accuracy, even though
the LU factorization in single precision might be of low quality. The error
analysis in [50, 51] was recently generalized in [5] to the case where the
GMRES solver and matrix-vector multiplications in line 5 may proceed in
independent precisions, yielding a possibility of using five precisions. On
the other hand, the overflow and underflow problems associated with mixed-
precision algorithms were investigated by Higham et al. [139].

The survey article [1] provides additional details and references on mixed-
precision techniques. Moreover, although not directly connected to GMRES, the
studies in [288] and [289] can deepen our insight into mixed-precision
orthogonalization algorithms.

##  4 Parallel algorithms

On modern computer architectures, communication is expensive relative to
computation. Sparse matrix-vector products (SpMVs) and dot products,
especially the latter, often limit the possible speedups for parallel
algorithms. Therefore, the high cost of global reduction operations in the
orthogonalization process makes GMRES hard to parallelize. To alleviate
performance bottleneck, much work has focused on ways to reduce or hide
communication costs. Early work concerning the parallelization of GMRES can be
found in, e.g., [277, 71, 232, 60]; see also [68] for a general discussion on
parallel algorithms. Here we focused on the $s$-step GMRES algorithm, and then
present the idea of pipelining. These approaches have gained popularity in the
last decade. Finally, some other potential strategies are included at the end
of this section.

###  4.1 Matrix powers and polynomial basis

In [69] and some subsequent publications [70, 187, 141, 21], a Krylov subspace
algorithm is thought to contain several “kernels”, each kernel being an
important and time-consuming building block of that algorithm. In GMRES,
matrix powers computation is thus considered as a kernel, which computes $s$
vectors $Av,\dots,A^{s}v$. The most trivial implementation requires $s$
messages between each processor and its neighbors, whereas a well-designed
matrix powers kernel (MPK) may require only $1+o(1)$ message at the cost of
redundant computation and storage consumption; see [69], in which data
dependencies and performance models are discussed in terms of both distributed
and shared memory systems. Nonetheless, this technique has several
disadvantages. First, MPK introduces additional overheads. Second, it’s still
difficult to apply effective preconditioners with MPK [141, 285]. Third, basis
vectors will become more and more linearly dependent [63].

The last drawback can be partially corrected by using polynomial basis, that
is,

$W_{s}=[w_{1},\dots,w_{s}]=[\phi_{0}(A)w,\phi_{1}(A)w,\dots,\phi_{s-1}(A)w],$

where $\phi_{j}(z)$ is a polynomial of degree $j$. Choosing

$\phi_{j}(z)=z^{j}$ (25)

yields the monomial basis, which converges to the eigenvector corresponding to
the dominant eigenvalue of $A$, and thus, as mentioned above, will rapidly
become numerically rank deficient. Two alternatives have been suggested. One
is to use Newton polynomials

$\phi_{0}(z)=1,\quad\phi_{j}(z)=\varrho_{j}(z-\theta_{j})\phi_{j-1}(z),$ (26)

as suggested by Bai et al. [14] in 1994, in which an $s$-step GMRES algorithm
was proposed. Here, $\varrho_{j}$ are scaling factors, which are generally
taken as $1$, and $\theta_{j}$ are eigenvalue estimates computed by a few
iterations of the Arnoldi process. The other is Chebyshev polynomials, which
can be defined by

$\phi_{0}(z)=1,\quad\phi_{1}(z)=\frac{1}{2\gamma}(z-\zeta)\phi_{0}(z),\quad\phi_{j}(z)=\frac{1}{\gamma}((z-\zeta)\phi_{j-1}(z)-\frac{\tau^{2}}{4\gamma}\phi_{j-2}(z)),$
(27)

where the spectrum is assumed to be circumscribed by the rectangle
$\\{z=z_{1}+\iota z_{2}:\left\lvert
z_{1}-\zeta\right\rvert\leq\xi_{1},\left\lvert
z_{2}\right\rvert\leq\xi_{2}\\}$ with $\iota$ the imaginary unit. Then, we
choose $\gamma=\max\\{\xi_{1},\xi_{2}\\}$ and
$\tau=\sqrt{\xi_{1}^{2}-\xi_{2}^{2}}$, so that $\zeta\pm\tau$ are the foci of
the ellipse for the spectrum. Chebyshev polynomials for GMRES were discussed
by Joubert and Carey [153, 154] in 1992; more details of Krylov subspace bases
can be found in [141, 210, 49].

###  4.2 $s$-step GMRES

The idea of performing $s$ iterations at once can date back to the early
1950s, as quoted in [102]. In 1990, Chronopoulos and Kim [58] developed an
$s$-step GMRES algorithm using monomial basis. A Newton basis $s$-step GMRES
algorithm was considered by Bai et al. [14], in which modified Leja ordering
and basis scaling were employed to enhance robustness. They also suggested the
use of parallel QR factorization, and this was implemented by Erhel [93] in
1995 by means of the so-called RODDEC algorithm. On the other hand, Joubert
and Carey [153, 154] in 1992 considered an $s$-step GMRES algorithm using
Chebyshev basis along with an MPK for 2D regular meshes; see also [66]. In
2009, Mohiyuddin et al. [187] introduced an improved variant, for which a new
kernel called tall skinny QR was derived to achieve the optimal reduction of
communication while maintaining the accuracy; see [67]. Then, the thesis by
Hoemmen [141] gave a complete treatise of $s$-step GMRES.

Let $\mkern 3.0mu\overline{\mkern-3.0muB\mkern-1.0mu}\mkern 1.0mu_{s}$ be an
$(s+1)\times s$ matrix such that

$AW_{s}=W_{s+1}\mkern 3.0mu\overline{\mkern-3.0muB\mkern-1.0mu}\mkern
1.0mu_{s}.$ (28)

It is called a basis conversion matrix, by which the three bases with respect
to (25)–(27) can be exploited in an elegant manner. It was shown in [141] that
the QR factorization of $W_{s}$ could produce the same orthogonal matrix as
that in (10). Then, it follows that

$W_{s}=V_{s}T_{s},$ (29)

where $T_{s}\in\mathbb{R}^{s\times s}$ is an upper triangular matrix. Along
with (10a) and (28), this implies

$\mkern 3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern
1.0mu_{s}=T_{s+1}\mkern 3.0mu\overline{\mkern-3.0muB\mkern-1.0mu}\mkern
1.0mu_{s}T_{s}^{-1}.$

The above formulation can potentially bring two benefits. First, one could
compute $W_{s+1}$ with less communication than in the classical case. Second,
the QR factorization $W_{s+1}=V_{s+1}T_{s+1}$ could be performed in parallel.

The tall skinny QR (TSQR) algorithm studied in [141, 67] is a good choice for
reducing communication when there are many more rows than columns. Let $W_{s}$
be vertically partitioned into $\nu_{0}$ submatrices such that the number of
rows of each submatrix is larger than $s$. Then, we perform the QR
factorization of each block entry, that is,

$W_{s}=\left(\begin{array}[]{c}W_{s,1}\\\ W_{s,2}\\\ \vdots\\\
W_{s,\nu_{0}}\end{array}\right)=\left(\begin{array}[]{c}Q_{1}^{(0)}R_{1}^{(0)}\\\
Q_{2}^{(0)}R_{2}^{(0)}\\\ \vdots\\\
Q_{\nu_{0}}^{(0)}R_{\nu_{0}}^{(0)}\end{array}\right)=\left(\begin{array}[]{cccc}Q_{1}^{(0)}&&&\\\
&Q_{2}^{(0)}&&\\\ &&\ddots&\\\
&&&Q_{\nu_{0}}^{(0)}\end{array}\right)\left(\begin{array}[]{c}R_{1}^{(0)}\\\
R_{2}^{(0)}\\\ \vdots\\\ R_{\nu_{0}}^{(0)}\end{array}\right).$

Once such initialization is complete, in each cycle we perform the QR
factorization of vertically stacked upper triangular submatrices

$\left(\begin{array}[]{c}R_{1}^{(0)}\\\
R_{2}^{(0)}\end{array}\right)=Q_{1}^{(1)}R_{1}^{(1)},\quad\left(\begin{array}[]{c}R_{3}^{(0)}\\\
R_{4}^{(0)}\end{array}\right)=Q_{2}^{(1)}R_{2}^{(1)},\quad\left(\begin{array}[]{c}R_{1}^{(1)}\\\
R_{2}^{(1)}\end{array}\right)=Q_{1}^{(2)}R_{1}^{(2)},\quad\dots,$

with $R_{i}^{(j)}\in\mathbb{R}^{s\times s}$, and define

$Q^{(j)}=\left(\begin{array}[]{cccc}Q_{1}^{(j)}&&&\\\ &Q_{2}^{(j)}&&\\\
&&\ddots&\\\ &&&Q_{\nu_{j}}^{(j)}\end{array}\right),\quad
R^{(j)}=\left(\begin{array}[]{c}R_{1}^{(j)}\\\ R_{2}^{(j)}\\\ \vdots\\\
R_{\nu_{j}}^{(j)}\end{array}\right).$

Here, $\nu_{j}=\lceil\nu_{j-1}/2\rceil$ for
$j=1,2,\dots,\lceil\log_{2}{\nu_{0}}\rceil$. If $\nu_{j-1}\bmod 2\neq 0$, then
$Q_{\nu_{j}}^{(j)}=I$ and $R_{\nu_{j}}^{(j)}=R_{\nu_{j-1}}^{(j-1)}$. Algorithm
11 shows this procedure.

Algorithm 11 TSQR

1:for $i=1,2,\dots,\nu_{0}$ do

2: Compute QR factorization $W_{s,i}=Q_{i}^{(0)}R_{i}^{(0)}$

3:end for

4:for $j=1,2,\dots,\lceil\log_{2}{\nu_{0}}\rceil$ do

5:
$\hat{\nu}_{j}=\lfloor\nu_{j-1}/2\rfloor,\quad\nu_{j}=\lceil\nu_{j-1}/2\rceil$

6: for $i=1,2,\dots,\hat{\nu}_{j}$ do

7: Compute QR factorization $\left(\begin{array}[]{c}R_{2i-1}^{(j-1)}\\\
R_{2i}^{(j-1)}\end{array}\right)=Q_{i}^{(j)}R_{i}^{(j)}$

8: end for

9: if $\nu_{j}\neq\hat{\nu}_{j}$ then

10: $Q_{\nu_{j}}^{(j)}=I,\quad R_{\nu_{j}}^{(j)}=R_{\nu_{j}}^{(j-1)}$

11: end if

12:end for

Notice that
$Q^{(\lceil\log_{2}{\nu_{0}}\rceil)}=Q_{1}^{(\lceil\log_{2}{\nu_{0}}\rceil)}$
and
$R^{(\lceil\log_{2}{\nu_{0}}\rceil)}=R_{1}^{(\lceil\log_{2}{\nu_{0}}\rceil)}$.
Then, along with (29), it follows that

$V_{s}=Q^{(0)}Q^{(1)}\dots Q^{(\lceil\log_{2}{\nu_{0}}\rceil)},\quad
T_{s}=R^{(\lceil\log_{2}{\nu_{0}}\rceil)},$

yielding the desired factorization of $W_{s}$.

A large $s$ may lead to stability problems even when a polynomial basis is
used. In [187, 141] it was thus suggested to combine TSQR with the block Gram-
Schmidt (BGS) algorithm, such that $s$ can be chosen independently of the
restart length $m$. As a result, one can pursue both numerical stability and
convergence rate. Now, assume that $m=s\cdot t$. The BGS process can be
written as

$\mathfrak{R}_{s}=\mathfrak{V}_{s(j-1)}^{\mathsf{T}}W_{s},\quad
W_{s}=W_{s}-\mathfrak{V}_{s(j-1)}\mathfrak{R}_{s},$

where $\mathfrak{V}_{s(j-1)}$ denotes the orthogonal matrix that contains all
previous orthogonal basis vectors established in this cycle. After performing
the TSQR procedure, we can get
$\mathfrak{V}_{sj}=[\mathfrak{V}_{s(j-1)},V_{s}]$. In addition, let
$\mathfrak{H}_{sj}$ be the upper Hessenberg matrix corresponding to
$\mathfrak{V}_{sj}$. Then the BGS-based $s$-step GMRES satisfies the relation

$A\mathfrak{V}_{sj}=\mathfrak{V}_{sj+1}\mkern
3.0mu\overline{\mkern-3.0mu\mathfrak{H}\mkern-1.0mu}\mkern 1.0mu_{sj},$

which is equivalent to (10a) if choosing $n=sj$. Here, BGS amounts to a block
version of CGS. For more details on BGS, we refer the reader to [52, 53] and
references therein. BGS and TSQR are both regarded as kernels (see [187,
141]).

Notice that the last basis vector of $\mathfrak{V}_{sj+1}$ can be used as the
starting vector for the next $W_{s}$. Hence, we use the acute accent to
indicate a shift of indices such as
$\acute{W}_{s}=[w_{2},w_{3},\dots,w_{s+1}]$, on which BGS and TSQR would be
performed.

Algorithm 12 $s$-step GMRES

1:$\beta=\left\lVert r_{0}\right\rVert,\quad w_{1}=r_{0}/\beta$

2:for $j=1,2,\dots,t$ do

3: Select $\mkern 3.0mu\overline{\mkern-3.0muB\mkern-1.0mu}\mkern 1.0mu_{s}$
and compute $\acute{W}_{s}=[w_{2},w_{3},\dots,w_{s+1}]$

4: if $j=1$ then

5: Compute QR factorization $W_{s+1}=V_{s+1}T_{s+1}$

6: $\mkern 3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern
1.0mu_{s}=T_{s+1}\mkern 3.0mu\overline{\mkern-3.0muB\mkern-1.0mu}\mkern
1.0mu_{s}T_{s}^{-1}$

7: $\mathfrak{V}_{s+1}=V_{s+1},\quad\mkern
3.0mu\overline{\mkern-3.0mu\mathfrak{H}\mkern-1.0mu}\mkern 1.0mu_{s}=\mkern
3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern 1.0mu_{s},\quad\mkern
3.0mu\overline{\mkern-3.0mu\mathfrak{B}\mkern-1.0mu}\mkern 1.0mu_{s}=\mkern
3.0mu\overline{\mkern-3.0muB\mkern-1.0mu}\mkern 1.0mu_{s}$

8: else

9:
$\acute{\mathfrak{R}}_{s}=\mathfrak{V}_{s(j-1)+1}^{\mathsf{T}}\acute{W}_{s},\quad\acute{W}_{s}=\acute{W}_{s}-\mathfrak{V}_{s(j-1)+1}\acute{\mathfrak{R}}_{s}$

10: Compute QR factorization $\acute{W}_{s}=\acute{V}_{s}\acute{T}_{s}$

11: $\mkern 3.0mu\overline{\mkern-3.0mu\mathfrak{B}\mkern-1.0mu}\mkern
1.0mu_{sj}=\left(\begin{array}[]{cc}\mathfrak{H}_{s(j-1)}&0\\\
\eta_{j-1}e_{1}e_{s(j-1)}^{\mathsf{T}}&\mkern
3.0mu\overline{\mkern-3.0muB\mkern-1.0mu}\mkern
1.0mu_{s}\end{array}\right),\quad\mathfrak{T}_{sj+1}=\left(\begin{array}[]{cc}I_{s(j-1)+1}&\acute{\mathfrak{R}}_{s}\\\
0&\acute{T}_{s}\end{array}\right)$

12: $\mkern 3.0mu\overline{\mkern-3.0mu\mathfrak{H}\mkern-1.0mu}\mkern
1.0mu_{sj}=\mathfrak{T}_{sj+1}\mkern
3.0mu\overline{\mkern-3.0mu\mathfrak{B}\mkern-1.0mu}\mkern
1.0mu_{sj}\mathfrak{T}_{sj}^{-1},\quad\mathfrak{V}_{sj+1}=[\mathfrak{V}_{s(j-1)+1},\acute{V}_{s}]$

13: end if

14:end for

15:Compute $y$ such that $\left\lVert\beta e_{1}-\mkern
3.0mu\overline{\mkern-3.0mu\mathfrak{H}\mkern-1.0mu}\mkern
1.0mu_{st}y\right\rVert$ is minimized

16:$x_{st}=x_{0}+\mathfrak{V}_{st}y$

Algorithm 12 illustrates the $s$-step GMRES algorithm in a restart cycle,
which is mathematically equivalent to Algorithm 2. Indeed, the lack of
consensus on the meaning of “$s$-step GMRES” often leads to confusion. In
[187, 141] the new formulation with MPK, BGS and TSQR was called
communication-avoiding GMRES (CA-GMRES). In our context, Algorithm 12 depicts
the $s$-step GMRES algorithm, whereas CA-GMRES refers to the MPK-based
$s$-step GMRES algorithm, namely, a matrix powers kernel is employed in line
3, where $s$ basis vectors can be computed for the same communication cost as
computing 1 basis vector. Here $\mathfrak{T}_{sj}$ is the $sj\times sj$
principal submatrix of $\mathfrak{T}_{sj+1}$. In line 11, $\eta_{j-1}$ denotes
the bottommost value of the $s(j-1)$th column of $\mkern
3.0mu\overline{\mkern-3.0mu\mathfrak{H}\mkern-1.0mu}\mkern 1.0mu_{s(j-1)}$.
For the sake of clarity, we denote by $I_{s}$ the identity matrix of size $s$.
Further, the upper Hessenberg matrix $\mkern
3.0mu\overline{\mkern-3.0mu\mathfrak{H}\mkern-1.0mu}\mkern 1.0mu_{sj}$ in line
12 can be updated in a more efficient way; see [141] for more details.

A full account of TSQR can be found in [67], while Hoemmen’s thesis [141]
contains further discussion of QR algorithms. Parallel version of deflated
GMRES was studied by Wakam and Erhel [276] in 2013, who combined Newton basis
$s$-step GMRES with adaptive augmentation. Similar topics were discussed
independently in [287]; we refer the reader to [276] and [287] for algorithms
and details. In 2017, Imberti and Erhel [143] proposed an $s$-step GMRES
algorithm where the block size $s$ is variable. They suggested to use a
predetermined sequence based on the Fibonacci numbers. Experiments on multiple
GPUs can be found in [284], with a particular focus on orthogonalization
strategies; see also the experiments in [285, 286]. Moreover, polynomial
preconditioning could be used with $s$-step GMRES; see [178, 177] for more
discussion.

###  4.3 Pipelined GMRES

In 2013, Ghysels et al. [113] proposed a pipelined variant of GMRES which only
requires a single nonblocking reduction per iteration and can interleave dot
products with SpMVs to hide communication latency. An important observation is
that (9) can be rewritten via the Pythagorean theorem as

$h_{j+1,j}=\sqrt{\left\lVert
Av_{j}\right\rVert^{2}-\sum_{i=1}^{j}h_{i,j}^{2}}.$ (30)

The CGS process involving (30) could be denoted by CGS-P (see, e.g., [246,
52]), in which one finds that only one reduction is required, whereas two are
required in the traditional CGS-Arnoldi algorithm. This results in the single-
reduction GMRES algorithm. At the $j$th iteration, the orthogonalization
process can be summarized as follows:

  1. 1.

$w_{j}=Av_{j}$;

  2. 2.

Compute $h_{i,j}=(w_{j},v_{i})$ for $i=1,2,\dots,j$ and $(w_{j},w_{j})$ by a
global reduction;

  3. 3.

$h_{j+1,j}=\sqrt{(w_{j},w_{j})-\sum_{i=1}^{j}h_{i,j}^{2}}$;

  4. 4.

$v_{j+1}=(w_{j}-\sum_{i=1}^{j}h_{i,j}v_{i})/h_{j+1,j}$.

In practice for stability reasons, a polynomial basis should be used, and a
reorthogonalization process could be helpful when a breakdown occurs in step 3
(see [113]).

It was observed that $w_{j+1}$ can be obtained by means of an updating
formula. Along with a subtle rearrangement of operations, we obtain a
communication-hiding variant shown in Algorithm 13.

Algorithm 13 Pipelined GMRES

1:$\beta=\left\lVert r_{0}\right\rVert,\quad v_{1}=r_{0}/\beta,\quad
w_{1}=(A-\theta I)v_{1}$

2:for $j=1,2,\dots,n$ do

3: Compute $h_{i,j}=(w_{j},v_{i})$ for $i=1,2,\dots,j$ and
$\varsigma=(w_{j},w_{j})$ using nonblocking reductions

4: $u=Aw_{j}$

5: $h_{j+1,j}=\sqrt{\varsigma-\sum_{i=1}^{j}h_{i,j}^{2}}$

6: $v_{j+1}=(w_{j}-\sum_{i=1}^{j}h_{i,j}v_{i})/h_{j+1,j}$

7: $h_{i,j}=h_{i,j}+\theta$

8: $w_{j+1}=(u-\sum_{i=1}^{j}h_{i,j}w_{i})/h_{j+1,j}$

9:end for

10:Compute $y$ such that $\left\lVert\beta e_{1}-\mkern
3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern 1.0mu_{n}y\right\rVert$ is
minimized

11:$x_{n}=x_{0}+V_{n}y$

The single-stage pipelined GMRES algorithm, here loosely called pipelined
GMRES, can achieve a speedup factor of at most 2. We can see in the first line
that the Newton basis is used to improve numerical stability. In the for-loop,
one finds that $w_{j+1}=(A-\theta I)v_{j+1}$ holds, and thus Newton
polynomials are implicitly involved. If a dot product entails higher latency
than an SpMV, then a greater speedup might be possible by overlapping a
nonblocking reduction with at most $l$ SpMVs, in other words, using deeper
pipelining to hide a dot-product latency of at most $l$ iterations; see [113]
for the algorithm and a detailed derivation. The $l$-stage variant, however,
does not decrease the number of global reductions and can introduce more
source of instability.

Note that pipelining allows us to easily use preconditioners, whereas $s$-step
algorithms would encounter some difficulties (see, e.g, [113, 285]). In 2016,
Sanan et al. [231] presented a pipelined variant of FGMRES, which was also
mentioned in [113], and expressed doubts about the combination of deeper
pipelining with nonlinear preconditioners. The idea of combining the strengths
of the pipelined and $s$-step algorithms was considered by Yamazaki et al.
[285]. It was argued in [178] that polynomial preconditioners would be well
suited to pipelined GMRES. Besides, pipelining was used in [190] to study the
performance variability on high-performance computing systems. We also refer
the reader to [87] for a detailed account of nonblocking Krylov solvers.

###  4.4 Low-sync GMRES

In 2007, Hernández et al. [138] developed three variants of the Arnoldi
process for distributed-memory platforms. In particular, the so-called Arnoldi
with delayed reorthogonalization (ADR) algorithm could reduce synchronization
points by rearranging and deferring operations. Numerical experiments reported
in [138] show that ADR provides a large performance improvement when the
number of processors was sufficiently increased, but suffers from numerical
instability; see [141] for a comparison of ADR and $s$-step strategies.

Recently, Swirydowicz et al. [256] suggested a new variant of MGS-GMRES that
requires only one reduction per iteration. The main idea proposed is to apply
the delayed normalization technique of ADR to the inverse compact WY
representation of MGS. Recall that MGS can be viewed as a composition of rank-
one projectors. Swirydowicz et al. [256] noted the inverse compact WY form of
MGS (ICWY-MGS), which was partly anticipated by Björck’s [32] more than five
decades ago, but laid dormant for a long time. Defining

$L_{n}=\left(\begin{array}[]{cc}L_{n-1}&0\\\
(V_{n-1}^{\mathsf{T}}v_{n})^{\mathsf{T}}&0\end{array}\right),\quad L_{1}=0,$

it follows that

$\hat{P}_{n}=(I-v_{n}v_{n}^{\intercal})\dots(I-v_{1}v_{1}^{\intercal})=I-V_{n}\hat{T}_{n}V_{n}^{\mathsf{T}},\quad\hat{T}_{n}=(I+L_{n})^{-1},$

where $L_{n}$ is a strictly lower triangular matrix in finite precision
arithmetic; see [256]. Only one reduction is needed for ICWY-MGS to perform
the projection. In general, line 7 in Algorithm 1 requires an extra global
reduction for the normalization. Now, let $l_{i,j}$ and $h_{i,j}$ denote
respectively the $(i,j)$ entries of $L_{n}$ and $\mkern
3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern 1.0mu_{n}$. The delayed
normalization technique described in [138] was employed in [256], resulting in
a low-sync GMRES algorithm as shown in Algorithm 14.

Algorithm 14 Low-sync GMRES

1:$\beta=\left\lVert r_{0}\right\rVert,\quad v_{1}=r_{0}/\beta,\quad
h_{1,1}=(w_{1},v_{1}),\quad w_{1}=w_{1}-h_{1,1}v_{1}$

2:for $j=2,3,\dots,n+1$ do

3: $w_{j}=Aw_{j-1}$

4: $l_{j,1:j-1}=(V_{j-1}^{\mathsf{T}}w_{j-1})^{\mathsf{T}},\quad
u=[V_{j-1},w_{j-1}]^{\mathsf{T}}w_{j},\quad h_{j,j-1}=\left\lVert
w_{j-1}\right\rVert$

5: $v_{j}=w_{j-1}/h_{j,j-1}$

6: $l_{j,1:j-1}=l_{j,1:j-1}/h_{j,j-1},\quad u=u/h_{j,j-1},\quad
u_{j}=u_{j}/h_{j,j-1},\quad w_{j}=w_{j}/h_{j,j-1}$

7: $h_{1:j,j}=(I+L_{j})^{-1}u$

8: $w_{j}=w_{j}-V_{j}h_{1:j,j}$

9:end for

10:Compute $y$ such that $\left\lVert\beta e_{1}-\mkern
3.0mu\overline{\mkern-3.0muH\mkern-1.0mu}\mkern 1.0mu_{n}y\right\rVert$ is
minimized

11:$x_{n}=x_{0}+V_{n}y$

We can see that delaying the normalization of $w_{j-1}$ brings additional
scaling (line 6), but can merge the normalization with dot products (line 4).
As a result, only one global reduction per iteration is required.

Several formulations of Gram-Schmidt were studied in [256] and [30] as they
can minimize the number of global reductions. For example, a single-reduction
CGS2 variant can also be derived by means of delayed normalization. Numerical
results reported promising performance of the low-sync Gram-Schmidt and GMRES
algorithms. A good survey of Gram-Schmidt up to 2010 was given by Leon et al.
[162] in 2013; we also refer the interested reader to [23, 22, 52, 53] and
their references for the discussion of block Gram-Schmidt algorithms. On the
other hand, a low-sync GMRES variant based on Householder QR was derived by
Walker [277] in 1988. The single-reduction orthogonalization procedures were
applied to $s$-step and pipelined algorithms by Yamazaki et al. [286],
together with numerical experiments on their performance. The recent survey
article [1] that were intended to provide a road map for mixed-precision
algorithms contains also a perspective on mixed-precision ICWY-MGS.

##  5 Extensions to other problems

This section provides a concise literature review on GMRES solvers for other
related problems. First, a more complicated but interesting problem is to
solve systems with multiple right-hand sides

$AX=B,$ (31)

where the $N\times S$ matrix $B$ is a collection of right-hand sides with
$S\ll N$. It is known that block Krylov algorithms can handle all columns at
once, and thus has a few practical advantages in terms of data locality and
commnication efficiency. Block GMRES was described by Vital [273] in 1990 and
then studied by Simoncini and Gallopoulos [239, 240]; see also [131, 220, 194,
213, 133, 43, 42, 2, 3, 108, 158, 257]. Block GMRES has also been extended to
handle the single right-hand side system (1); see, e.g., [55, 15, 130]. As an
alternative, Jbilou et al. [146] in 1999 proposed a global GMRES algorithm for
solving (31); see also [36, 85].

In 2006, Parks et al. [207] proposed a recycled GMRES algorithm, called GCRO-
DR, for solving a sequence of linear systems

$A^{(i)}x^{(i)}=b^{(i)},$

where the matrix changes slowly from one system to the next, but the
cumulative change over many iterations can be large. On the other hand,
Frommer and Glässner [106] in 1998 considered the shifted coefficient matrix
$A+\sigma I$ and proposed a shifted GMRES algorithm for this case. The
solution of sequences of shifted linear systems

$(A+\sigma_{i}I)x^{(i)}=b$ (32)

was studied in [250, 247, 248] by combining shifted GMRES with subspace
recycling. We refer to the survey by Soodhalter et al. [249] for more details
on recycling strategies. A multi-preconditioned GMRES algorithm was proposed
in [130], which when combined with the flexible approach developed in [230],
could be applied to (32); see [18]. GMRES could have broader applications for
more general shifted systems; see [62, 250, 247, 255, 108] for more details.
Now consider the Sylvester equation

$AX+XD=B,$ (33)

where $A\in\mathbb{R}^{N\times N}$, $D\in\mathbb{R}^{S\times S}$, and
$B\in\mathbb{R}^{N\times S}$. A block GMRES algorithm for solving (33) was
described by Robbé and Sadkane [212] in 2002; see also [238, 292]. It is known
that a family of shifted systems can be interpreted as a Sylvester equation of
the form (33) where $D$ is a diagonal matrix; see [247, 84]. GMRES has also
been extended to solve least square problems [137], tensor equations [57], and
quaternion linear systems [148].

Going back to the original system (1), we point the reader to [101, 29, 107,
252], where the coefficient matrix $A$ has a $2\times 2$ block structure.
There are practical cases where the matrix $A$ is singular or nearly singular.
The article [40] by Brown and Walker in 1997 seems to be the first to study
the behavior of GMRES for solving singular or ill-conditioned systems. Later,
the Drazin-inverse solution of singular linear systems was discussed by Ipsen
and Meyer [145] and Sidi [234]; see also [46, 211, 86, 125, 196, 111]. In
parallel to these continuing efforts, a range restricted GMERS algorithm
geared again toward singular systems was developed by Calvetti et al. [44] in
2000. Many publications on this subject followed; see, e.g., [45, 211, 13,
26].

##  6 Conclusions

The efforts to develop GMRES algorithms span now three and a half decades. In
this paper we have given an overview of GMRES, including basic algorithms,
convergence results, acceleration strategies, parallel algorithms, and a brief
account of extensions. An important area which was not addressed is the
performance of GMRES in real-world applications; see, e.g., [282, 39, 203,
166, 161, 109], a thorough discussion of which is beyond the scope of this
article.

At present, we can see that efficient parallelization of GMRES represents the
frontier for linear systems. Rigorous stability studies of parallel and block
algorithms are still lacking; we refer to [52] as a recent work in this
direction. Furthermore, for the acceleration strategies mentioned above,
preconditioning may offer significant advantages in a variety of contexts, and
thus should play a role in solving challenging problems. There is strong
evidence that much effort should be spent in developing mixed-precision
parallel GMRES algorithms, along with scalable preconditioners suitable for
highly parallel environments. Describing their finite-precision behavior would
be a quite difficult task.

Another area of interest in numerical linear algebra is the design of
randomized algorithms; see [179] for a detailed survey. We highlight the work
by Balabanov and Grigori [19], where a randomized Gram-Schmidt algorithm was
developed and applied to GMRES.

## References

  * [1] A. Abdelfattah, H. Anzt, E. G. Boman, E. C. Carson, T. Cojean, J. Dongarra, A. Fox, M. Gates, N. J. Higham, X. S. Li, J. Loe, P. Luszczek, S. Pranesh, S. Rajamanickam, T. Ribizel, B. F. Smith, K. Swirydowicz, S. Thomas, S. Tomov, Y. M. Tsai, and U. M. Yang.  A survey of numerical linear algebra methods utilizing mixed-precision arithmetic.  Int. J. High Perform. Comput. Appl., 35(4):344–369, 2021.
  * [2] E. Agullo, L. Giraud, and Y.-F. Jing.  Block GMRES method with inexact breakdowns and deflated restarting.  SIAM J. Matrix Anal. Appl., 35(4):1625–1651, 2014.
  * [3] H. Al Daas, L. Grigori, P. Hénon, and P. Ricoux.  Enlarged GMRES for solving linear systems with one or multiple right-hand sides.  IMA J. Numer. Anal., 39(4):1924–1956, 2019.
  * [4] J. I. Aliaga, H. Anzt, T. Grützmacher, E. S. Quintana-Ortí, and A. E. Tomás.  Compressed basis GMRES on high performance GPUs, 2020.
  * [5] P. Amestoy, A. Buttari, N. J. Higham, J.-Y. L’Excellent, T. Mary, and B. Vieublé.  Five-precision GMRES-based iterative refinement, 2021.
  * [6] M. Arioli and I. S. Duff.  Using FGMRES to obtain backward stability in mixed precision.  Electron. Trans. Numer. Anal., 33:31–44, 2009.
  * [7] M. Arioli, I. S. Duff, S. Gratton, and S. Pralet.  A note on GMRES preconditioned by a perturbed $LDL^{T}$ decomposition with static pivoting.  SIAM J. Sci. Comput., 29(5):2024–2044, 2007.
  * [8] M. Arioli, V. Pták, and Z. Strakoš.  Krylov sequences of maximal length and convergence of GMRES.  BIT, 38(4):636–643, 1998.
  * [9] W. E. Arnoldi.  The principle of minimized iterations in the solution of the matrix eigenvalue problem.  Q. Appl. Math., 9(1):17–29, 1951.
  * [10] O. Axelsson.  Conjugate gradient type methods for unsymmetric and inconsistent systems of linear equations.  Linear Algebra Appl., 29:1–16, 1980.
  * [11] E. H. Ayachour.  A fast implementation for GMRES method.  J. Comput. Appl. Math., 159(2):269–283, 2003.
  * [12] J. Baglama, D. Calvetti, G. H. Golub, and L. Reichel.  Adaptively preconditioned GMRES algorithms.  SIAM J. Sci. Comput., 20(1):243–269, 1998.
  * [13] J. Baglama and L. Reichel.  Augmented GMRES-type methods.  Numer. Linear Algebra Appl., 14(4):337–350, 2007.
  * [14] Z. Bai, D. Hu, and L. Reichel.  A Newton basis GMRES implementation.  IMA J. Numer. Anal., 14(4):563–581, 1994.
  * [15] A. H. Baker, J. M. Dennis, and E. R. Jessup.  On improving linear solver performance: A block variant of GMRES.  SIAM J. Sci. Comput., 27(5):1608–1626, 2006.
  * [16] A. H. Baker, E. R. Jessup, and T. V. Kolev.  A simple strategy for varying the restart parameter in GMRES($m$).  J. Comput. Appl. Math., 230(2):751–761, 2009.
  * [17] A. H. Baker, E. R. Jessup, and T. A. Manteuffel.  A technique for accelerating the convergence of restarted GMRES.  SIAM J. Matrix Anal. Appl., 26(4):962–984, 2005.
  * [18] T. Bakhos, P. K. Kitanidis, S. Ladenheim, A. K. Saibaba, and D. B. Szyld.  Multipreconditioned GMRES for shifted systems.  SIAM J. Sci. Comput., 39(5):S222–S247, 2017.
  * [19] O. Balabanov and L. Grigori.  Randomized Gram-Schmidt process with application to GMRES, 2020\.
  * [20] J. Ballani and L. Grasedyck.  A projection method to solve linear systems in tensor format.  Numer. Linear Algebra Appl., 20(1):27–43, 2013.
  * [21] G. Ballard, E. C. Carson, J. W. Demmel, M. Hoemmen, N. Knight, and O. Schwartz.  Communication lower bounds and optimal algorithms for numerical linear algebra.  Acta Numer., 23:1–155, 2014.
  * [22] J. L. Barlow.  Block modified Gram-Schmidt algorithms and their analysis.  SIAM J. Matrix Anal. Appl., 40(4):1257–1290, 2019.
  * [23] J. L. Barlow and A. Smoktunowicz.  Reorthogonalized block classical Gram-Schmidt.  Numer. Math., 123(3):395–423, 2013.
  * [24] B. Beckermann.  Numerical range, GMRES and Faber polynomials.  C. R. Math., 340(11):855–860, 2005.  (in French).
  * [25] B. Beckermann, S. A. Goreinov, and E. E. Tyrtyshnikov.  Some remarks on the Elman estimate for GMRES.  SIAM J. Matrix Anal. Appl., 27(3):772–778, 2005.
  * [26] M. Bellalij, L. Reichel, and H. Sadok.  Some properties of range restricted GMRES methods.  J. Comput. Appl. Math., 290:310–318, 2015.
  * [27] M. Benzi.  Preconditioning techniques for large linear systems: A survey.  J. Comput. Phys., 182(2):418–477, 2002.
  * [28] M. Benzi.  Some uses of the field of values in numerical analysis.  Boll. Unione. Mat. Ital., 14(1):159–177, 2021.
  * [29] M. Benzi, G. H. Golub, and J. Liesen.  Numerical solution of saddle point problems.  Acta Numer., 14:1–137, 2005.
  * [30] D. Bielich, J. Langou, S. Thomas, K. Swirydowicz, I. Yamazaki, and E. G. Boman.  Low-synch Gram-Schmidt with delayed reorthogonalization for Krylov solvers, 2021.
  * [31] D. Bindel, J. W. Demmel, W. Kahan, and O. Marques.  On computing Givens rotations reliably and efficiently.  ACM Trans. Math. Softw., 28(2):206–238, 2002.
  * [32] Å. Björck.  Solving linear least squares problems by Gram-Schmidt orthogonalization.  BIT, 7(1):1–21, 1967.
  * [33] J. Blechta.  Stability of linear GMRES convergence with respect to compact perturbations.  SIAM J. Matrix Anal. Appl., 42(1):436–447, 2021.
  * [34] A. Bouras and V. Frayssé.  A relaxation strategy for inexact matrix-vector products for Krylov methods.  Technical report, CERFACS, 2000.
  * [35] A. Bouras and V. Frayssé.  Inexact matrix-vector products in Krylov methods for solving linear systems: A relaxation strategy.  SIAM J. Matrix Anal. Appl., 26(3):660–678, 2005.
  * [36] R. Bouyouli, K. Jbilou, R. Sadaka, and H. Sadok.  Convergence properties of some block Krylov subspace methods for multiple linear systems.  J. Comput. Appl. Math., 196(2):498–511, 2006.
  * [37] C. Brezinski, M. Redivo-Zaglia, and Y. Saad.  Shanks sequence transformations and Anderson acceleration.  SIAM Rev., 60(3):646–669, 2018.
  * [38] P. N. Brown.  A theoretical comparison of the Arnoldi and GMRES algorithms.  SIAM J. Sci. Stat. Comput., 12(1):58–78, 1991.
  * [39] P. N. Brown and A. C. Hindmarsh.  Reduced storage matrix methods in stiff ODE systems.  Appl. Math. Comput., 31:40–91, 1989.
  * [40] P. N. Brown and H. F. Walker.  GMRES on (nearly) singular systems.  SIAM J. Matrix Anal. Appl., 18(1):37–51, 1997.
  * [41] K. Burrage and J. Erhel.  On the performance of various adaptive preconditioned GMRES strategies.  Numer. Linear Algebra Appl., 5(2):101–121, 1998.
  * [42] H. Calandra, S. Gratton, R. Lago, X. Vasseur, and L. M. Carvalho.  A modified block flexible GMRES method with deflation at each iteration for the solution of non-Hermitian linear systems with multiple right-hand sides.  SIAM J. Sci. Comput., 35(5):S345–S367, 2013.
  * [43] H. Calandra, S. Gratton, J. Langou, X. Pinel, and X. Vasseur.  Flexible variants of block restarted GMRES methods with application to geophysics.  SIAM J. Sci. Comput., 34(2):A714–A736, 2012.
  * [44] D. Calvetti, B. Lewis, and L. Reichel.  GMRES-type methods for inconsistent systems.  Linear Algebra Appl., 316(1):157–169, 2000.
  * [45] D. Calvetti, B. Lewis, and L. Reichel.  On the choice of subspace for iterative methods for linear discrete ill-posed problems.  Int. J. Appl. Math. Comput. Sci., 11(5):1069–1092, 2001.
  * [46] D. Calvetti, B. Lewis, and L. Reichel.  On the regularizing properties of the GMRES method.  Numer. Math., 91(4):605–625, 2002.
  * [47] S. L. Campbell, I. C. F. Ipsen, C. T. Kelley, and C. D. Meyer.  GMRES and the minimal polynomial.  BIT, 36(4):664–675, 1996.
  * [48] Z.-H. Cao and X.-Y. Yu.  A note on weighted FOM and GMRES for solving nonsymmetric linear systems.  Appl. Math. Comput., 151(3):719–727, 2004.
  * [49] E. C. Carson.  Communication-Avoiding Krylov Subspace Methods in Theory and Practice.  PhD thesis, UC Berkeley, 2015.
  * [50] E. C. Carson and N. J. Higham.  A new analysis of iterative refinement and its application to accurate solution of ill-conditioned sparse linear systems.  SIAM J. Sci. Comput., 39(6):A2834–A2856, 2017.
  * [51] E. C. Carson and N. J. Higham.  Accelerating the solution of linear systems by iterative refinement in three precisions.  SIAM J. Sci. Comput., 40(2):A817–A847, 2018.
  * [52] E. C. Carson, K. Lund, and M. Rozložník.  The stability of block variants of classical Gram-Schmidt.  SIAM J. Matrix Anal. Appl., 42(3):1365–1380, 2021.
  * [53] E. C. Carson, K. Lund, M. Rozložník, and S. Thomas.  An overview of block Gram-Schmidt methods and their stability properties, 2020.
  * [54] E. C. Carson and Z. Strakoš.  On the cost of iterative computations.  Philos. Trans. R. Soc. A, 378(2166):20190050, 2020.
  * [55] A. Chapman and Y. Saad.  Deflated and augmented Krylov subspace techniques.  Numer. Linear Algebra Appl., 4(1):43–66, 1997.
  * [56] G. Chen and Z. Jia.  Theoretical and numerical comparisons of GMRES and WZ-GMRES.  Comput. Math. Appl., 47(8):1335–1350, 2004.
  * [57] Z. Chen and L. Lu.  A projection method and Kronecker product preconditioner for solving Sylvester tensor equations.  Sci. China Math., 55(6):1281–1292, 2012.
  * [58] A. T. Chronopoulos and S. K. Kim.  $s$-step Orthomin and GMRES implemented on parallel computers.  Technical report, University of Minnesota, 1990.
  * [59] J. K. Cullum and A. Greenbaum.  Relations between Galerkin and norm-minimizing iterative methods for solving linear systems.  SIAM J. Matrix Anal. Appl., 17(2):223–247, 1996.
  * [60] R. D. da Cunha and T. Hopkins.  A parallel implementation of the restarted GMRES iterative algorithm for nonsymmetric systems of linear equations.  Adv. Comput. Math., 2(3):261–277, 1994.
  * [61] J. W. Daniel, W. B. Gragg, L. Kaufman, and G. W. Stewart.  Reorthogonalization and stable algorithms for updating the Gram-Schmidt QR factorization.  Math. Comput., 30(136):772–795, 1976.
  * [62] D. Darnell, R. B. Morgan, and W. Wilcox.  Deflated GMRES for systems with multiple shifts and multiple right-hand sides.  Linear Algebra Appl., 429(10):2415–2434, 2008.
  * [63] E. de Sturler.  A parallel variant of GMRES($m$).  In Proceedings of the 13th IMACS World Congress on Computation and Applied Mathematics, pages 682–683, Dublin, Ireland, 1991. Criterion Press.
  * [64] E. de Sturler.  Nested Krylov methods based on GCR.  J. Comput. Appl. Math., 67(1):15–41, 1996.
  * [65] E. de Sturler.  Truncation strategies for optimal Krylov subspace methods.  SIAM J. Numer. Anal., 36(3):864–889, 1999.
  * [66] E. de Sturler and H. A. van der Vorst.  Reducing the effect of global communication in GMRES($m$) and CG on parallel distributed memory computers.  Appl. Numer. Math., 18(4):441–459, 1995.
  * [67] J. W. Demmel, L. Grigori, M. Hoemmen, and J. Langou.  Communication-optimal parallel and sequential QR and LU factorizations.  SIAM J. Sci. Comput., 34(1):A206–A239, 2012.
  * [68] J. W. Demmel, M. T. Heath, and H. A. van der Vorst.  Parallel numerical linear algebra.  Acta Numer., 2:111–197, 1993.
  * [69] J. W. Demmel, M. Hoemmen, M. Mohiyuddin, and K. Yelick.  Avoiding communication in computing Krylov subspaces.  Technical report, UC Berkeley, 2007.
  * [70] J. W. Demmel, M. Hoemmen, M. Mohiyuddin, and K. Yelick.  Avoiding communication in sparse matrix computations.  In Proceedings of the International Parallel and Distributed Processing Symposium, pages 1–12, Miami, FL, USA, 2008. IEEE.
  * [71] G. R. di Brozolo and Y. Robert.  Parallel conjugate gradient-like algorithms for solving sparse nonsymmetric linear systems on a vector multiprocessor.  Parallel Comput., 11(2):223–239, 1989.
  * [72] S. V. Dolgov.  TT-GMRES: Solution to a linear system in the structured tensor format.  Russ. J. Numer. Anal. Math. Model., 28(2):149–172, 2013.
  * [73] J. Drkošová, A. Greenbaum, M. Rozložník, and Z. Strakoš.  Numerical stability of GMRES.  BIT, 35(3):309–330, 1995.
  * [74] K. Du, J. Duintjer Tebbens, and G. Meurant.  Any admissible harmonic Ritz value set is possible for GMRES.  Electron. Trans. Numer. Anal., 47:37–56, 2017.
  * [75] J. Duintjer Tebbens and G. Meurant.  Any Ritz value behavior is possible for Arnoldi and for GMRES.  SIAM J. Matrix Anal. Appl., 33(3):958–978, 2012.
  * [76] J. Duintjer Tebbens and G. Meurant.  Prescribing the behavior of early terminating GMRES and Arnoldi iterations.  Numer. Algorithms, 65(1):69–90, 2014.
  * [77] J. Duintjer Tebbens and G. Meurant.  On the residual norms, the Ritz values and the harmonic Ritz values that can be generated by restarted GMRES.  Numer. Algorithms, 84(4):1329–1352, 2020.
  * [78] J. Duintjer Tebbens, G. Meurant, H. Sadok, and Z. Strakoš.  On investigating GMRES convergence using unitary matrices.  Linear Algebra Appl., 450:83–107, 2014.
  * [79] M. Eiermann.  Fields of values and iterative methods.  Linear Algebra Appl., 180:167–197, 1993.
  * [80] M. Eiermann and O. G. Ernst.  Geometric aspects of the theory of Krylov subspace methods.  Acta Numer., 10:251–312, 2001.
  * [81] M. Eiermann, O. G. Ernst, and O. Schneider.  Analysis of acceleration strategies for restarted minimal residual methods.  J. Comput. Appl. Math., 123(1):261–292, 2000.
  * [82] T. Eirola and O. Nevanlinna.  Accelerating with rank-one updates.  Linear Algebra Appl., 121:511–520, 1989.
  * [83] S. C. Eisenstat, H. C. Elman, and M. H. Schultz.  Variational iterative methods for nonsymmetric systems of linear equations.  SIAM J. Numer. Anal., 20(2):345–357, 1983.
  * [84] L. Elbouyahyaoui, M. Heyouni, A. Tajaddini, and F. Saberi-Movahed.  On restarted and deflated block FOM and GMRES methods for sequences of shifted linear systems.  Numer. Algorithms, 87(3):1257–1299, 2021.
  * [85] L. Elbouyahyaoui, A. Messaoudi, and H. Sadok.  Algebraic properties of the block GMRES and block Arnoldi methods.  Electron. Trans. Numer. Anal., 33:207–220, 2009.
  * [86] L. Eldén and V. Simoncini.  Solving ill-posed linear systems with GMRES and a singular preconditioner.  SIAM J. Matrix Anal. Appl., 33(4):1369–1394, 2012.
  * [87] P. R. Eller.  Scalable Non-blocking Krylov Solvers for Extreme-scale Computing.  PhD thesis, University of Illinois at Urbana-Champaign, 2019.
  * [88] H. C. Elman.  Iterative Methods for Large, Sparse, Nonsymmetric Systems of Linear Equations.  PhD thesis, Technische Universität Hamburg, 1982.
  * [89] M. Embree.  How descriptive are GMRES convergence bounds?  Technical report, Oxford University Computing Laboratory, 1999.
  * [90] M. Embree.  The tortoise and the hare restart GMRES.  SIAM Rev., 45(2):259–266, 2003.
  * [91] M. Embree, J. A. Loe, and R. B. Morgan.  Polynomial preconditioned Arnoldi with stability control.  SIAM J. Sci. Comput., 43(1):A1–A25, 2021.
  * [92] M. Embree, R. B. Morgan, and H. V. Nguyen.  Weighted inner products for GMRES and GMRES-DR.  SIAM J. Sci. Comput., 39(5):S610–S632, 2017.
  * [93] J. Erhel.  A parallel GMRES version for general sparse matrices.  Electron. Trans. Numer. Anal., 3:160–176, 1995.
  * [94] J. Erhel, K. Burrage, and B. Pohl.  Restarted GMRES preconditioned by deflation.  J. Comput. Appl. Math., 69(2):303–318, 1996.
  * [95] Y. A. Erlangga and R. Nabben.  Multilevel projection-based nested Krylov iteration for boundary value problems.  SIAM J. Sci. Comput., 30(3):1572–1595, 2008.
  * [96] Y. A. Erlangga and R. Nabben.  On a multilevel Krylov method for the Helmholtz equation preconditioned by shifted Laplacian.  Electron. Trans. Numer. Anal., 31:403–424, 2008.
  * [97] A. Essai.  Weighted FOM and GMRES for solving nonsymmetric linear systems.  Numer. Algorithms, 18(3-4):277–292, 1998.
  * [98] V. Faber, A. Greenbaum, and D. E. Marshall.  The polynomial numerical hulls of Jordan blocks and related matrices.  Linear Algebra Appl., 374:231–246, 2003.
  * [99] V. Faber, W. D. Joubert, E. Knill, and T. A. Manteuffel.  Minimal residual method stronger than polynomial preconditioning.  SIAM J. Matrix Anal. Appl., 17(4):707–729, 1996.
  * [100] V. Faber, J. Liesen, and P. Tichý.  Properties of worst-case GMRES.  SIAM J. Matrix Anal. Appl, 34(4):1500–1519, 2013.
  * [101] B. Fischer, A. Ramage, D. J. Silvester, and A. J. Wathen.  Minimum residual methods for augmented systems.  BIT, 38(3):527–543, 1998.
  * [102] G. E. Forsythe.  On the asymptotic directions of the $s$-dimensional optimum gradient method.  Numer. Math., 11(1):57–76, 1968.
  * [103] V. Frayssé, L. Giraud, and S. Gratton.  Algorithm 881: A set of flexible GMRES routines for real and complex arithmetics on high-performance computers.  ACM Trans. Math. Softw., 35(2):13:1–13:12, 2008.
  * [104] V. Frayssé, L. Giraud, S. Gratton, and J. Langou.  Algorithm 842: A set of GMRES routines for real and complex arithmetics on high performance computers.  ACM Trans. Math. Softw., 31(2):228–238, 2005.
  * [105] R. W. Freund, G. H. Golub, and N. M. Nachtigal.  Iterative solution of linear systems.  Acta Numer., 1:57–100, 1992.
  * [106] A. Frommer and U. Glässner.  Restarted GMRES for shifted linear systems.  SIAM J. Sci. Comput., 19(1):15–26, 1998.
  * [107] A. Frommer, B. Jacob, K. Kahl, C. Wyss, and I. Zwaan.  Krylov type methods for linear systems exploiting properties of the quadratic numerical range.  Electron. Trans. Numer. Anal., 53:541–561, 2020.
  * [108] A. Frommer, K. Lund, and D. B. Szyld.  Block Krylov subspace methods for functions of matrices II: Modified block FOM.  SIAM J. Matrix Anal. Appl., 41(2):804–837, 2020.
  * [109] M. J. Gander, I. G. Graham, and E. A. Spence.  Applying GMRES to the Helmholtz equation with shifted Laplacian preconditioning: What is the largest shift for which wavenumber-independent convergence is guaranteed?  Numer. Math., 131(3):567–614, 2015.
  * [110] A. Gaul, M. H. Gutknecht, J. Liesen, and R. Nabben.  A framework for deflated and augmented Krylov subspace methods.  SIAM J. Matrix Anal. Appl., 34(2):495–518, 2013.
  * [111] S. Gazzola, S. Noschese, P. Novati, and L. Reichel.  Arnoldi decomposition, GMRES, and preconditioning for linear discrete ill-posed problems.  Appl. Numer. Math., 142:102–121, 2019.
  * [112] A. Ghai, C. Lu, and X. Jiao.  A comparison of preconditioned Krylov subspace methods for large-scale nonsymmetric linear systems.  Numer. Linear Algebra Appl., 26(1):e2215, 2019.
  * [113] P. Ghysels, T. J. Ashby, K. Meerbergen, and W. Vanroose.  Hiding global communication latency in the GMRES algorithm on massively parallel machines.  SIAM J. Sci. Comput., 35(1):C48–C71, 2013.
  * [114] L. Giraud, S. Gratton, and J. Langou.  Convergence in backward error of relaxed GMRES.  SIAM J. Sci. Comput., 29(2):710–728, 2007.
  * [115] L. Giraud, S. Gratton, X. Pinel, and X. Vasseur.  Flexible GMRES with deflated restarting.  SIAM J. Sci. Comput., 32(4):1858–1878, 2010.
  * [116] L. Giraud, J. Langou, and M. Rozložník.  The loss of orthogonality in the Gram-Schmidt orthogonalization process.  Comput. Math. Appl., 50(7):1069–1075, 2005.
  * [117] L. Giraud, J. Langou, M. Rozložník, and J. van den Eshof.  Rounding error analysis of the classical Gram-Schmidt orthogonalization process.  Numer. Math., 101(1):87–100, 2005.
  * [118] G. H. Golub and C. F. Van Loan.  Matrix Computations.  Johns Hopkins University Press, Baltimore, MD, 4th edition, 2013.
  * [119] S. Goossens and D. Roose.  Ritz and harmonic Ritz values and the convergence of FOM and GMRES.  Numer. Linear Algebra Appl., 6(4):281–293, 1999.
  * [120] S. Gratton, E. Simon, D. Titley-Peloquin, and P. Toint.  Exploiting variable precision in GMRES, 2019.
  * [121] A. Greenbaum.  Iterative Methods for Solving Linear Systems.  SIAM, Philadelphia, PA, 1997.
  * [122] A. Greenbaum.  Generalizations of the field of values useful in the study of polynomial functions of a matrix.  Linear Algebra Appl., 347(1):233–249, 2002.
  * [123] A. Greenbaum.  Some theoretical results derived from polynomial numerical hulls of Jordan blocks.  Electron. Trans. Numer. Anal., 18:81–90, 2004.
  * [124] A. Greenbaum and L. Gurvits.  Max-min properties of matrix factor norms.  SIAM J. Sci. Comput., 15(2):348–358, 1994.
  * [125] A. Greenbaum, F. Kyanfar, and A. Salemi.  On the convergence rate of DGMRES.  Linear Algebra Appl., 552:219–238, 2018.
  * [126] A. Greenbaum, V. Pták, and Z. Strakoš.  Any nonincreasing convergence curve is possible for GMRES.  SIAM J. Matrix Anal. Appl., 17(3):465–469, 1996.
  * [127] A. Greenbaum, M. Rozložník, and Z. Strakoš.  Numerical behaviour of the modified Gram-Schmidt GMRES implementation.  BIT, 37(3):706–719, 1997.
  * [128] A. Greenbaum and Z. Strakoš.  Matrices that generate the same Krylov residual spaces.  In Recent Advances in Iterative Methods, pages 95–118, New York, NY, 1994. Springer.
  * [129] A. Greenbaum and L. N. Trefethen.  GMRES/CR and Arnoldi/Lanczos as matrix approximation problems.  SIAM J. Sci. Comput., 15(2):359–368, 1994.
  * [130] C. Greif, T. Rees, and D. B. Szyld.  GMRES with multiple preconditioners.  SeMA J., 74(2):213–231, 2017.
  * [131] G.-D. Gu and Z.-H. Cao.  A block GMRES method augmented with eigenvectors.  Appl. Math. Comput., 121(2):271–289, 2001.
  * [132] M. H. Gutknecht.  Spectral deflation in Krylov solvers: A theory of coordinate space based methods.  Electron. Trans. Numer. Anal., 39:156–185, 2012.
  * [133] M. H. Gutknecht and T. Schmelzer.  The block grade of a block Krylov space.  Linear Algebra Appl., 430(1):174–185, 2009.
  * [134] S. Güttel and J. Pestana.  Some observations on weighted GMRES.  Numer. Algorithms, 67(4):733–752, 2014.
  * [135] R. Haelterman, J. Degroote, D. Van Heule, and J. Vierendeels.  On the similarities between the quasi-Newton inverse least squares method and GMRes.  SIAM J. Numer. Anal., 47(6):4660–4679, 2010.
  * [136] R. Haelterman, B. Lauwens, F. Van Utterbeeck, H. Bruyninckx, and J. Vierendeels.  On the similarities between the quasi-Newton least squares method and GMRes.  J. Comput. Appl. Math., 273(6):25–28, 2015.
  * [137] K. Hayami, J.-F. Yin, and T. Ito.  GMRES methods for least squares problems.  SIAM J. Matrix Anal. Appl., 31(5):2400–2430, 2010.
  * [138] V. Hernández, J. E. Román, and A. Tomás.  Parallel Arnoldi eigensolvers with enhanced scalability via global communications rearrangement.  Parallel Comput., 33(7):521–540, 2007.
  * [139] N. J. Higham, S. Pranesh, and M. Zounon.  Squeezing a matrix into half precision, with an application to solving linear systems.  SIAM J. Sci. Comput., 41(4):A2536–A2551, 2019.
  * [140] M. Hochbruck and C. Lubich.  Error analysis of Krylov methods in a nutshell.  SIAM J. Sci. Comput., 19(2):695–701, 1998.
  * [141] M. Hoemmen.  Communication-Avoiding Krylov Subspace Methods.  PhD thesis, UC Berkeley, 2010.
  * [142] A. Imakura, R.-C. Li, and S.-L. Zhang.  Locally optimal and heavy ball GMRES methods.  Jpn. J. Ind. Appl. Math., 33(2):471–499, 2016.
  * [143] D. Imberti and J. Erhel.  Varying the $s$ in your $s$-step GMRES.  Electron. Trans. Numer. Anal., 47:206–230, 2017.
  * [144] I. C. F. Ipsen.  Expressions and bounds for the GMRES residual.  BIT, 40(3):524–535, 2000.
  * [145] I. C. F. Ipsen and C. D. Meyer.  The idea behind Krylov methods.  Am. Math. Mon., 105(10):889–899, 1998.
  * [146] K. Jbilou, A. Messaoudi, and H. Sadok.  Global FOM and GMRES algorithms for matrix equations.  Appl. Numer. Math., 31(1):49–63, 1999.
  * [147] K. C. Jea and D. M. Young.  On the simplification of generalized conjugate-gradient methods for nonsymmetrizable linear systems.  Linear Algebra Appl., 52-53:399–417, 1983.
  * [148] Z. Jia and M. K. Ng.  Structure preserving quaternion generalized minimal residual method.  SIAM J. Matrix Anal. Appl., 42(2):616–634, 2021.
  * [149] P. Jiránek and M. Rozložník.  Adaptive version of simpler GMRES.  Numer. Algorithms, 53(1):93–112, 2010.
  * [150] P. Jiránek, M. Rozložník, and M. H. Gutknecht.  How to make simpler GMRES and GCR more stable.  SIAM J. Matrix Anal. Appl., 30(4):1483–1499, 2008.
  * [151] W. D. Joubert.  On the convergence behavior of the restarted GMRES algorithm for solving nonsymmetric linear systems.  Numer. Linear Algebra Appl., 1(5):427–447, 1994.
  * [152] W. D. Joubert.  A robust GMRES-based adaptive polynomial preconditioning algorithm for nonsymmetric linear systems.  SIAM J. Sci. Comput., 15(2):427–439, 1994.
  * [153] W. D. Joubert and G. F. Carey.  Parallelizable restarted iterative methods for nonsymmetric linear systems. Part I: Theory.  Int. J. Comput. Math., 44(1-4):243–267, 1992.
  * [154] W. D. Joubert and G. F. Carey.  Parallelizable restarted iterative methods for nonsymmetric linear systems. Part II: Parallel implementation.  Int. J. Comput. Math., 44(1-4):269–290, 1992.
  * [155] W. D. Joubert and T. A. Manteuffel.  Iterative methods for nonsymmetric linear systems.  In Iterative Methods for Large Linear Systems, pages 149–171, Austin, TX, USA, 1990. Academic Press.
  * [156] I. M. Khabaza.  An iterative least-square method suitable for solving large sparse matrices.  Comput. J., 6(2):202–206, 1963.
  * [157] S. A. Kharchenko and A. Y. Yeremin.  Eigenvalue translation based preconditioners for the GMRES($k$) method.  Numer. Linear Algebra Appl., 2(1):51–77, 1995.
  * [158] M. Kubínová and K. M. Soodhalter.  Admissible and attainable convergence behavior of block Arnoldi and GMRES.  SIAM J. Matrix Anal. Appl., 41(2):464–486, 2020.
  * [159] C. Le Calvez and B. Molina.  Implicitly restarted and deflated GMRES.  Numer. Algorithms, 21(1-4):261–285, 1999.
  * [160] R. B. Lehoucq and D. C. Sorensen.  Deflation techniques for an implicitly restarted Arnoldi iteration.  SIAM J. Matrix Anal. Appl., 17(4):789–821, 1996.
  * [161] J.-F. Lemieux, B. Tremblay, S. Thomas, J. Sedláček, and L. A. Mysak.  Using the preconditioned Generalized Minimum RESidual (GMRES) method to solve the sea-ice momentum equation.  J. Geophys. Res.-Oceans, 113(C10):C10004, 2008.
  * [162] S. J. Leon, Å. Björck, and W. Gander.  Gram-Schmidt orthogonalization: 100 years and more.  Numer. Linear Algebra Appl., 20(3):492–532, 2013.
  * [163] J. Liesen.  Computable convergence bounds for GMRES.  SIAM J. Matrix Anal. Appl., 21(3):882–903, 2000.
  * [164] J. Liesen, M. Rozložník, and Z. Strakoš.  Least squares residuals and minimal residual methods.  SIAM J. Sci. Comput., 23(5):1503–1525, 2002.
  * [165] J. Liesen and Z. Strakoš.  Convergence of GMRES for tridiagonal Toeplitz matrices.  SIAM J. Matrix Anal. Appl., 26(1):233–251, 2004.
  * [166] J. Liesen and Z. Strakoš.  GMRES convergence analysis for a convection-diffusion model problem.  SIAM J. Sci. Comput., 26(6):1989–2009, 2005.
  * [167] J. Liesen and Z. Strakoš.  Krylov Subspace Methods: Principles and Analysis.  Oxford University Press, Oxford, 2013.
  * [168] J. Liesen and P. Tichý.  Convergence analysis of Krylov subspace methods.  GAMM-Mitteilungen, 27(2):153–173, 2004.
  * [169] J. Liesen and P. Tichý.  The worst-case GMRES for normal matrices.  BIT, 44(1):79–98, 2004.
  * [170] J. Liesen and P. Tichý.  On best approximations of polynomials in matrices in the matrix 2-norm.  SIAM J. Matrix Anal. Appl., 31(2):853–863, 2009.
  * [171] J. Liesen and P. Tichý.  The field of values bound on ideal GMRES, 2012.
  * [172] N. Lindquist, P. Luszczek, and J. Dongarra.  Improving the performance of the GMRES method using mixed-precision techniques.  In The Smoky Mountains Computational Sciences and Engineering Conference, pages 51–66, Oak Ridge, TN, USA, 2020. Springer.
  * [173] N. Lindquist, P. Luszczek, and J. Dongarra.  Accelerating restarted GMRES with mixed precision arithmetic.  IEEE Trans. Parallel Distrib. Syst., (in press).
  * [174] Q. Liu, R. B. Morgan, and W. Wilcox.  Polynomial preconditioned GMRES and GMRES-DR.  SIAM J. Sci. Comput., 37(5):S407–S428, 2015.
  * [175] J. A. Loe, C. A. Glusa, I. Yamazaki, E. G. Boman, and S. Rajamanickam.  Experimental evaluation of multiprecision strategies for GMRES on GPUs, 2021.
  * [176] J. A. Loe and R. B. Morgan.  Toward efficient and stable polynomial preconditioning for GMRES, 2019\.
  * [177] J. A. Loe, H. K. Thornquist, and E. G. Boman.  Polynomial preconditioned GMRES to reduce communication in parallel computing, 2019.
  * [178] J. A. Loe, H. K. Thornquist, and E. G. Boman.  Polynomial preconditioned GMRES in Trilinos: Practical considerations for high-performance computing.  In Proceedings of the SIAM Conference on Parallel Processing for Scientific Computing, pages 35–45, Seattle, WA, USA, 2020. SIAM.
  * [179] P.-G. Martinsson and T. Joel A.  Randomized numerical linear algebra: Foundations and algorithms.  Acta Numer., 29:403–572, 2020.
  * [180] M. Matinfar, H. Zareamoghaddam, M. Eslami, and M. Saeidy.  GMRES implementations and residual smoothing techniques for solving ill-posed linear systems.  Comput. Math. Appl., 63(1):1–13, 2012.
  * [181] L. C. McInnes, B. F. Smith, H. Zhang, and R. T. Mills.  Hierarchical Krylov and nested Krylov methods for extreme-scale computing.  Parallel Comput., 40(1):17–31, 2014.
  * [182] G. Meurant.  The complete stagnation of GMRES for $n\leq 4$.  Electron. Trans. Numer. Anal., 39:75–101, 2012.
  * [183] G. Meurant.  GMRES and the Arioli, Pták, and Strakoš parametrization.  BIT, 52(3):687–702, 2012.
  * [184] G. Meurant.  Necessary and sufficient conditions for GMRES complete and partial stagnation.  Appl. Numer. Math., 75:100–107, 2014.
  * [185] G. Meurant.  The coefficients of the FOM and GMRES residual polynomials.  SIAM J. Matrix Anal. Appl., 38(1):96–117, 2017.
  * [186] G. Meurant and J. Duintjer Tebbens.  Krylov Methods for Nonsymmetric Linear Systems: From Theory to Computations.  Springer, Cham, 2020.
  * [187] M. Mohiyuddin, M. Hoemmen, J. W. Demmel, and K. Yelick.  Minimizing communication in sparse matrix solvers.  In Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, pages 1–12, Portland, OR, USA, 2009. IEEE.
  * [188] C. B. Moler.  Iterative refinement in floating point.  J. ACM, 14(2):316–321, 1967.
  * [189] I. Moret.  A note on the superlinear convergence of GMRES.  SIAM J. Numer. Anal., 34(2):513–516, 1997.
  * [190] H. Morgan, P. Sanan, M. Knepley, and R. T. Mills.  Understanding performance variability in standard and pipelined parallel Krylov solvers.  Int. J. High Perform. Comput. Appl., 35(1):47–59, 2021.
  * [191] R. B. Morgan.  A restarted GMRES method augmented with eigenvectors.  SIAM J. Matrix Anal. Appl., 16(4):1154–1171, 1995.
  * [192] R. B. Morgan.  Implicitly restarted GMRES and Arnoldi methods for nonsymmetric systems of equations.  SIAM J. Matrix Anal. Appl., 21(4):1112–1135, 2000.
  * [193] R. B. Morgan.  GMRES with deflated restarting.  SIAM J. Sci. Comput., 24(1):20–37, 2002.
  * [194] R. B. Morgan.  Restarted block-GMRES with deflation of eigenvalues.  Appl. Numer. Math., 54(2):222–236, 2005.
  * [195] R. B. Morgan, Z. Yang, and B. Zhong.  Pseudoeigenvector bases and deflated GMRES for highly nonnormal matrices.  Numer. Linear Algebra Appl., 23(6):1032–1045, 2016.
  * [196] K. Morikuni and M. Rozložník.  On GMRES for singular EP and GP systems.  SIAM J. Matrix Anal. Appl., 39(2):1033–1048, 2018.
  * [197] K. Moriya and T. Nodera.  The DEFLATED-GMRES($m$, $k$) method with switching the restart frequency dynamically.  Numer. Linear Algebra Appl., 7(7-8):569–584, 2000.
  * [198] N. M. Nachtigal, S. C. Reddy, and L. N. Trefethen.  How fast are nonsymmetric matrix iterations?  SIAM J. Matrix Anal. Appl., 13(3):778–795, 1992.
  * [199] N. M. Nachtigal, L. Reichel, and L. N. Trefethen.  A hybrid GMRES algorithm for nonsymmetric linear systems.  SIAM J. Matrix Anal. Appl., 13(3):796–825, 1992.
  * [200] O. Nevanlinna.  Convergence of Iterations for Linear Equations.  Birkhäuser, Basel, 1993.
  * [201] N. Nigro, M. Storti, S. Idelsohn, and T. Tezduyar.  Physics based GMRES preconditioner for compressible and incompressible navier-stokes equations.  Comput. Meth. Appl. Mech. Eng., 154(3):203–228, 1998.
  * [202] Q. Niu, L. Lu, and J. Zhou.  Accelerate weighted GMRES by augmenting error approximations.  Int. J. Comput. Math., 87(9):2101–2112, 2010.
  * [203] T. Ohtsuka.  A continuation/GMRES method for fast computation of nonlinear receding horizon control.  Automatica, 40(4):563–574, 2004.
  * [204] C. C. Paige, M. Rozlozník, and Z. Strakoš.  Modified Gram-Schmidt (MGS), least squares, and backward stability of MGS-GMRES.  SIAM J. Matrix Anal. Appl., 28(1):264–284, 2006.
  * [205] C. C. Paige and M. A. Saunders.  Solution of sparse indefinite systems of linear equations.  SIAM J. Numer. Anal., 12(4):617–629, 1975.
  * [206] C. C. Paige and Z. Strakoš.  Residual and backward error bounds in minimum residual Krylov subspace methods.  SIAM J. Sci. Comput., 23(6):1898–1923, 2002.
  * [207] M. L. Parks, E. de Sturler, G. Mackey, D. D. Johnson, and S. Maiti.  Recycling Krylov subspaces for sequences of linear systems.  SIAM J. Sci. Comput., 28(5):1651–1674, 2006.
  * [208] J. W. Pearson and J. Pestana.  Preconditioners for Krylov subspace methods: An overview.  GAMM-Mitteilungen, 43(4):e202000015, 2020.
  * [209] J. Pestana and A. J. Wathen.  On the choice of preconditioner for minimum residual methods for non-Hermitian matrices.  J. Comput. Appl. Math., 249:57–68, 2013.
  * [210] B. Philippe and L. Reichel.  On the generation of Krylov subspace bases.  Appl. Numer. Math., 62(9):1171–1186, 2012.
  * [211] L. Reichel and Q. Ye.  Breakdown-free GMRES for singular systems.  SIAM J. Matrix Anal. Appl., 26(4):1001–1021, 2005.
  * [212] M. Robbé and M. Sadkane.  A convergence analysis of GMRES and FOM methods for Sylvester equations.  Numer. Algorithms, 30(1):71–89, 2002.
  * [213] M. Robbé and M. Sadkane.  Exact and inexact breakdowns in the block GMRES method.  Linear Algebra Appl., 419(1):265–285, 2006.
  * [214] S. Röllin and W. Fichtner.  Improving the accuracy of GMRes with deflated restarting.  SIAM J. Sci. Comput., 30(1):232–245, 2008.
  * [215] Y. Saad.  Variations on Arnoldi’s method for computing eigenelements of large unsymmetric matrices.  Linear Algebra Appl., 34:269–295, 1980.
  * [216] Y. Saad.  Krylov subspace methods for solving large unsymmetric linear systems.  Math. Comput., 37(155):105–126, 1981.
  * [217] Y. Saad.  A flexible inner-outer preconditioned GMRES algorithm.  SIAM J. Sci. Comput., 14(2):461–469, 1993.
  * [218] Y. Saad.  Analysis of augmented Krylov subspace methods.  SIAM J. Matrix Anal. Appl., 18(2):435–449, 1997.
  * [219] Y. Saad.  Further analysis of minimum residual iterations.  Numer. Linear Algebra Appl., 7(2):67–93, 2000.
  * [220] Y. Saad.  Iterative Methods for Sparse Linear Systems.  SIAM, Philadelphia, PA, 2nd edition, 2003.
  * [221] Y. Saad.  Numerical Methods for Large Eigenvalue Problems.  SIAM, Philadelphia, PA, 2nd edition, 2011.
  * [222] Y. Saad.  Iterative methods for linear systems of equations: A brief historical journey.  Contemp. Math., 754:197–215, 2020.
  * [223] Y. Saad and M. H. Schultz.  GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems.  SIAM J. Sci. Stat. Comput., 7(3):856–869, 1986.
  * [224] Y. Saad and K. Wu.  DQGMRES: a direct quasi-minimal residual algorithm based on incomplete orthogonalization.  Numer. Linear Algebra Appl., 3(4):329–343, 1996.
  * [225] H. Saberi Najafi and H. Zareamoghaddam.  A new computational GMRES method.  Appl. Math. Comput., 199(2):527–534, 2008.
  * [226] G. Sacchi and V. Simoncini.  A GMRES convergence analysis for localized invariant subspace ill-conditioning.  SIAM J. Matrix Anal. Appl., 40(2):542–563, 2019.
  * [227] H. Sadok.  CMRH: A new method for solving nonsymmetric linear systems based on the Hessenberg reduction algorithm.  Numer. Algorithms, 20(4):303–321, 1999.
  * [228] H. Sadok.  Analysis of the convergence of the minimal and the orthogonal residual methods.  Numer. Algorithms, 40(2):201–216, 2005.
  * [229] H. Sadok and D. B. Szyld.  A new look at CMRH and its relation to GMRES.  BIT, 52(2):485–501, 2012.
  * [230] A. K. Saibaba, T. Bakhos, and P. K. Kitanidis.  A flexible Krylov solver for shifted systems with application to oscillatory hydraulic tomography.  SIAM J. Sci. Comput., 35(6):A3001–A3023, 2013.
  * [231] P. Sanan, S. M. Schnepp, and D. A. May.  Pipelined, flexible Krylov subspace methods.  SIAM J. Sci. Comput., 38(5):C441–C470, 2016.
  * [232] J. N. Shadid and R. S. Tuminaro.  Sparse iterative algorithm software for large-scale MIMD machines: An initial discussion and implementation.  Concurr. Pract. Exp., 4(6):481–497, 1992.
  * [233] A. Sidi.  Extrapolation vs. projection methods for linear systems of equations.  J. Comput. Appl. Math., 22(1):71–88, 1988.
  * [234] A. Sidi.  DGMRES: A GMRES-type algorithm for Drazin-inverse solution of singular non-symmetric linear systems.  Linear Algebra Appl., 335(1):189–204, 2001.
  * [235] R. B. Sidje and N. Winkles.  Evaluation of the performance of inexact GMRES.  J. Comput. Appl. Math., 235(8):1956–1975, 2011.
  * [236] J. A. Sifuentes, M. Embree, and R. B. Morgan.  GMRES convergence for perturbed coefficient matrices, with application to approximate deflation preconditioning.  SIAM J. Matrix Anal. Appl., 34(3):1066–1088, 2013.
  * [237] V. Simoncini.  On a non-stagnation condition for GMRES and application to saddle point matrices.  Electron. Trans. Numer. Anal., 37:202–213, 2010.
  * [238] V. Simoncini.  Computational methods for linear matrix equations.  SIAM Rev., 58(3):377–441, 2016.
  * [239] V. Simoncini and E. Gallopoulos.  An iterative method for nonsymmetric systems with multiple right-hand sides.  SIAM J. Sci. Comput., 16(4):917–933, 1995.
  * [240] V. Simoncini and E. Gallopoulos.  Convergence properties of block GMRES and matrix polynomials.  Linear Algebra Appl., 247:97–119, 1996.
  * [241] V. Simoncini and D. B. Szyld.  Theory of inexact Krylov subspace methods and applications to scientific computing.  SIAM J. Sci. Comput., 25(2):454–477, 2003.
  * [242] V. Simoncini and D. B. Szyld.  On the occurrence of superlinear convergence of exact and inexact Krylov subspace methods.  SIAM Rev., 47(2):247–272, 2005.
  * [243] V. Simoncini and D. B. Szyld.  Recent computational developments in Krylov subspace methods for linear systems.  Numer. Linear Algebra Appl., 14(1):1–59, 2007.
  * [244] V. Simoncini and D. B. Szyld.  New conditions for non-stagnation of minimal residual methods.  Numer. Math., 109(3):477–487, 2008.
  * [245] D. A. Smith, W. F. Ford, and A. Sidi.  Extrapolation methods for vector sequences.  SIAM Rev., 29(2):199–233, 1987.
  * [246] A. Smoktunowicz, J. L. Barlow, and J. Langou.  A note on the error analysis of classical Gram-Schmidt.  Numer. Math., 105(2):299–313, 2006.
  * [247] K. M. Soodhalter.  Block Krylov subspace recycling for shifted systems with unrelated right-hand sides.  SIAM J. Sci. Comput., 38(1):A302–A324, 2016.
  * [248] K. M. Soodhalter.  Two recursive GMRES-type methods for shifted linear systems with general preconditioning.  Electron. Trans. Numer. Anal., 45:499–523, 2016.
  * [249] K. M. Soodhalter, E. de Sturler, and M. E. Kilmer.  A survey of subspace recycling iterative methods.  GAMM-Mitteilungen, 43(4):e202000016, 2020.
  * [250] K. M. Soodhalter, D. B. Szyld, and F. Xue.  Krylov subspace recycling for sequences of shifted linear systems.  Appl. Numer. Math., 81:105–118, 2014.
  * [251] M. Sosonkina, L. T. Watson, R. K. Kapania, and H. F. Walker.  A new adaptive GMRES algorithm for achieving high accuracy.  Numer. Linear Algebra Appl., 5(4):275–297, 1998.
  * [252] B. S. Southworth, A. A. Sivas, and S. Rhebergen.  On fixed-point, Krylov, and $2\times 2$ block preconditioners for nonsymmetric problems.  SIAM J. Matrix Anal. Appl., 41(2):871–900, 2020.
  * [253] G. Starke.  Field-of-values analysis of preconditioned iterative methods for nonsymmetric elliptic problems.  Numer. Math., 78(1):103–117, 1997.
  * [254] G. Starke and R. S. Varga.  A hybrid Arnoldi-Faber iterative method for nonsymmetric systems of linear equations.  Numer. Math., 64(1):213–240, 1993.
  * [255] D.-L. Sun, T.-Z. Huang, Y.-F. Jing, and B. Carpentieri.  A block GMRES method with deflated restarting for solving linear systems with multiple shifts and multiple right-hand sides.  Numer. Linear Algebra Appl., 25(5):e2148, 2018.
  * [256] K. Swirydowicz, J. Langou, S. Ananthan, U. Yang, and S. Thomas.  Low synchronization Gram-Schmidt and generalized minimal residual algorithms.  Numer. Linear Algebra Appl., 28(2):e2343, 2021.
  * [257] A. Tajaddini, G. Wu, F. Saberi-Movahed, and N. Azizizadeh.  Two new variants of the simpler block GMRES method with vector deflation and eigenvalue deflation for multiple linear systems.  J. Sci. Comput., 86(1):9:1–9:33, 2021.
  * [258] P. Tichý, J. Liesen, and V. Faber.  On worst-case GMRES, ideal GMRES, and the polynomial numerical hull of a Jordan block.  Electron. Trans. Numer. Anal., 26:453–473, 2007.
  * [259] D. Titley-Peloquin, J. Pestana, and A. J. Wathen.  GMRES convergence bounds that depend on the right-hand-side vector.  IMA J. Numer. Anal., 34(2):462–479, 2014.
  * [260] K.-C. Toh.  GMRES vs. ideal GMRES.  SIAM J. Matrix Anal. Appl., 18(1):30–36, 1997.
  * [261] L. N. Trefethen.  Approximation Theory and Numerical Linear Algebra, pages 336–360.  Chapman and Hall, London, 1990.
  * [262] L. N. Trefethen and M. Embree.  Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators.  Princeton University Press, Princeton, 2005.
  * [263] K. Turner and H. F. Walker.  Efficient high accuracy solutions with GMRES($m$).  SIAM J. Sci. Stat. Comput., 13(3):815–825, 1992.
  * [264] J. van den Eshof and G. L. G. Sleijpen.  Inexact Krylov subspace methods for linear systems.  SIAM J. Matrix Anal. Appl., 26(1):125–153, 2004.
  * [265] A. van der Sluis and H. A. van der Vorst.  The rate of convergence of conjugate gradients.  Numer. Math., 48(5):543–560, 1986.
  * [266] H. A. van der Vorst.  Iterative Krylov Methods for Large Linear Systems.  Cambridge University Press, Cambridge, 2003.
  * [267] H. A. van der Vorst and C. Vuik.  The superlinear convergence behaviour of GMRES.  J. Comput. Appl. Math., 48(3):327–341, 1993.
  * [268] H. A. van der Vorst and C. Vuik.  GMRESR: A family of nested GMRES methods.  Numer. Linear Algebra Appl., 1(4):369–386, 1994.
  * [269] M. B. van Gijzen.  A polynomial preconditioner for the GMRES algorithm.  J. Comput. Appl. Math., 59(1):91–107, 1995.
  * [270] E. Vecharynski and J. Langou.  The cycle-convergence of restarted GMRES for normal matrices is sublinear.  SIAM J. Sci. Comput., 32(1):186–196, 2010.
  * [271] E. Vecharynski and J. Langou.  Any admissible cycle-convergence behavior is possible for restarted GMRES at its initial cycles.  Numer. Linear Algebra Appl., 18(3):499–511, 2011.
  * [272] P. K. W. Vinsome.  Orthomin, an iterative method for solving sparse sets of simultaneous linear equations.  In Proceedings of the Fourth Symposium on Reservoir Simulation, pages 149–159. Society of Petroleum Engineers of AIME, Los Angeles, CA, USA, 1976\.
  * [273] B. Vital.  Etude de Quelques Méthodes de Résolution de Problèmes Linéaires de Grande Taille sur Multiprocesseur.  PhD thesis, Université de Rennes 1, 1990.
  * [274] C. Vuik.  New insights in GMRES-like methods with variable preconditioners.  J. Comput. Appl. Math., 61(2):189–204, 1995.
  * [275] C. Vuik and H. A. van der Vorst.  A comparison of some GMRES-like methods.  Linear Algebra Appl., 160:131–162, 1992.
  * [276] D. N. Wakam and J. Erhel.  Parallelism and robustness in GMRES with a newton basis and deflated restarting.  Electron. Trans. Numer. Anal., 40:381–406, 2013.
  * [277] H. F. Walker.  Implementation of the GMRES method using Householder transformations.  SIAM J. Sci. Stat. Comput., 9(1):152–163, 1988.
  * [278] H. F. Walker.  Implementations of the GMRES method.  Comput. Phys. Commun., 53(1-3):311–320, 1989.
  * [279] H. F. Walker and P. Ni.  Anderson acceleration for fixed-point iterations.  SIAM J. Numer. Anal., 49(4):1715–1735, 2011.
  * [280] H. F. Walker and L. Zhou.  A simpler GMRES.  Numer. Linear Algebra Appl., 1(6):571–581, 1994.
  * [281] A. J. Wathen.  Preconditioning.  Acta Numer., 24:329–376, 2015.
  * [282] L. Wigton, N. Yu, and D. Young.  GMRES acceleration of computational fluid dynamics codes.  In Proceedings of the 7th Computational Physics Conference, pages 67–74, Cincinnati, OH, USA, 1985. AIAA.
  * [283] K. Wu and H. Simon.  Thick-restart Lanczos method for large symmetric eigenvalue problems.  SIAM J. Matrix Anal. Appl., 22(2):602–616, 2000.
  * [284] I. Yamazaki, H. Anzt, S. Tomov, M. Hoemmen, and J. Dongarra.  Improving the performance of CA-GMRES on multicores with multiple GPUs.  In Proceedings of the International Parallel and Distributed Processing Symposium, pages 382–391, Phoenix, AZ, USA, 2014. IEEE.
  * [285] I. Yamazaki, M. Hoemmen, P. Luszczek, and J. Dongarra.  Improving performance of GMRES by reducing communication and pipelining global collectives.  In Proceedings of the International Parallel and Distributed Processing Symposium, pages 1118–1127, Lake Buena Vista, FL, USA, 2017. IEEE.
  * [286] I. Yamazaki, S. Thomas, M. Hoemmen, E. G. Boman, K. Świrydowicz, and J. J. Elliott.  Low-synchronization orthogonalization schemes for $s$-step and pipelined Krylov solvers in Trilinos.  In Proceedings of the SIAM Conference on Parallel Processing for Scientific Computing, pages 118–128, Seattle, WA, USA, 2020. SIAM.
  * [287] I. Yamazaki, S. Tomov, and J. Dongarra.  Deflation strategies to improve the convergence of communication-avoiding GMRES.  In Proceedings of the 5th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems, pages 39–46, New Orleans, LA, USA, 2014. IEEE.
  * [288] I. Yamazaki, S. Tomov, J. Kurzak, J. Dongarra, and J. L. Barlow.  Mixed-precision block Gram Schmidt orthogonalization.  In Proceedings of the 6th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems, pages 2:1–2:8, New York, NY, USA, 2015. ACM.
  * [289] L. M. Yang, A. Fox, and G. Sanders.  Rounding error analysis of mixed precision block Householder QR algorithms.  SIAM J. Sci. Comput., 43(3):A1723–A1753, 2021.
  * [290] X. Ye, Y. Xi, and Y. Saad.  Proxy-GMRES: Preconditioning via GMRES in polynomial space.  SIAM J. Matrix Anal. Appl., 42(3):1248–1267, 2021.
  * [291] D. M. Young and K. C. Jea.  Generalized conjugate-gradient acceleration of nonsymmetrizable iterative methods.  Linear Algebra Appl., 34:159–194, 1980.
  * [292] N. A. Zadeh, A. Tajaddini, and G. Wu.  Weighted and deflated global GMRES algorithms for solving large Sylvester matrix equations.  Numer. Algorithms, 82(1):155–181, 2019.
  * [293] I. Zavorin, D. P. O’Leary, and H. Elman.  Complete stagnation of GMRES.  Linear Algebra Appl., 367:165–183, 2003.
  * [294] J. Zítko.  Generalization of convergence conditions for a restarted GMRES.  Numer. Linear Algebra Appl., 7(3):117–131, 2000.
  * [295] J. Zítko.  Some remarks on the restarted and augmented GMRES method.  Electron. Trans. Numer. Anal., 31:221–227, 2008.