# Multistage Mixed Precision Iterative Refinement111We acknowledge funding
from the Charles University PRIMUS project No. PRIMUS/19/SCI/11, Charles
University Research Program No. UNCE/SCI/023, and the Exascale Computing
Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy
Office of Science and the National Nuclear Security Administration.

Eda Oktay and Erin Carson

##### Abstract.

Low precision arithmetic, in particular half precision (16-bit) floating point
arithmetic, is now available in commercial hardware. It is anticipated that
exascale machines, capable of $10^{18}$ floating point operations per second,
will expose a range of hardware precisions through GPUs and other
accelerators. Using lower precision can offer significant savings in
computation and communication costs with proportional savings in energy.
Motivated by this, there have recently emerged a number of new iterative
refinement schemes for solving linear systems $Ax=b$, both based on standard
LU factorization and GMRES solvers, that exploit multiple different
precisions. Each particular algorithm and each combination of precisions leads
to different condition number-based constraints for convergence of the
backward and forward errors, and each has different performance costs. Given
that the user may not necessarily know the condition number of their matrix a
priori, it may be difficult to select the optimal variant for their problem.

In this work, we develop a three-stage mixed precision iterative refinement
solver which aims to combine existing mixed precision approaches to balance
performance and accuracy and improve usability. For a given combination of
precisions, the algorithm begins with the least expensive approach and
convergence is monitored via inexpensive computations with quantities produced
during the iteration. If slow convergence or divergence is detected using
particular stopping criteria, the algorithm switches to use more expensive,
but more reliable GMRES-based refinement approaches. After presenting the
algorithm and its details, we perform extensive numerical experiments on a
variety of random dense problems and problems from real applications. Our
experiments demonstrate that the theoretical constraints derived in the
literature are often overly strict in practice, further motivating the need
for a multistage approach.

##  1 Introduction

Iterative refinement (IR), shown in Algorithm 1, is frequently used in solving
linear systems $Ax=b$ to improve the accuracy of the computed approximate
solution $\hat{x}$. Typically, one computes an initial approximate solution
$\hat{x}_{0}$ using Gaussian elimination with partial pivoting (GEPP), saving
the approximate factorization $A\approx\hat{L}\hat{U}$. After computing the
residual $\hat{r}=b-A\hat{x}$ (potentially in higher precision), one reuses
$\hat{L}$ and $\hat{U}$ to solve the system $A\hat{d}=\hat{r}$. The original
approximate solution is subsequently refined by adding the corrective term,
$\hat{x}=\hat{x}+\hat{d}$. This process can be repeated until either the
desired accuracy is reached or the iterative refinement process converges to
an approximate solution. Iterative refinement algorithms that exploit multiple
different precisions have gained much attention recently due to the emergence
of mixed precision hardware, and form the basis for the recently developed
HPL-AI benchmark, on which today’s top supercomputers exceed exascale
performance; see e.g., [12, 15, 19].

Algorithm 1 General Iterative Refinement Scheme

1:$n\times n$ matrix $A$; right-hand side $b$; maximum number of refinement
steps $i_{max}$.

2:Approximate solution $x_{i+1}$ to $Ax=b$.

3:Compute LU factorization $A=LU$ in precision $u_{f}$.

4:Solve $Ax_{0}=b$ by substitution in precision $u_{f}$; store $x_{0}$ in
precision $u$.

5:for $i=0$: $i_{max}-1$ do

6: Compute $r_{i}=b-Ax_{i}$ in precision $u_{r}$; store in precision $u$.

7: Solve $Ad_{i+1}=r_{i}$ in precision $u_{s}$; store $d_{i+1}$ in precision
$u$.

8: Compute $x_{i+1}=x_{i}+d_{i+1}$ in precision $u$.

9: if converged then return $x_{i+1}$. end if

Algorithm 1 shows a general mixed precision iterative refinement scheme.
Following the work of [4], there are four different precisions specified:
$u_{f}$ denotes the factorization precision, $u$ denotes the working
precision, $u_{r}$ denotes the precision for the residual computation, and
$u_{s}$ denotes the effective precision with which the correction equation is
solved. This final precision, $u_{s}$, is not a hardware precision but will
rather depend on the particular solver used in line 7 and the particular
precision(s) used within the solver (which may be different from $u$, $u_{f}$,
and $u_{r}$).

Given the choice of precisions $u$, $u_{f}$, and $u_{r}$, and the choice of a
solver to use in line 7, there are many potential variants of iterative
refinement. We refer to any variant that solves the correction equation in
line 7 using triangular solves with the computed LU factors as “standard
iterative refinement” (SIR). For SIR, the effective precision of the
correction solve is limited by the precision with which the factorization is
computed, and thus we have $u_{s}=u_{f}$. The most commonly-used SIR variant
is what we call “traditional” iterative refinement, in which $u_{f}=u$ and
$u_{r}=u^{2}$. Traditional iterative refinement was used already by Wilkinson
in 1948 [20], and was analyzed in fixed point arithmetic by Wilkinson in 1963
[21] and in floating point arithmetic by Moler in 1967 [17]. There have also
been analyses of fixed precision iterative refinement, in which
$u_{f}=u=u_{r}$ [14], as well as low precision factorization iterative
refinement, in which $u_{f}^{2}=u=u_{r}$ [16]. Motivated by the trend of low
and mixed precision capabilities in hardware, the authors in [4] developed and
analyzed a three-precision iterative refinement scheme, in which $u_{f}$, $u$,
and $u_{r}$ may differ, which generalizes (and in some cases, improves the
bounds for) many existing variants. For references to analyses of variants of
iterative refinement, see [4, Table 1.1].

If $A$ is very ill conditioned or badly scaled, SIR can fail, i.e., the error
is not eventually bounded by a small multiple of machine precision. In extreme
cases of ill conditioning, the error can grow with each refinement step. In
[3] and [4], the authors developed a mixed precision GMRES-based iterative
refinement scheme (GMRES-IR), shown in Algorithm 2. The only difference with
SIR is the way in which the correction equation is solved. The idea is that
instead of using the LU factors to solve for the correction in each step, one
can instead use these factors as (left) preconditioners for a preconditioned
GMRES method which solves for the correction. In this way, the effective solve
precision becomes $u_{s}=u$, and more ill-conditioned problems can be handled
relative to SIR. The GMRES-based refinement approaches have seen much success
in practice. Current GMRES-based iterative refinement variants are implemented
in the MAGMA library (2.5.0) and in the NVIDIA cuSOLVER library. GMRES-IR also
forms the basis for the new HPL-AI benchmark, used to rank computers in the
TOP500 list [12]. Experiments detailing the performance benefits of three-
precision GMRES-IR and SIR can be found in [9].

Algorithm 2 GMRES-IR [3]

1:$n\times n$ matrix $A$; right-hand side $b$; maximum number of refinement
steps $i_{max}$; GMRES convergence tolerance $\tau$.

2:Approximate solution $\hat{x}$ to $Ax=b$.

3:Compute LU factorization $A=LU$ in precision $u_{f}$.

4:Solve $Ax_{0}=b$ by substitution in precision $u_{f}$; store $x_{0}$ in
precision $u$.

5:for $i=0$: $i_{max}-1$ do

6: Compute $r_{i}=b-Ax_{i}$ in precision $u_{r}$; store in precision $u$.

7: Solve $U^{-1}L^{-1}Ad_{i+1}=U^{-1}L^{-1}r_{i}$ by GMRES in working
precision $u$, with matrix-vector products with $\tilde{A}=U^{-1}L^{-1}A$
computed at precision $u^{2}$; store $d_{i+1}$ in precision $u$.

8: Compute $x_{i+1}=x_{i}+d_{i+1}$ in precision $u$.

9: if converged then return $x_{i+1}$. end if

Although GMRES-IR can succeed where SIR fails, each step of GMRES-IR is
potentially more expensive than each step of SIR (albeit still potentially
less expensive than running the entire process in double the working
precision). Whereas each SIR refinement step involves only two triangular
solves in precision $u$, each GMRES-IR step involves two triangular solves in
precision $u^{2}$ in _each_ GMRES iteration (in addition to a matrix-vector
product in precision $u^{2}$ as well as multiple vector operations in
precision $u$). The convergence behavior of each GMRES solve thus plays a
large role, and it is important for the performance of GMRES-IR that each call
to GMRES converges relatively quickly. It is additionally a performance
concern that each GMRES iteration requires computations in higher precision;
in order to obtain the needed bound on backward error of the GMRES solve, the
authors in [3, 4] required that the preconditioned system $U^{-1}L^{-1}A$ (not
formed explicitly) is applied to a vector in precision $u^{2}$.

