Bounded mixed-precision iterative refinement

Input nonsingular real $A\in {\mathbb{R}}^{n\times n}$, $n\ge 1$, exact data $b\in {\mathbb{R}}^n$, and an initial iterate $x_0$ stored at solve/update precision $u_w$. Also input stored triangular factors $\widehat{L},\widehat{U}$ formed at factorization precision $u_f$, where $\widehat{L}$ is lower triangular, $\widehat{U}$ is upper triangular, both have nonzero diagonals, and $\widehat{A}:=\widehat{L}\widehat{U}$ is the exact real product; thus $\widehat{A}$ is nonsingular. Input residual precision $u_r$, solve/update precision $u_w$, finite update budget $N\ge 0$, tolerance $\tau \ge 0$, a stagnation window $K_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}}\ge 1$, nonnegative constants ${\eta }_r,{\xi }_r,{\eta }_w,{\xi }_w$, and a proposed contraction upper bound $0\le q_{\mathrm{c}\mathrm{e}\mathrm{r}\mathrm{t}}<1$. Let $x_{\ast }=A^{-1}b$. Use the operationwise model of Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add at each stated precision. Entries of the exact data $A,b$ are available as real operands to the residual operations; the stored factors, iterates, residuals, and corrections are used as their exact represented values in those operations.

Certificates have the following fixed format: one real matrix $B\in {\mathbb{R}}^{n\times n}$ before iteration, nonnegative real scalars ${\rho }_k$ for $0\le k\le N$, and nonnegative real scalars ${\omega }_k$ for $0\le k<N$. Each entry may instead be marked missing. All certificate checks use exact real addition, subtraction, multiplication, division by a positive number, nonnegative square roots, and comparisons; they are not rounded operations at $u_r$ or $u_w$. Evaluate finite sums and products in increasing index order. In particular, evaluate $\mathrm{\parallel }v{\mathrm{\parallel }}_2=\sqrt{{\sum }_i{v}_i^2}$ and $\mathrm{\parallel }M{\mathrm{\parallel }}_F=\sqrt{{\sum }_{i,j}{M}_{ij}^2}$ by these exact operations. This is an ideal exact-real certificate verifier, not a claim that arbitrary real inputs admit a finite-bit implementation.

Before iterating, check triangularity and nonzero diagonals of the factors; on failure return inconclusive with $x_0$. A missing or incorrectly shaped $B$ returns inconclusive with $x_0$. Form $\widehat{A}=\widehat{L}\widehat{U}$ exactly and check $B\widehat{A}=I$ entrywise. If this fails, return inconclusive with $x_0$. Set $a_{+}=\mathrm{\parallel }A{\mathrm{\parallel }}_F,a_{-}=a_{+}\mathrm{/}\sqrt{n},h=\mathrm{\parallel }B{\mathrm{\parallel }}_F,d=\mathrm{\parallel }I-BA{\mathrm{\parallel }}_F\mathrm{.}$ Accept the initial certificate precisely when $d+h((1+{\eta }_w){\eta }_r+{\eta }_w)a_{+}\le q_{\mathrm{c}\mathrm{e}\mathrm{r}\mathrm{t}}<1;$ otherwise return inconclusive with $x_0$. Here $B={\widehat{A}}^{-1}$, and the inequalities $a_{-}\le \mathrm{\parallel }A{\mathrm{\parallel }}_2\le a_{+}$ and $\mathrm{\parallel }M{\mathrm{\parallel }}_2\le \mathrm{\parallel }M{\mathrm{\parallel }}_F$ from Spectral and Frobenius norms are unitarily invariant, are given by singular values, and satisfy the sharp rank comparison make this a sufficient certificate for $q=\mathrm{\parallel }I-{\widehat{A}}^{-1}A{\mathrm{\parallel }}_2+\mathrm{\parallel }{\widehat{A}}^{-1}{\mathrm{\parallel }}_2((1+{\eta }_w){\eta }_r+{\eta }_w)\mathrm{\parallel }A{\mathrm{\parallel }}_2\le q_{\mathrm{c}\mathrm{e}\mathrm{r}\mathrm{t}}\mathrm{.}$ These conservative tests may reject data for which sharper norm certificates exist.

Before entering the iteration loop, initialize the consecutive-equality counter $c_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}}=0$ and mark the preceding residual as undefined. Retain this counter and preceding residual across iterations; do not reinitialize them inside the loop. The counter records the number of consecutive equalities between adjacent successfully certified computed residual vectors.

