Accelerating the Solution of Linear Systems by Iterative Refinement in Three Precisions

Abstract

Abstract. We propose a general algorithm for solving an n\times n nonsingular linear system Ax = b based on iterative refinement with three precisions. The working precision is combined with possibly different precisions for solving for the correction term and for computing the residuals. Via rounding error analysis of the algorithm we derive sufficient conditions for convergence and bounds for the attainable forward error and normwise and componentwise backward errors. Our results generalize and unify many existing rounding error analyses for iterative refinement. With single precision as the working precision, we show that by using LU factorization in IEEE half precision as the solver and calculating the residuals in double precision it is possible to solve Ax = b to full single precision accuracy for \infty -norm condition numbers \kappa \infty (A) \leq 10 4 , with the O(n 3 ) part of the computations carried out entirely in half precision. We show further that by solving the correction equations by GMRES preconditioned by the LU factors the restriction on the condition number can be weakened to \kappa \infty (A) \leq 10 8 , although in general there is no guarantee that GMRES will converge quickly. Taking for comparison a standard Ax = b solver that uses LU factorization in single precision, these results suggest that on architectures for which half precision is efficiently implemented it will be possible to solve certain linear systems Ax = b up to twice as fast and to greater accuracy. Analogous results are given with double precision as the working precision.Key words. iterative refinement, linear system, multiple precision, mixed precision, rounding error analysis, backward error, forward error, GMRES, preconditioning AMS subject classifications. 65G50, 65F10 DOI. 10.1137/17M11408191. Introduction. Iterative refinement is a method for improving an approximate solution y to a linear system Ax = b by computing the residual r = b -Ay, solving the correction equation Ad = r, forming the update y \leftarr y + d, and repeating these steps as necessary. We consider a general iterative refinement algorithm that includes a variety of existing ones as special cases. The algorithm contains three precisions:\bullet u is the precision at which the data A, b and the solution x are stored (the working precision), \bullet u f is the precision at which the factorization of A is computed, \bullet u r is the precision at which residuals are computed. The precisions are assumed to satisfy

Harnessing GPU Tensor Cores for Fast FP16 Arithmetic to Speed up Mixed-Precision Iterative Refinement Solvers

“…Replacing the direct triangular solves of the correction equation with an iterative method, as suggested in [4] in a mixed precision context, leads to "nesting" of two iterative methods, which in general are called "inner-outer" iterations, the latter having been studied both theoretically and computationally [9], [21], [23], including in mixed-precision computation scenarios [2]. Recently, Carson and Higham [4], [5] analyzed the convergence property of a three precision iterative refinement scheme (factorization precision, working precision, residual precision) and concluded that if the condition number of A is not too large, κ ∞ (A) = A ∞ A −1 ∞ < 10 4 , then using FP16 for the O(n 3 ) portion (the LU factorization) and (FP32, FP64) or (FP64, FP128) as the (working, residual) precision for the O(n 2 ) portion (refinement loop), one can expect to achieve forward error and backward error on the order of 10 −8 and 10 −16 respectively. We note that, ifx is the solution of Ax = b the forward error is defined by x − x ∞ / x ∞ and the backward error is defined by r 2 / A 2 x 2 where r = b−Ax.…”

“…The convergence tolerance is chosen of the order of the unit roundoff of the low precision arithmetic used during the factorization (e.g., we use 10 −4 or 10 −8 for when the LU is in FP16 or FP32 respectively). Since this paper focuses on practical usage and possible performance gains rather than error analysis, we point the reader to [4], [5] for detailed error analysis of the IR and IRGM techniques.…”

Adaptive precision in block‐Jacobi preconditioning for iterative sparse linear system solvers

“…preconditioner memory transfers (11) data transfers (from memory) per iteration, where fpxx i denotes the precision format selected for the ith diagonal block of the preconditioner. The data transfer volume of the block-Jacobi preconditioner thus depends on the format employed to store the block inverse.…”

“…To avoid the previous two pitfalls, in our final experiment, we compute the total data transfers of a single iteration of the PCG method with the block-Jacobi preconditioner stored in fp64, fp32, fp16, or adaptive precision, see Equation (11). To obtain an estimated total data transfer volume, we then combine the data transfer volume per iteration with the number of iterations needed to reach convergence in each case, ignoring those cases for which half precision does not converge.…”

“…Carson and Higham provides a detailed error analysis of LU‐based mixed refinement approaches for ill‐conditioned systems. In another work of the aforementioned authors, the same authors go as far as using half precision for computing an LU preconditioner that is used in the solution process of a GMRES solver that is part of a mixed precision iterative refinement process.…”