The authors of [3, 4] suggested that this required use of extra precision was
likely to be overly strict in many practical scenarios. Indeed, most practical
implementations of GMRES-based iterative refinement do not use this extra
precision [2]. This motivated Amestoy et al. [2] to develop an extension of
the GMRES-IR algorithm, called GMRES-IR5. The authors in [2] revisit the proof
of backward stability for GMRES from [18] and develop a bound on the backward
error for a mixed precision GMRES algorithm, where $u_{g}$ is the working
precision and $u_{p}$ is the precision in which the preconditioned matrix is
applied to a vector. This enables a five-precision GMRES-IR scheme, where
$u,u_{f},u_{r},u_{g}$, and $u_{p}$ may take on different values. A particular
improvement over the GMRES-IR scheme of [3, 4] is that the analysis provides
bounds on forward and backward errors for a variant of GMRES-IR in which the
entire GMRES iteration is carried out in a uniform precision, i.e., the
analysis does not require that extra precision is used in applying the
preconditioned matrix to a vector. This particular variant of the algorithm is
thus cheaper than the former GMRES-IR in terms of both time and memory. The
drawback is that it is only theoretically applicable to a smaller set of
problems due to a tighter constraint on condition number. We call the
particular instance of GMRES-IR5 where $u=u_{g}=u_{p}$ “SGMRES-IR” (for
“simpler”), shown in Algorithm 3.

Algorithm 3 SGMRES-IR (a particular variant of GMRES-IR5 [2])

1:$n\times n$ matrix $A$; right-hand side $b$; maximum number of refinement
steps $i_{max}$; GMRES convergence tolerance $\tau$.

2:Approximate solution $x_{i+1}$ to $Ax=b$.

3:Compute LU factorization $A=LU$ in precision $u_{f}$.

4:Solve $Ax_{0}=b$ by substitution in precision $u_{f}$; store $x_{0}$ in
precision $u$.

5:for $i=0$: $i_{max}-1$ do

6: Compute $r_{i}=b-Ax_{i}$ in precision $u_{r}$; store in precision $u$.

7: Solve $U^{-1}L^{-1}Ad_{i+1}=U^{-1}L^{-1}r_{i}$ by GMRES in working
precision $u$, with matrix-vector products with $\tilde{A}=U^{-1}L^{-1}A$
computed at precision $u$; store $d_{i+1}$ in precision $u$.

8: Compute $x_{i+1}=x_{i}+d_{i+1}$ in precision $u$.

9: if converged then return $x_{i+1}$. end if

In terms of both cost and range of condition numbers to which the algorithm
can be applied, we expect SGMRES-IR to be, in general, somewhere between SIR
and GMRES-IR. For example, if we use single precision for $u_{f}$, double
precision for $u$, and quadruple precision for $u_{r}$, SIR is guaranteed to
converge to the level of double precision in both forward and backward errors
as long as the infinity-norm condition number of the matrix $A$,
$\kappa_{\infty}(A)=\|A\|_{\infty}\|A^{-1}\|_{\infty}$, is less than $2\cdot
10^{7}$. For SGMRES-IR, this constraint becomes $\kappa_{\infty}(A)\leq
10^{10}$. GMRES-IR, on the other hand, only requires $\kappa_{\infty}(A)\leq
2\cdot 10^{15}$. Thus going from SIR to SGMRES-IR to GMRES-IR, we expect that
these algorithms will be increasingly expensive but also expect that they will
converge for increasingly ill-conditioned matrices. We note this may not
always be the case. For some precision combinations, SGMRES-IR has a tighter
constraint on condition number than SIR (see [2]), although this may be an
artifact of the analysis. Also, the metric of relative “cost” is difficult to
determine a priori, in particular between SGMRES-IR and GMRES-IR, since it
depends on the number of GMRES iterations required in each refinement step.

It is thus difficult to choose a priori which particular variant of iterative
refinement is the most appropriate for a particular problem. It is often the
case in practice that a user does not know the condition number of their
matrix. Further, even if the matrix condition number meets the constraint for
convergence for the chosen iterative refinement algorithm, the convergence
rate may be unacceptably slow, or each step may require so many GMRES
iterations that it becomes impractical.

In this work, we aim to solve this problem through the development of a three-
stage, three-precision iterative refinement scheme, which we abbreviate TSIR.
Our approach automatically switches between variants if slow convergence is
detected using stopping criteria based on the work in [7]. In Section 2, we
first present the TSIR algorithm and give a motivating numerical example. We
then give details of the stopping criteria used and summarize the analysis for
each algorithm variant from [4] and [2]. In Section 3 we present more thorough
numerical experiments on both random dense matrices and matrices from the
SuiteSparse collection [6]. We conclude and discuss future extensions in
Section 4.

##  2 The TSIR algorithm

In order to balance reliability and cost, we develop a three-stage iterative
refinement algorithm (TSIR), presented in Algorithm 4. The algorithm starts
with three-precision SIR (as in [4]) and, using the stopping criteria
developed in [7], switches to SGMRES-IR if the algorithm is not converging at
an acceptable rate (or not converging at all). Then, using the same stopping
criteria along with an additional constraint on the number of GMRES iterations
per refinement step, the algorithm may choose to switch a second time to the
GMRES-IR algorithm (which uses higher precision in applying the preconditioned
matrix to a vector within GMRES). Switching between different variants in this
way, TSIR is potentially less expensive than GMRES-IR for well-conditioned
matrices and more reliable than SIR for ill-conditioned matrices.

Algorithm 4 Three-Stage Iterative Refinement (TSIR)

1:$n\times n$ matrix $A$; right-hand-side $b$; maximum number of refinement
steps of each type $i_{max}$; GMRES convergence tolerance $\tau$; stopping
criteria parameter $\rho_{thresh}$; maximum GMRES iterations $k_{max}$.

2:Approximate solution $x_{i+1}$ to $Ax=b$, boolean cged.

3:Compute LU factorization $A=LU$ in precision $u_{f}$.

4:Solve $Ax_{0}=b$ by substitution in precision $u_{f}$; store $x_{0}$ in
precision $u$.

5:Initialize: $d_{0}=\infty$; alg $=$ SIR; iter $=0$; $i=0$; converged $=0$;
$\rho_{max}=0$.

6:while not cged do

7: Compute $r_{i}=b-Ax_{i}$ in precision $u_{r}$; Scale
$r_{i}=r_{i}/\|r_{i}\|_{\infty}$; store in precision $u$.

8: if alg $=$ SIR then

9: iter $=$ iter $+1$

10: Compute $d_{i+1}=U^{-1}(L^{-1}r_{i})$ in precision $u_{f}$; store
$d_{i+1}$ in precision $u$.

11: if $d_{i+1}$ contains Inf or NaN then alg $=$ SGMRES-IR; iter $=0$; break.
end if

12: Compute $x_{i+1}=x_{i}+\|r_{i}\|_{\infty}d_{i+1}$ at precision $u$.

13: $z=\|d_{i+1}\|_{\infty}/\|x_{i}\|_{\infty}$;
$v=\|d_{i+1}\|_{\infty}/\|d_{i}\|_{\infty}$; $\rho_{max}=\max(\rho_{max},v)$;
$\phi_{i}=z/(1-\rho_{max})$

14: if $z\leq u$ or $v\geq\rho_{thresh}$ or iter $>i_{max}$ or
$\phi_{i}\leq\sqrt{n}u$ then

15: if not converged then

16: alg $=$ SGMRES-IR; iter $=0$.

17: if $\phi_{i}>\phi_{0}$ then $x_{i+1}=x_{0}$ end if

18: else

19: cged $=1$

20: else if alg $=$ SGMRES-IR then

21: iter $=$ iter $+1$

22: Solve $U^{-1}L^{-1}Ad_{i+1}=U^{-1}L^{-1}r_{i}$ by GMRES in precision $u$
with matrix-vector products with $\tilde{A}=U^{-1}L^{-1}A$ computed at
precision $u$; store $d_{i+1}$ in precision $u$.

23: Compute $x_{i+1}=x_{i}+\|r_{i}\|_{\infty}d_{i+1}$ at precision $u$.

24: $z=\|d_{i+1}\|_{\infty}/\|x_{i}\|_{\infty}$;
$v=\|d_{i+1}\|_{\infty}/\|d_{i}\|_{\infty}$; $\rho_{max}=\max(\rho_{max},v)$;
$\phi_{i}=z/(1-\rho_{max})$

25: if $z\leq u$ or $v\geq\rho_{thresh}$ or $\texttt{iter}>i_{max}$ or
$k_{GMRES}>k_{max}$ or $\phi_{i}\leq\sqrt{n}u$ then

26: if not converged then

27: alg $=$ GMRES-IR; iter $=0$.

28: if $\phi_{i}>\phi_{0}$ then $x_{i+1}=x_{0}$ end if

29: else

30: cged $=1$

31: else if alg $=$ GMRES-IR then