For each $k=0,\dots ,N$, compute the residual by rows in increasing $i$: initialize $s_i=0$, then for $j=1,\dots ,n$ in increasing order compute $p_{ij}={\mathrm{fl}}_{u_r}(A_{ij}{x}_{k,j})$ and $s_i={\mathrm{fl}}_{u_r}(s_i+p_{ij})$; finally set ${\widehat{r}}_{k,i}={\mathrm{fl}}_{u_r}(b_i-s_i)$. Any operationwise range failure returns inconclusive with the current iterate and the first failing row, column, and operation. Define the exact residual defect ${\delta }_{r,k}={\widehat{r}}_k-(b-A{x}_k)$, so ${\widehat{r}}_k=A(x_{\ast }-x_k)+{\delta }_{r,k}$. Check that ${\rho }_k$ is present and nonnegative and that ${\sum }_i{\delta }_{r,k,i}^2\le {\rho }_k^2,{\rho }_k\le {\eta }_r(a_{-}\mathrm{\parallel }{x}_k{\mathrm{\parallel }}_2+\mathrm{\parallel }b{\mathrm{\parallel }}_2)+{\xi }_r\mathrm{.}$ A missing scalar or failed comparison returns inconclusive with the current iterate. These exact tests imply $\mathrm{\parallel }{\delta }_{r,k}{\mathrm{\parallel }}_2\le {\eta }_r(\mathrm{\parallel }A{\mathrm{\parallel }}_2\mathrm{\parallel }{x}_k{\mathrm{\parallel }}_2+\mathrm{\parallel }b{\mathrm{\parallel }}_2)+{\xi }_r$. Use the retained consecutive-equality counter and preceding residual. Use the exact Euclidean norm of the stored residual vector for the status comparison; this is a comparison of the already computed vector, not an additional rounded arithmetic operation. If $\mathrm{\parallel }{\widehat{r}}_k{\mathrm{\parallel }}_2\le \tau$, return residual-tolerance with $x_k$. If a preceding residual exists, compare its components exactly with those of ${\widehat{r}}_k$: if all are equal, increment $c_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}}$ by one; otherwise reset it to zero. If $c_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}}=K_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}}$, return successful-stagnation with $x_k$. This status certifies only that the last $K_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}}$ successful updates produced identical computed residual vectors, with all range and contraction certificates accepted; it certifies neither a small exact error nor convergence. Before any further update, store ${\widehat{r}}_k$ as the preceding residual.

If $k=N$ and neither status above was returned, return budget-exhausted with $x_k$. Otherwise solve $\widehat{L}{y}_k={\widehat{r}}_k$ by forward substitution, then $\widehat{U}{\widehat{z}}_k=y_k$ by backward substitution, using precision $u_w$. For forward substitution, visit rows $i=1,\dots ,n$ in increasing order, initialize $s_i=0$, and for $j=1,\dots ,i-1$ in increasing order compute $p_{ij}={\mathrm{fl}}_{u_w}({\widehat{L}}_{ij}{y}_{k,j})$ and $s_i={\mathrm{fl}}_{u_w}(s_i+p_{ij})$, then $t_i={\mathrm{fl}}_{u_w}({\widehat{r}}_{k,i}-s_i)$ and $y_{k,i}={\mathrm{fl}}_{u_w}(t_i\mathrm{/}{\widehat{L}}_{ii})$. For backward substitution, visit rows $i=n,\dots ,1$ in decreasing order, initialize $s_i=0$, and for $j=i+1,\dots ,n$ in increasing order compute $p_{ij}={\mathrm{fl}}_{u_w}({\widehat{U}}_{ij}{\widehat{z}}_{k,j})$ and $s_i={\mathrm{fl}}_{u_w}(s_i+p_{ij})$, then $t_i={\mathrm{fl}}_{u_w}(y_{k,i}-s_i)$ and ${\widehat{z}}_{k,i}={\mathrm{fl}}_{u_w}(t_i\mathrm{/}{\widehat{U}}_{ii})$. Update components in increasing order by $x_{k+1,i}={\mathrm{fl}}_{u_w}(x_{k,i}+{\widehat{z}}_{k,i})$. A failed operation returns inconclusive with the current iterate and first failing operation location.

After a successful update, define the exact residual-equivalent update defect ${\delta }_{w,k}:=\widehat{A}(x_{k+1}-x_k)-{\widehat{r}}_k,$ so $x_{k+1}=x_k+{\widehat{A}}^{-1}({\widehat{r}}_k+{\delta }_{w,k})$. Check that ${\omega }_k$ is present and nonnegative and that ${\sum }_i{\delta }_{w,k,i}^2\le {\omega }_k^2$ and ${\omega }_k\le {\eta }_w\mathrm{\parallel }{\widehat{r}}_k{\mathrm{\parallel }}_2+{\xi }_w$; otherwise return inconclusive with $x_{k+1}$. These tests imply $\mathrm{\parallel }{\delta }_{w,k}{\mathrm{\parallel }}_2\le {\eta }_w\mathrm{\parallel }{\widehat{r}}_k{\mathrm{\parallel }}_2+{\xi }_w$. Continue until one of the stated statuses is returned. All loops are finite. The procedure does not claim exact termination or unconditional convergence.