This package is an implementation of multi-precision iterative refinement for certain dense-matrix linear algebra problems.

The purpose of iterative refinement (IR) is to improve the accuracy of a solution. If x is the exact solution of A*x=b, a simple solve of the form y = A \ b will have a relative forward error (norm(y-x)/norm(x)) of approximately ϵ * O(n) * cond(A) where ϵ is the unit roundoff error in the standard precision. Iterative refinement with higher precision residuals can reduce this to ϵ * O(n), as long as the matrix A is not very badly conditioned relative to ϵ.

Why not do everything in high precision? The factorization step is typically very expensive (O(n^3)) in high precision, whereas the residual computation is relatively cheap (O(n^2)). Furthermore, IR schemes often provide useful error bounds.

For typical use, one would have a basic working precision of Float64 (ϵ = 2.2e-16), so that fast LAPACK/BLAS routines dominate the runtime. rfldiv will then (by default) use BigFloat for residuals. One might alternatively use Double64 from DoubleFloats.jl or Float128 from Quadmath.jl

The most mature part of the package provides a function rfldiv, which handles linear matrix-vector problems of the form

A x = b.

julia> using LinearAlgebra, IterativeRefinement
julia> x, bnorm, bcomp = rfldiv(A,b)

This provides an accurate solution vector x and estimated bounds on norm-wise and component-wise relative error. By default LU decomposition is used.

See the function docstring for details.

If one has several right-hand-sides, one can equilibrate and factor A in advance; see the tests for an example.

J.Demmel et al., "Error bounds from extra precise iterative refinement," LAPACK Working Note Nr. 165 (2005), also published as ACM TOMS, 32, 325 (2006). The work described therein eventually turned into a collection of subroutines included in some versions of LAPACK. This implementation is based on the paper; minor modifications were introduced based on experimentation. To be precise, this package implements Algorithm 3.

Additional methods (rfeigen) are provided for improving estimates of eigenvalue/subspace pairs of the form

A X = X λ.

For a simple eigenvalue, X is the corresponding eigenvector, and the user provides coarse estimates of both. In the case of multiple or defective eigenvalues, columns of X are generators for the corresponding invariant subspace, and the user provides a Schur decomposition with a list of indices for the cluster of interest.

Problem-specific error bound estimates are not yet provided for eigensystems.

J.J.Dongarra, C.B.Moler, and J.H.Wilkinson, "Improving the accuracy of computed eigenvalues and eigenvectors," SIAM J. Numer. Anal. 20, 23-45 (1983).