32: iter $=$ iter $+1$

33: Solve $U^{-1}L^{-1}Ad_{i+1}=U^{-1}L^{-1}r_{i}$ by GMRES in precision $u$
with matrix-vector products with $\tilde{A}=U^{-1}L^{-1}A$ computed at
precision $u^{2}$; store $d_{i+1}$ in precision $u$.

34: Compute $x_{i+1}=x_{i}+\|r_{i}\|_{\infty}d_{i+1}$ at precision $u$.

35: $z=\|d_{i+1}\|_{\infty}/\|x_{i}\|_{\infty}$;
$v=\|d_{i+1}\|_{\infty}/\|d_{i}\|_{\infty}$; $\rho_{max}=\max(\rho_{max},v)$;
$\phi_{i}=z/(1-\rho_{max})$

36: if iter $>i_{max}$ or $\phi_{i}\leq\sqrt{n}u$ then

37: if converged then cged $=1$ end if

38: return

39: $i=i+1$

Before explaining the details of the algorithm, we begin with a brief
motivating example illustrating how TSIR works. We use precisions
$(u_{f},u,u_{r})$ = (single, double, quad) and generate random dense matrices
using the MATLAB command `gallery(’randsvd’,n,kappa(i),2)`, where $n=100$ and
`kappa` is the array of the desired 2-norm condition numbers. Here we test
2-norm condition numbers $10^{1},10^{5}$, and $10^{16}$. We set $b$ to be a
vector of normally distributed numbers generated by the MATLAB command
`randn`. For reproducibility, we use the MATLAB command `rng(1)` to seed the
random number generator before generating each linear system. The GMRES
convergence tolerance $\tau$, which also appears in Algorithms 2 and 3,
dictates the stopping criterion for the inner GMRES iterations. The algorithm
is considered to be converged if the relative (preconditioned) residual norm
drops below $\tau$. In all tests here we use $\tau=10^{-6}$.

Figures 1, 2, and 3 show results for condition numbers $10^{1},10^{5}$, and
$10^{16}$, respectively. The red, blue, and green lines in the figures show
the behavior of the forward error `ferr` (red), normwise relative backward
error `nbe` (blue), and componentwise relative backward error `cbe` (green).
The dotted black line shows the value of the working precision $u$. If there
is a switch in TSIR, the step at which it happened is marked with a star. If
TSIR uses both SGMRES-IR and GMRES-IR, then there are two stars: the first one
is for SGMRES-IR, the latter one is for GMRES-IR.

The number of refinement steps performed by SIR, SGMRES-IR, GMRES-IR, and TSIR
for these matrices are presented in Table 1. For SIR, the number in each row
gives the number of refinement steps. For SGMRES-IR and GMRES-IR, the number
of parenthetical elements gives the total number of refinement steps, and
element $i$ in the list gives the number of GMRES iterations performed in
refinement step $i$. For example, $(3,4)$ indicates that 2 refinement steps
were performed, the first of which took 3 GMRES iterations and the second of
which took 4. For the TSIR column, the data for SIR, SGMRES-IR, and GMRES-IR
is comma-separated. For example, $2,(2)$ indicates that 2 SIR steps were
performed, the algorithm then switched to SGMRES-IR, and performed 1 SGMRES-IR
step which required 2 GMRES iterations. Since there is no second set of
parentheses, this indicates that there was no switch to GMRES-IR. In all
columns, a dash denotes that the algorithm failed to converge using the
maximum number of refinement steps $i_{max}=10$.

Table 1: Number of SIR, SGMRES-IR, GMRES-IR, and TSIR steps with the number of GMRES iterations for each SGMRES-IR and GMRES-IR step for $(u_{f},u,u_{r})$ = (single,double,quad). $\kappa_{\infty}(A)$ | $\kappa_{2}(A)$ | SIR | SGMRES-IR | GMRES-IR | TSIR
---|---|---|---|---|---
$2\cdot 10^{2}$ | $10^{1}$ | 2 | (2) | (2) | 2
$2\cdot 10^{6}$ | $10^{5}$ | 5 | (2,3) | (2) | 2, (2)
$2\cdot 10^{17}$ | $10^{16}$ | - | - | (3,4) | 2, (3,3), (3,4)

Figure 1: Convergence of errors for a $100\times 100$ random dense matrix with
$\kappa_{2}(A)=10^{1}$ using SIR (top left), SGMRES-IR (top right), GMRES-IR
(bottom left), and TSIR (bottom right), with $(u_{f},u,u_{r})$ = (single,
double, quad). Note that for SGMRES-IR and GMRES-IR, the forward error is
measured as 0 (in double precision) after one refinement step.

Figure 2: Convergence of errors for a $100\times 100$ random dense matrix with
$\kappa_{2}(A)=10^{5}$ using SIR (top left), SGMRES-IR (top right), GMRES-IR
(bottom left), and TSIR (bottom right), with $(u_{f},u,u_{r})$ = (single,
double, quad).

Figure 3: Convergence of errors for a $100\times 100$ random dense matrix with
$\kappa_{2}(A)=10^{16}$ using SIR (top left), SGMRES-IR (top right), GMRES-IR
(bottom left), and TSIR (bottom right), with $(u_{f},u,u_{r})$ = (single,
double, quad).

From Figure 1, we can see that for well-conditioned matrices, SIR quickly
converges to the solution. When the condition number increases closer to
$u_{f}^{-1}$ as in Figure 2, however, SIR convergence begins to slow down, and
TSIR chooses to switch to SGMRES-IR, which then converges in one step. On the
other hand, when the matrix becomes very ill-conditioned as in Figure 3, then
neither SIR nor SGMRES-IR converge in a reasonable number of iterations, and
thus TSIR makes the second switch to GMRES-IR.

This illustrates the benefit of the multistage approach. In the case that the
problem is well conditioned enough that SIR suffices, TSIR will only use SIR.
For the $\kappa_{2}(A)=10^{1}$ case, TSIR behaves exactly the same as SIR, and
is thus likely to be less expensive than SGMRES-IR and GMRES-IR since a single
GMRES iteration is more expensive than an SIR step. In the extremely ill-
conditioned case, both SIR and SGMRES-IR fail to converge. TSIR however does
converge, although at roughly double the cost of GMRES-IR. This is the
inherent tradeoff. Compared to GMRES-IR, we expect TSIR to be less expensive
when the problem is well or reasonably well conditioned and a constant factor
more expensive when problems are extremely ill conditioned relative to the
working precision. Compared to SIR and SGMRES-IR, the benefit is clear: when
SIR and/or SGMRES-IR converges, TSIR will also converge at roughly the same
cost, but TSIR is more reliable and can converge for problems where these
variants don’t.

###  2.1 Algorithm details

We now discuss the TSIR algorithm (Algorithm 4) in more detail. The algorithm
uses $i$ to denote the global number of refinement steps of any type. The iter
variable counts the number of refinement steps of the current refinement
variant; e.g., it will first count the number of SIR steps, and is then reset
to 0 if the algorithm switches to SGMRES-IR. Unlike the single stage
algorithms, here the input parameter $i_{max}$ specifies the maximum number of
refinement steps of each type, meaning we will do up to $i_{max}$ SIR steps
before switching to SGMRES-IR. This means that in theory the maximum number of
steps performed will be $3\cdot i_{max}$.

####  2.1.1 Stopping criteria and convergence detection

In [7], Demmel et al. analyze an iterative refinement scheme in which extra
precision is used in select computations and provide reliable normwise and
componentwise error bounds for the computed solution as well as stopping
criteria. They devise a new approach that adaptively changes the precision
with which the approximate solution is stored based on monitoring the
convergence; if consecutive corrections to the approximate solution are not
decreasing at a sufficient rate, the precision of the approximate solution is
increased. This has the effect of improving the componentwise accuracy.

Following this strategy for monitoring the behavior of iterative refinement
from [7], the TSIR algorithm will switch to the next variant if any of the
following conditions applies:

  1. 1.

$\frac{\|d_{i+1}\|_{\infty}}{\|x_{i}\|_{\infty}}\leq u$ (the correction
$d_{i+1}$ changes solution $x_{i}$ too little),

  2. 2.

$\frac{\|d_{i+1}\|_{\infty}}{\|d_{i}\|_{\infty}}\geq\rho_{thresh}$
(convergence slows down sufficiently),

  3. 3.

$\texttt{iter}\geq i_{max}$ (too many iterations of a particular variant have
been performed), and, for SGMRES-IR,

  4. 4.

$k_{GMRES}\geq k_{max}$ (too many GMRES iterations are performed in one step
of SGMRES-IR)

where $d_{i}$ is the correction of the solution $x_{i}$, $\rho_{thresh}<1$ is
a threshold for convergence, $i_{max}$ is the maximum number of iterations,
$k_{max}$ is the maximum number of GMRES iterations performed per SGMRES-IR
step, and $u$ is the working precision. We reiterate that the fourth condition
above only applies to the switch from SGMRES-IR to GMRES-IR. We note that as
in line 11 in Algorithm 4, in case SIR produces a correction $d_{i+1}$
containing Inf or NaN, we immediately switch to SGMRES-IR without performing
the solution update.

Moreover, the convergence detection in lines 15, 26, and 37 in Algorithm 4 can
be performed using the normwise relative error estimate discussed in [7],

$max\left\\{\dfrac{\frac{\|d_{i+1}\|_{\infty}}{\|x_{i}\|_{\infty}}}{1-\rho_{max}},\gamma
u\right\\}\approx\dfrac{\|x_{i}-x\|_{\infty}}{\|x\|_{\infty}},$ (1)

where $\gamma=max(10,\sqrt{n})$, $n$ is the size of the matrix $A$, and
$\rho_{max}:=max_{j\leq i}\frac{\|d_{j+1}\|_{\infty}}{\|d_{j}\|_{\infty}}$ is
the maximum ratio of successive corrections. The quantity
$(\|d_{i+1}\|_{\infty}/\|x_{i}\|_{\infty})/(1-\rho_{max})$ is stored in
Algorithm 4 as the quantity $\phi_{i}$. If the condition in lines 15, 26, and
37 is triggered, one could then test for convergence by checking whether both
$\phi_{i}\leq\sqrt{n}u$ and $\phi_{i}\geq 0$, and if so, declare convergence
(note that $\phi_{i}<0$ indicates divergence). If convergence is not detected,
we must also decide whether to keep the current approximate solution or reset
the approximate solution to $x_{0}$, since this will be input to the next
stage. We test if $\phi_{i}>\phi_{0}$, and if so, we reset the initial
solution. We note that other measures, such as for the componentwise error,
could also be used to detect convergence.

We note that the criteria used to determine whether the current algorithm
should quit iterating are one iteration behind; in other words, we can only
detect convergence (or non-convergence) in step $i$ after computing the
correction term in step $i+1$. For this reason, the TSIR algorithm will
generally compute at least two steps of a particular variant before deciding
to switch. The two cases where fewer than 2 steps of a variant will occur are
1) when the computed update contains Infs or NaNs (as in line 11 in Algorithm
4), in which case we switch without updating the current solution, and 2)
when, for SGMRES-IR, the number of GMRES iterations exceeds $k_{max}$, since
this is detectable immediately in the current step.

In [7], to determine $\rho_{thresh}$, the authors define ‘cautious’ and
‘aggressive’ settings. Cautious settings produce maximally reliable error
bounds for well-conditioned problems, whereas aggressive ones will lead to
more steps on the hardest problems and usually, but not always, give error
bounds within a factor of 100 of the true error. For cautious mode, the
authors suggest $\rho_{thresh}$ should be set to $0.5$, which was used by
Wilkinson [21], and for aggressive mode, $0.9$. In our experiments we always
use the cautious setting, but we note that $\rho_{thresh}$, $i_{max}$, and
$k_{max}$ should be set according to the relative costs of the different
refinement schemes in practice.

####  2.1.2 Scaling

Because of the smaller range of half precision, simply rounding higher
precision quantities to a lower precision can cause overflow, underflow, or
the introduction of subnormal numbers. In [11], the authors develop a new
algorithm for converting single and double precision quantities to half
precision. This algorithm involves performing a two-sided scaling for
equilibration and then an additional scaling is performed to make full use of
the range of half precision before finally rounding to half precision. In our
algorithm, when half precision is used for the factorization, we first attempt
an LU factorization without scaling. We then test whether the resulting L and
U factors contain Inf or NaN; if so, we retry the LU factorization in line 3
using the two-sided scaling algorithm of [11]. In all cases, regardless of
what precision is used for the factorization, after computing $x_{0}$ in line
4 we test whether the initial solution contains Inf or NaN; if so, we simply
use the zero vector as the initial approximate solution and proceed.

We also incorporate scaling in each refinement step. After computing the
residual $r_{i}$, we scale the result to obtain
$r_{i}=r_{i}/\|r_{i}\|_{\infty}$ (line 7 in Algorithm 4). This scaling is then
undone when we update the approximate solution, via
$x_{i+1}=x_{i}+\|r_{i}\|_{\infty}d_{i+1}$ (lines 12, 23, and 34 in Algorithm
4). As long as $1/\|A\|_{\infty}$ does not underflow and $\|A^{-1}\|_{\infty}$
does not overflow, then this scaling avoids the largest element of $d_{i+1}$
overflowing or underflowing.

###  2.2 Error bounds for different variants

There are various scenarios on the usage of precisions which will yield
different error bounds and different constraints on condition number. Besides
three-precision variants, two precisions or a fixed precision can be used in
the algorithms discussed in this study. We restrict ourselves here to the case
where three different precisions are used, but refer the reader to [4] and [2]
for more general bounds. In particular, we assume that $u_{r}\leq u^{2}$ and
$u\leq u_{f}^{2}$. We also restrict ourselves to IEEE precisions (see Table
2), although we note that alternative formats like bfloat16 [13] could also be
used.

For SIR, both the relative forward and backward errors are guaranteed to
converge to the level of the working precision when
$\kappa_{\infty}(A)<u_{f}^{-1}$. For SGMRES-IR and GMRES-IR, the constraints
for convergence of forward and backward errors differ. Summarizing the
analysis in [2], for our particular SGMRES-IR variant, the constraint on
convergence of the backward error is $\kappa_{\infty}(A)<u^{-1/3}u_{f}^{-1/3}$
and the constraint on the convergence of the forward error is
$\kappa_{\infty}(A)<u^{-1/3}u_{f}^{-2/3}$. It is interesting to note that, as
an artifact of the analysis, the constraint for convergence of the backward
error for this variant is more strict than that for the forward error.
However, since the backward error is bounded by the forward error, the
constraint for both backward and forward error to converge to the level $u$
can be given as $\kappa_{\infty}(A)<u^{-1/3}u_{f}^{-2/3}$.

From [4], for GMRES-IR the constraint for convergence of the backward error to
level $u$ is $\kappa_{\infty}(A)<u^{-1}$. The constraint for the convergence
of the forward error to this level is actually more strict for the GMRES-IR
variant: $\kappa_{\infty}(A)<u^{-1/2}u_{f}^{-1}$. Thus we can take this
tighter constraint for the forward error as our constraint for both forward
and backward errors.

In Table 3 we quantify these constraints on $\kappa_{\infty}(A)$ required for
convergence of the normwise relative backward and forward errors to the level
of the working precision for various precision combinations. To compute the
constraints on $\kappa_{\infty}(A)$, we have used the unit roundoff values in
Table 2.

Table 2: Various IEEE precisions and their units roundoff.  Precision | Unit Roundoff
---|---
fp16 (half) | $4.88\cdot 10^{-4}$
fp32 (single) | $5.96\cdot 10^{-8}$
fp64 (double) | $1.11\cdot 10^{-16}$
fp128 (quad) | $9.63\cdot 10^{-35}$
Table 3: Constraints on $\kappa_{\infty}(A)$ for which the relative forward and normwise backward errors are guarantee to converge to the level $u$ for a given combination of precisions for the different variants of IR.  $u_{f}$ | $u$ | $u_{r}$ | SIR | SGMRES-IR | GMRES-IR
---|---|---|---|---|---
half | single | double | $2\cdot 10^{3}$ | $4\cdot 10^{4}$ | $8\cdot 10^{6}$
half | double | quad | $2\cdot 10^{3}$ | $3\cdot 10^{7}$ | $2\cdot 10^{11}$
single | double | quad | $2\cdot 10^{7}$ | $1\cdot 10^{10}$ | $2\cdot 10^{15}$

##  3 Numerical experiments

In this section we present numerical experiments performed in MATLAB on a
number of synthetic and real-world matrices from SuiteSparse [6] for TSIR and
other single-stage iterative refinement variants in three precisions. For
quadruple precision, we use the Advanpix multiprecision computing toolbox for
MATLAB [1]. To simulate half precision floating point arithmetic, we use the
chop library and associated functions from [10], available at
https://github.com/higham/chop and
https://github.com/SrikaraPranesh/LowPrecision_Simulation. For single and
double precision we use the built-in MATLAB datatypes. All experiments were
performed on a computer with an Intel Core i7-9750H processor with 12 CPUs and
16 GB RAM with OS system Ubuntu 20.04.1 using MATLAB 2020a. Our TSIR algorithm
and associated functions are available through the repository
https://github.com/edoktay/tsir, which includes scripts for generating the
data and plots in this work.

In all experiments, we use $i_{max}=10$ and $\rho_{thresh}=0.5$. For GMRES
convergence tolerance, we use $\tau=10^{-6}$ when the working precision is
single and $\tau=10^{-10}$ when the working precision is double. We set
$k_{max}=0.1n$. For purposes of discerning the actual attainable accuracy of
TSIR, in the experiments we explicitly compute and plot the computed forward
and backward errors and use these as stopping criteria rather than the error
estimate given by (1). For a fair numerical comparison, we also apply the
scalings discussed in Section 2.1.2 in SIR, SGMRES-IR, and GMRES-IR.

For each variant, we compare the number of refinement steps required for
forward and backward errors to reach the level of accuracy corresponding to
the working precision. For GMRES-based approaches and TSIR, the number of
GMRES iterations per step is given parenthetically (see the explanation in
Section 2). A dash (-) in the tables indicates that the forward and/or
backward errors did not converge within $i_{max}$ iterations. Select
convergence plots for TSIR are presented in figures while the number of steps
and iterations for all variants are shown in tables.

Again, the red, blue, and green lines in the figures show the behavior of the
forward error `ferr` (red), normwise relative backward error `nbe` (blue), and
componentwise relative backward error `cbe` (green). The dotted black line
shows the value of the working precision $u$. If there is a switch in TSIR,
the step at which it happened is marked with a star. If TSIR uses both SGMRES-
IR and GMRES-IR, then there are two stars: the first one denoting the switch
from SIR to SGMRES-IR and the second one denoting the switch from SGMRES-IR to
GMRES-IR.

###  3.1 Random dense matrices

We first test our algorithm on random dense matrices generated using the
MATLAB commands `gallery(’randsvd’,n,kappa(i),2)` (in Section 3.1.1) and
`gallery(’randsvd’,n,kappa(i),3)` (in Section 3.1.2), where $n=100$ and
`kappa` is the array of the desired 2-norm condition numbers
$\kappa_{2}(A)=${$10^{1}$, $10^{2}$, $10^{4}$, $10^{5}$, $10^{7}$, $10^{9}$,
$10^{11}$, $10^{14}$}, and `2` and `3` stand for the modes which dictate the
singular value distribution. Mode 2 generates matrices with one singular value
equal to $1/\kappa_{2}(A)$ and the rest of the singular values equal to 1.
Mode 3 generates matrices with geometrically distributed singular values. Mode
3 represents a difficult case for GMRES-based iterative refinement methods. It
is known that a cluster of singular values near zero causes stagnation of
GMRES, and a low precision preconditioner may fail to effectively shift this
cluster away from the origin. We test the algorithms using precision
combinations $(u_{f},u,u_{r})$ = (single, double, quad), $(u_{f},u,u_{r})$ =
(half, single, double), and $(u_{f},u,u_{r})$ = (half, double, quad). For
simplicity, $b$ is chosen to be a vector of normally distributed numbers
generated by the MATLAB command `randn` for all experiments in this section.
For reproducibility, we use the MATLAB command `rng(1)` before generating each
linear system.

####  3.1.1 Random dense matrices with one small singular value

Table 4 shows the convergence behavior of SIR, SGMRES-IR, GMRES-IR, and TSIR
for matrices generated via the MATLAB command
`gallery(’randsvd’,n,kappa(i),2)` using precisions $(u_{f},u,u_{r})$ =
(single, double, quad). From Table 4, we can see that SIR no longer converges
in $i_{max}$ iterations once $\kappa_{\infty}(A)$ exceeds $10^{7}$ (note that
for the test problem with $\kappa_{2}(A)=10^{7}$, SIR does eventually
converge, but requires 19 iterations to do so). As $\kappa(A)$ increases, the
convergence rate of SIR slows down; in cases where 5 SIR steps are required,
notice that the TSIR algorithm deems this too slow and begins switching to
SGMRES-IR after 2 SIR steps. SGMRES-IR converges in just a few refinement
steps with a small numbers of GMRES iterations per step except for the most
ill-conditioned problem. For the most ill-conditioned matrix, SGMRES-IR alone
requires 7 steps to converge. TSIR in this case switches twice: from SIR to
SGMRES-IR, and then finally to GMRES-IR. We plot the convergence trajectories
of TSIR for the $\kappa_{2}(A)=10^{4}$ and $\kappa_{2}(A)=10^{14}$ problems in
Figure 4.

Table 4: Number of SIR, SGMRES-IR, GMRES-IR, and TSIR steps with the number of GMRES iterations for each SGMRES-IR and GMRES-IR step for random dense matrices (mode 2) with various condition numbers $\kappa_{\infty}(A)$, $\kappa_{2}(A)$, using $(u_{f},u,u_{r})$ = (single, double, quad). $\kappa_{\infty}(A)$ | $\kappa_{2}(A)$ | SIR | SGMRES-IR | GMRES-IR | TSIR
---|---|---|---|---|---
$2\cdot 10^{2}$ | $10^{1}$ | 2 | (2) | (2) | 2
$2\cdot 10^{3}$ | $10^{2}$ | 3 | (2) | (2) | 2, (2)
$2\cdot 10^{5}$ | $10^{4}$ | 5 | (2,3) | (2) | 2, (2)
$2\cdot 10^{6}$ | $10^{5}$ | 5 | (2,3) | (2) | 2, (2)
$2\cdot 10^{8}$ | $10^{7}$ | - | (2,3) | (2,3) | 2, (2,3)
$2\cdot 10^{10}$ | $10^{9}$ | - | (2,2) | (2,3) | 2, (2,3)
$2\cdot 10^{12}$ | $10^{11}$ | - | (3,3,4) | (3,4) | 2, (3,3), (3)
$2\cdot 10^{15}$ | $10^{14}$ | - | (3,3,3,4,4,4,4) | (3,4) | 2, (3,3), (3,4)

Figure 4: Convergence of errors in TSIR for random dense matrices (mode 2)
with $\kappa_{2}(A)=10^{4}$ (left), and $\kappa_{2}(A)=10^{14}$ (right) for
$(u_{f},u,u_{r})$ = (single, double, quad); see also Table 4.

In Table 5, we show data for experiments with $(u_{f},u,u_{r})$ = (half,
single, double). In this case, as predicted in Table 3, SIR fails to converge
once $\kappa_{\infty}(A)>2\cdot 10^{3}$. Moreover, it is seen that TSIR
switches to SGMRES-IR even for relatively well-conditioned matrices for which
SIR still converges (i.e., for $\kappa_{2}(A)=10^{2}$). Although this switch
is technically not necessary, we argue that the difference in cost versus SIR
is minimal (two additional triangular solves at the working precision).
SGMRES-IR works well up to condition number $\kappa_{\infty}(A)=2\cdot
10^{6}$, and thus TSIR only performs one switch in these cases.

It is interesting to notice that the theory from [2] says that we can only
expect SGMRES-IR to converge as long as $\kappa_{\infty}(A)<4\cdot 10^{4}$;
see also Table 3. However, as this and other experiments in this work show,
SGMRES-IR as well as GMRES-IR actually often perform beyond the limits given
by the analysis. This serves as further motivation for the multistage TSIR
approach since we cannot rely entirely on the existing theoretical bounds to
determine whether or not an algorithm will be effective.

For the condition number $\kappa_{\infty}(A)=2\cdot 10^{8}$, SGMRES-IR
convergence eventually slows down, and as a result TSIR switches a second time
to GMRES-IR. Beyond this limit, surprisingly, GMRES-IR and thus TSIR still
converge (despite that the theory only guarantees that GMRES-IR and thus TSIR
will work up to condition number $8\cdot 10^{6}$). However, GMRES-IR is
performing up to $n$ iterations per refinement step, so this is likely not a
good option from a performance perspective. We show plots of TSIR convergence
behavior for the cases $\kappa_{2}(A)\in\\{10^{1},10^{2},10^{7},10^{14}\\}$ in
Figure 5.

Table 5: Number of SIR, SGMRES-IR, GMRES-IR, and TSIR steps with the number of GMRES iterations for each SGMRES-IR and GMRES-IR step for random dense matrices (mode 2) with various condition numbers $\kappa_{\infty}(A)$, $\kappa_{2}(A)$, using $(u_{f},u,u_{r})$ = (half, single, double). $\kappa_{\infty}(A)$ | $\kappa_{2}(A)$ | SIR | SGMRES-IR | GMRES-IR | TSIR
---|---|---|---|---|---
$2\cdot 10^{2}$ | $10^{1}$ | 3 | (3) | (3) | 3
$2\cdot 10^{3}$ | $10^{2}$ | 3 | (3) | (3) | 2, (3)
$2\cdot 10^{5}$ | $10^{4}$ | - | (3,4) | (3,4) | 2, (3,4)
$2\cdot 10^{6}$ | $10^{5}$ | - | (3,4) | (3,4) | 2, (3,4)
$2\cdot 10^{8}$ | $10^{7}$ | - | (5,12,9,17,5,5,5) | (5,5) | 2, (5,12), (17)
$2\cdot 10^{10}$ | $10^{9}$ | - | - | (100,100,100,22,11,3) | 2, (5,9), (100,100,100,100,7,15)
$2\cdot 10^{12}$ | $10^{11}$ | - | - | (100,63,10) | 2, (100), (100,100,14)
$2\cdot 10^{15}$ | $10^{14}$ | - | - | (100,63,10) | 2, (100), (100,100,14)

Figure 5: Convergence of errors in TSIR for random dense matrices (mode 2)
with $\kappa_{2}(A)=10^{1}$ (top left), $\kappa_{2}(A)=10^{2}$ (top right),
$\kappa_{2}(A)=10^{7}$ (bottom left), and $\kappa_{2}(A)=10^{14}$ (bottom
right) for $(u_{f},u,u_{r})$ = (half, single, double); see also Table 5.

Finally, in Table 6, we test $(u_{f},u,u_{r})$ = (half, double, quad). Except
for SIR, which again fails to converge once $\kappa_{\infty}(A)>2\cdot
10^{3}$, all algorithms converge for $\kappa_{\infty}(A)<10^{15}$. Even in
cases where SIR converges, the SIR convergence is slow enough that TSIR always
switches to SGMRES-IR. TSIR only performs the second switch to GMRES-IR for
the two most ill-conditioned problems (although this switch is actually
unnecessary at least for $\kappa_{2}(A)=10^{11}$; see also the right plot in
Figure 6). Again, this emphasizes the need for a multistage, adaptive
approach; if we were to go by the theoretical bounds, we would not expect
SGMRES-IR to work beyond $\kappa_{\infty}(A)=3\cdot 10^{7}$, but here it works
well even for a matrix with a condition number five orders of magnitude
greater. We plot convergence trajectories for TSIR for the problems with
$\kappa_{2}(A)=10^{1}$ and $\kappa_{2}(A)=10^{11}$ in Figure 6.

Table 6: Number of SIR, SGMRES-IR, GMRES-IR, and TSIR steps with the number of GMRES iterations for each SGMRES-IR and GMRES-IR step for random dense matrices (mode 2) with various condition numbers $\kappa_{\infty}(A)$, $\kappa_{2}(A)$, using $(u_{f},u,u_{r})$ = (half, double, quad). $\kappa_{\infty}(A)$ | $\kappa_{2}(A)$ | SIR | SGMRES-IR | GMRES-IR | TSIR
---|---|---|---|---|---
$2\cdot 10^{2}$ | $10^{1}$ | 7 | (5,5) | (5,5) | 2, (5)
$2\cdot 10^{3}$ | $10^{2}$ | 9 | (5,5) | (5,5) | 2, (5)
$2\cdot 10^{5}$ | $10^{4}$ | - | (5,6) | (5,6) | 2, (5,6)
$2\cdot 10^{6}$ | $10^{5}$ | - | (5,6) | (5,6) | 2, (5,6)
$2\cdot 10^{8}$ | $10^{7}$ | - | (6,7) | (6,7) | 2, (6,7)
$2\cdot 10^{10}$ | $10^{9}$ | - | (7,8) | (7,8) | 2, (7,8)
$2\cdot 10^{12}$ | $10^{11}$ | - | (8,8,9) | (8,9) | 2, (8,8), (9)
$2\cdot 10^{15}$ | $10^{14}$ | - | (100,100,100,25,10,10) | (100,32,10) | 2, (100), (100,10)

Figure 6: Convergence of errors in TSIR for random dense matrices (mode 2)
with $\kappa_{2}(A)=10^{1}$ (left) and $\kappa_{2}(A)=10^{11}$ (right) for
$(u_{f},u,u_{r})$ = (half, double, quad); see also Table 6.

####  3.1.2 Random dense matrices with geometrically distributed singular
values

We now test the convergence behavior of the IR variants for matrices generated
via the MATLAB command `gallery(’randsvd’,n,kappa(i),3)`. As mentioned, we
expect that these are more challenging problems for the GMRES-based iterative
refinement schemes since they contain clusters of eigenvalues close to the
origin which may fail to be effectively shifted away from the origin by a low
precision preconditioner.

In Table 7 we test precisions $(u_{f},u,u_{r})$ = (single, double, quad).
Here, SIR fails to converge within $i_{max}$ steps once $\kappa_{\infty}(A)$
exceeds $10^{7}$ (SIR does eventually converge for the problem with
$\kappa_{\infty}(A)=7\cdot 10^{7}$, but it requires 13 steps). Once
$\kappa_{\infty}(A)$ exceeds $10^{9}$, the number of GMRES iterations required
per SGMRES-IR and GMRES-IR step begins to increase dramatically. This is
expected since there is a cluster of eigenvalues remaining close to the origin
even after low precision preconditioning. Notice that in these cases, TSIR
performs the second switch to GMRES-IR after just one SGMRES-IR step, since
the number of GMRES iterations per SGMRES-IR step exceeds our specified value
of $k_{max}=0.1n$. This has a significant benefit versus SGMRES-IR for the
most ill-conditioned problem, since TSIR does many fewer GMRES iterations. We
show plots for $\kappa_{2}(A)\in\\{10^{1},10^{4},10^{7},10^{14}\\}$ in Figure
7.

Table 7: Number of SIR, SGMRES-IR, GMRES-IR, and TSIR steps with the number of GMRES iterations for each SGMRES-IR and GMRES-IR step for random dense matrices having geometrically distributed singular values (mode 3) with various condition numbers $\kappa_{\infty}(A)$, $\kappa_{2}(A)$, using $(u_{f},u,u_{r})$ = (single, double, quad). $\kappa_{\infty}(A)$ | $\kappa_{2}(A)$ | SIR | SGMRES-IR | GMRES-IR | TSIR
---|---|---|---|---|---
$2\cdot 10^{2}$ | $10^{1}$ | 2 | (2) | (2) | 2
$10^{3}$ | $10^{2}$ | 2 | (2) | (2) | 2
$9\cdot 10^{4}$ | $10^{4}$ | 3 | (3) | (3) | 2, (3)
$8\cdot 10^{5}$ | $10^{5}$ | 4 | (4) | (4) | 2, (4)
$7\cdot 10^{7}$ | $10^{7}$ | - | (7,7) | (7,7) | 2, (6,7)
$6\cdot 10^{9}$ | $10^{9}$ | - | (22,23,25) | (22,26) | 2, (17), (24,26)
$6\cdot 10^{11}$ | $10^{11}$ | - | (53,54,55) | (53,55) | 2, (53), (54)
$5\cdot 10^{14}$ | $10^{14}$ | - | (84,84,84,85,85) | (84,88) | 2, (84), (84,85)

Figure 7: Convergence of errors in TSIR for random dense matrices having
geometrically distributed singular values (mode 3) with $\kappa_{2}(A)=10^{1}$
(top left), $\kappa_{2}(A)=10^{4}$ (top right), $\kappa_{2}(A)=10^{7}$ (bottom
left), and $\kappa_{2}(A)=10^{14}$ (bottom right) for $(u_{f},u,u_{r})$ =
(single, double, quad); see also Table 7.

The story for $(u_{f},u,u_{r})$ = (half, single, double), shown in Table 8, is
similar. Once the condition number exceeds around $\sqrt{u}$, the number of
GMRES iterations required per refinement step increases significantly. Again,
this means that TSIR performs the second switch to GMRES-IR since $k_{max}$ is
exceeded bin SGMRES-IR. And again, for the most ill-conditioned problem for
which SGMRES-IR converges, TSIR does significantly fewer GMRES iterations in
total than SGMRES-IR (although again we note that cases where $n$ GMRES
iterations are required per step are not practical from a performance
perspective). Notice also that for $\kappa_{2}(A)\geq 10^{5}$, TSIR performs
no SIR steps since an Inf or NaN is detected in the first correction term. We
plot TSIR convergence for $\kappa_{2}(A)\in\\{10^{1},10^{4},10^{7},10^{11}\\}$
in Figure 8.

Table 8: Number of SIR, SGMRES-IR, GMRES-IR, and TSIR steps with the number of GMRES iterations for each SGMRES-IR and GMRES-IR step for random dense matrices having geometrically distributed singular values (mode 3) with various condition numbers $\kappa_{\infty}(A)$, $\kappa_{2}(A)$, using $(u_{f},u,u_{r})$ = (half, single, double). $\kappa_{\infty}(A)$ | $\kappa_{2}(A)$ | SIR | SGMRES-IR | GMRES-IR | TSIR
---|---|---|---|---|---
$2\cdot 10^{2}$ | $10^{1}$ | 3 | (3) | (3) | 2, (3)
$10^{3}$ | $10^{2}$ | 4 | (4) | (4) | 2, (4)
$9\cdot 10^{4}$ | $10^{4}$ | - | (13,14) | (13,14) | 2, (11), (13)
$8\cdot 10^{5}$ | $10^{5}$ | - | (38,40) | (38,41) | 0, (38), (40)
$7\cdot 10^{7}$ | $10^{7}$ | - | (100,100,100,100,100,83) | (100,100) | 0, (100), (100,100)
$6\cdot 10^{9}$ | $10^{9}$ | - | - | (100,100,100,100) | 0, (100), (100,100,100,100)
$6\cdot 10^{11}$ | $10^{11}$ | - | - | - | -
$5\cdot 10^{14}$ | $10^{14}$ | - | - | - | -

Figure 8: Convergence of errors in TSIR for random dense matrices having
geometrically distributed singular values (mode 3) with $\kappa_{2}(A)=10^{1}$
(top left), $\kappa_{2}(A)=10^{4}$ (top right), $\kappa_{2}(A)=10^{7}$ (bottom
left), and $\kappa_{2}(A)=10^{11}$ (bottom right) for $(u_{f},u,u_{r})$ =
(half, single, double); see also Table 8.

In Table 9 we show results for $(u_{f},u,u_{r})$ = (half, double, quad). Here
the increase in the number of GMRES iterations per GMRES-based refinement step
as $\kappa(A)$ increases is also dramatic. One interesting thing to notice
here and even moreso in the previous examples with mode 3 matrices is that, in
cases where both SGMRES-IR and GMRES-IR converge, the extra precision used in
GMRES-IR does not help improve the convergence rate of GMRES (nearly the same
number of GMRES iterations per refinement step are required). The aspect that
improves is the number of refinement steps. This is likely related to the
GMRES convergence trajectories for this class of problems; the residual will
nearly stagnate (or at least not make much progress) for a number of
iterations and then convergence will happen suddenly. The difference is likely
that when this convergence does happen, the extra precision in GMRES-IR
results in the approximate solution being found to greater accuracy than in
the uniform precision iterations in SGMRES-IR.

Note also that for the most ill-conditioned matrix, no single-stage algorithm
converges, but mysteriously, TSIR does. Looking at the convergence
trajectories for SGMRES-IR and GMRES-IR, the forward error stagnates just
above $u$, so it is likely that TSIR was just lucky. In any case, based on the
bounds in Table 3 we do not expect any of the variants to converge for this
problem. We show TSIR convergence for $\kappa_{2}(A)=10^{1}$ and
$\kappa_{2}(A)=10^{14}$ in Figure 9.

Table 9: Number of SIR, SGMRES-IR, GMRES-IR, and TSIR steps with the number of GMRES iterations for each SGMRES-IR and GMRES-IR step for random dense matrices having geometrically distributed singular values (mode 3) with various condition numbers $\kappa_{\infty}(A)$, $\kappa_{2}(A)$, using $(u_{f},u,u_{r})$ = (half, double, quad). $\kappa_{\infty}(A)$ | $\kappa_{2}(A)$ | SIR | SGMRES-IR | GMRES-IR | TSIR
---|---|---|---|---|---
$2\cdot 10^{2}$ | $10^{1}$ | 7 | (5,5) | (5,5) | 2, (5)
$10^{3}$ | $10^{2}$ | 9 | (6,6) | (6,6) | 2, (6,6)
$9\cdot 10^{4}$ | $10^{4}$ | - | (19,20) | (19,20) | 2, (19), (20)
$8\cdot 10^{5}$ | $10^{5}$ | - | (43,46) | (43,46) | 0, (43), (46)
$7\cdot 10^{7}$ | $10^{7}$ | - | (83,85) | (83,85) | 0, (83), (85)
$6\cdot 10^{9}$ | $10^{9}$ | - | (100,100) | (100,100) | 0, (100), (100)
$6\cdot 10^{11}$ | $10^{11}$ | - | (100,100,100) | (100,100) | 0, (100), (100,100)
$5\cdot 10^{14}$ | $10^{14}$ | - | - | - | 0, (100), (100,100,100)

Figure 9: Convergence of errors in TSIR for random dense matrices having
geometrically distributed singular values (mode 3) with $\kappa_{2}(A)=10^{1}$
(left), and $\kappa_{2}(A)=10^{14}$ (right) for $(u_{f},u,u_{r})$ = (half,
double, quad).

###  3.2 SuiteSparse Matrices

We now test a subset of matrices taken from the SuiteSparse Collection [6]
whose properties are listed in Table 10. As before, we test SIR, SGMRES-IR,
GMRES-IR, and TSIR with precisions $(u_{f},u,u_{r})$ = (single, double, quad),
$(u_{f},u,u_{r})$ = (half, single, double), and $(u_{f},u,u_{r})$ = (half,
double, quad). For these problems, we always set the right-hand side $b$ to be
the vector of ones. Tables 11-13 show the convergence behavior of the
different algorithmic variants. For each precision combination, we show a few
interesting convergence trajectories in Figures 10-12.

In Table 11 for $(u_{f},u,u_{r})$ = (single, double, quad), it is observed
that SIR converges for all matrices except ww_36_pmec_36 (the most ill-
conditioned one). Since SIR converges in few steps, TSIR does not switch to
SGMRES-IR. For ww_36_pmec_36, however, TSIR detects nonconvergence of SIR and
switches to SGMRES-IR, which then converges in 3 steps.

Table 10: Matrices from [6] used for numerical experiments and their
properties.

Name | Size | Nnz | $\kappa_{\infty}$ | Group | Kind
---|---|---|---|---|---
cage6 | 93 | 785 | 2.34E+01 | vanHeukelum | Directed Weighted Graph
tols90 | 90 | 1746 | 3.14E+04 | Bai | Computational Fluid Dynamics Problem
bfwa62 | 62 | 450 | 1.54E+03 | Bai | Electromagnetics Problem
cage5 | 37 | 233 | 2.91E+01 | vanHeukelum | Directed Weighted Graph
d_dyn | 87 | 230 | 8.71E+06 | Grund | Chemical Process Simulation Problem
d_ss | 53 | 144 | 6.02E+08 | Grund | Chemical Process Simulation Problem
Hamrle1 | 32 | 98 | 5.51E+05 | Hamrle | Circuit Simulation Problem
ww_36_pmec_36 | 66 | 1194 | 4.283E+11 | Rommes | Eigenvalue/Model Reduction Problem
steam3 | 80 | 314 | 7.64E+10 | HB | Computational Fluid Dynamics Problem

Table 11: Number of SIR, SGMRES-IR, GMRES-IR, and TSIR steps with the number of GMRES iterations for each SGMRES-IR and GMRES-IR step for real matrices with $(u_{f},u,u_{r})$ = (single, double, quad). Matrix | SIR | SGMRES-IR | GMRES-IR | TSIR
---|---|---|---|---
cage6 | 2 | (2) | (2) | 2
tols90 | 2 | (2) | (2) | 2
bfwa62 | 2 | (2) | (2) | 2
cage5 | 2 | (2) | (2) | 2
Hamrle1 | 2 | (2) | (2) | 2
d_ss | 2 | (2) | (2) | 2
d_dyn | 2 | (2) | (2) | 2
ww_36_pmec_36 | - | (2,3,3) | (2,3,3) | 2, (2,3,3)
steam3 | 2 | (2) | (2) | 2

Figure 10: Convergence of errors in TSIR for cage6 (left) and ww_36_pmec_36
(right) using $(u_{f},u,u_{r})$ = (single, double, quad).

In Table 12 for $(u_{f},u,u_{r})$ = (half, single, double), it is seen that
SIR converges for the matrices with $\kappa_{\infty}(A)<10^{6}$. It is also
observed that SGMRES-IR fails to converge for the extremely ill-conditioned
matrix ww_36_pmec_36. TSIR switches to GMRES-IR for the two most ill-
conditioned matrices, although for the steam3 case, the second switch is not
technically necessary. It is not clear why for the steam3 matrix, TSIR
performs two SGMRES-IR steps and one GMRES-IR step when SGMRES-IR alone only
requires one step.

Table 12: Number of SIR, SGMRES-IR, GMRES-IR, and TSIR steps with the number of GMRES iterations for each SGMRES-IR and GMRES-IR step for real matrices with $(u_{f},u,u_{r})$ = (half, single, double). Matrix | SIR | SGMRES-IR | GMRES-IR | TSIR
---|---|---|---|---
cage6 | 2 | (3) | (3) | 2
tols90 | 2 | (2) | (2) | 2
bfwa62 | 4 | (3) | (3) | 2, (3)
cage5 | 2 | (3) | (3) | 2
Hamrle1 | 2 | (2) | (2) | 2
d_ss | - | (3,3) | (3,3) | 0, (3,3)
d_dyn | - | (2,2) | (2,2) | 0, (2,2)
ww_36_pmec_36 | - | - | (66,66,7) | 0, (66), (66,66,7)
steam3 | - | (2) | (2) | 2, (2,2), (2)

Figure 11: Convergence of errors in TSIR for cage5 (top left), bwfa62 (top
right), steam3 (bottom left), and ww_36_pmec_36 (bottom right) using
$(u_{f},u,u_{r})$ = (half, single, double).

In Table 13 using $(u_{f},u,u_{r})$ = (half, double, quad), we see that TSIR
is able to converge for all test problems and switches at least to SGMRES-IR
in every case due to slow convergence of SIR. For the matrix ww_36_pmec_36,
TSIR makes the second switch to GMRES-IR after one step due to the constraint
on the number of GMRES iterations per step.

Overall, we note that for these real-world sparse problems, to an even greater
extent than in the dense test cases, SGMRES-IR and GMRES-IR often still
perform well for matrices with condition numbers greatly exceeding constraints
given in Table 3. We again stress that this motivates our multistage approach,
since even if we did know the condition number of our matrix, the theoretical
analysis cannot necessarily give indication of which algorithm variant will be
the best choice.

Table 13: Number of SIR, SGMRES-IR, GMRES-IR, and TSIR steps with the number of GMRES iterations for each SGMRES-IR and GMRES-IR step for real matrices with $(u_{f},u,u_{r})$ = (half, double, quad).  Matrix | SIR | SGMRES-IR | GMRES-IR | TSIR
---|---|---|---|---
cage6 | 5 | (4,4) | (4,4) | 2, (4)
tols90 | 5 | (3,3) | (3,3) | 2, (3)
bfwa62 | 9 | (4,5) | (4,5) | 3, (4)
cage5 | 5 | (4,4) | (4,4) | 2, (3)
Hamrle1 | 5 | (3,3) | (3,3) | 2, (3)
d_ss | - | (4,5) | (4,5) | 0, (4,5)
d_dyn | - | (4,4) | (4,4) | 0, (4,4)
ww_36_pmec_36 | - | (8,8) | (8,9) | 0, (8), (8)
steam3 | - | (3,3) | (3,3) | 2, (3,3,3)

Figure 12: Convergence of errors in TSIR for cage6 (left) and ww_36_pmec_36
(right) using $(u_{f},u,u_{r})$ = (half, double, quad).

##  4 Conclusions and future work

Multiprecision iterative refinement, and in particular, GMRES-based iterative
refinement, have seen great use recently. The freedom to choose different
precision combinations as well as a particular solver leads to an explosion of
potential iterative refinement variants, each of which has different
performance costs and different condition number constraints under which
convergence is guaranteed. In practice, it is often the case that a user does
not know the condition number of their matrix, and thus selecting the best
iterative refinement variant to use may require much trial and error.
Motivated by improving usability and reliability, we have developed a
multistage mixed precision iterative refinement approach, called TSIR. The
TSIR algorithm switches between three different iterative refinement variants,
in order of increasing cost, according to criteria adapted from [7]. We
discuss particulars of the algorithm and perform extensive numerical
experiments on both random dense matrices and matrices from SuiteSparse [6]
using a variety of precision combinations. Our experiments demonstrate that
the algorithmic variants often outperform what is dictated by the theoretical
condition number constraints. This further motivates a multistage approach,
since even if the user does know something about the conditioning of their
matrix, it is not necessarily obvious which variant will be optimal.

There are a number of potential extensions to the work presented here. First,
while we have only performed experiments for the case where the factorization
precision, working precision, and residual precision are different, it is
worthwhile to experiment for the two-precision or fixed-precision case as
well. It is also worthwhile to perform experiments using non-IEEE floating
point formats such as bfloat16 [13]. Our multistage approach is applicable to
these cases. We also plan to extend the multistage approach to least squares
problems, which can be accomplished by combining the error bounds derived in
[8] and the three-precision iterative refinement schemes for least squares
problems in [5]. Finally, while we have roughly used the number of triangular
solves in a given precision as a metric for discussing probable performance
characteristics, this may not be a good indication of performance in practice.
Thorough high-performance experiments on modern GPUs are needed to
appropriately gauge the overhead of the TSIR approach and the algorithm’s
suitability in practice.

## References

  * [1] Advanpix LLC., Multiprecision computing toolbox for MATLAB, http://www.advanpix.com/.
  * [2] P. Amestoy, A. Buttari, N. J. Higham, J.-Y. L’Excellent, T. Mary, and B. Vieublé, Five-precision GMRES-based iterative refinement, MIMS EPrint 2021.5, Manchester Institute for Mathematical Sciences, The University of Manchester, Manchester, UK, Apr. 2021.
  * [3] E. Carson and N. J. Higham, A new analysis of iterative refinement and its application to accurate solution of ill-conditioned sparse linear systems, SIAM Journal on Scientific Computing, 39 (2017), pp. A2834–A2856.
  * [4] E. Carson and N. J. Higham, Accelerating the solution of linear systems by iterative refinement in three precisions, SIAM Journal on Scientific Computing, 40 (2018), pp. A817–A847.
  * [5] E. Carson, N. J. Higham, and S. Pranesh, Three-precision GMRES-based iterative refinement for least squares problems, SIAM Journal on Scientific Computing, 42 (2020), pp. A4063–A4083.
  * [6] T. A. Davis and Y. Hu, The University of Florida Sparse Matrix Collection, ACM Transactions on Mathematical Software, 38 (2011).
  * [7] J. Demmel, Y. Hida, W. Kahan, X. S. Li, S. Mukherjee, and E. J. Riedy, Error bounds from extra-precise iterative refinement, ACM Trans. Math. Softw., 32 (2006), p. 325–351.
  * [8] J. Demmel, Y. Hida, E. J. Riedy, and X. S. Li, Extra-precise iterative refinement for overdetermined least squares problems, ACM Transactions on Mathematical Software (TOMS), 35 (2009), pp. 1–32.
  * [9] A. Haidar, S. Tomov, J. Dongarra, and N. J. Higham, Harnessing GPU tensor cores for fast FP16 arithmetic to speed up mixed-precision iterative refinement solvers, in SC18: International Conference for High Performance Computing, Networking, Storage and Analysis, 2018, pp. 603–613.
  * [10] N. J. Higham and S. Pranesh, Simulating low precision floating-point arithmetic, SIAM Journal on Scientific Computing, 41 (2019), pp. C585–C602.
  * [11] N. J. Higham, S. Pranesh, and M. Zounon, Squeezing a matrix into half precision, with an application to solving linear systems, SIAM Journal on Scientific Computing, 41 (2019), pp. A2536–A2551.
  * [12] HPL-AI mixed-precision benchmark.  https://icl.bitbucket.io/hpl-ai/, November 2019.
  * [13] Intel Corporation, Bfloat16 – hardware numerics definition, Tech. Report 338302-001US, Revision 1.0, Intel, November 2018, https://software.intel.com/en-us/download/bfloat16-hardware-numerics-definition.
  * [14] M. Jankowski and H. Woźniakowski, Iterative refinement implies numerical stability, BIT Numerical Mathematics, 17 (1977), pp. 303–311.
  * [15] S. Kudo, K. Nitadori, T. Ina, and T. Imamura, Prompt report on exa-scale HPL-AI benchmark, in 2020 IEEE Int. Conf. Cluster Comput. (CLUSTER), IEEE, 2020, pp. 418–419.
  * [16] J. Langou, J. Langou, P. Luszczek, J. Kurzak, A. Buttari, and J. Dongarra, Exploiting the performance of 32 bit floating point arithmetic in obtaining 64 bit accuracy (revisiting iterative refinement for linear systems), in SC’06: Proceedings of the 2006 ACM/IEEE conference on Supercomputing, IEEE, 2006, pp. 50–50.
  * [17] C. B. Moler, Iterative refinement in floating point, J. Assoc. Comput. Mach., 14 (1967), pp. 316–321.
  * [18] C. C. Paige, M. Rozložník, and Z. Strakoš, Modified Gram-Schmidt (MGS), least squares, and backward stability of MGS-GMRES, SIAM Journal on Matrix Analysis and Applications, 28 (2006), pp. 264–284.
  * [19] Top500.  Online, June 2021, https://www.top500.org/.
  * [20] J. H. Wilkinson, Progress report on the automatic computing engine, Tech. Report MA/17/1024, Mathematics Division, Department of Scientific and Industrial Research, National Physical Laboratory, Teddington, UK, 1948.
  * [21] J. H. Wilkinson, Rounding errors in algebraic processes, Notes on Applied Science, 32 (1963).