|
| 1 | +# IterativeRefinement |
| 2 | + |
| 3 | +This package is an implementation of multi-precision iterative refinement for |
| 4 | +certain dense-matrix linear algebra problems. |
| 5 | + |
| 6 | +## Background |
| 7 | +The purpose of iterative refinement (IR) is to improve the accuracy of a |
| 8 | +solution. If `x` is the exact solution of `A*x=b`, a simple solve of |
| 9 | +the form `y = A \ b` will have a relative forward error |
| 10 | +(`norm(y-x)/norm(x)`) of approximately `ϵ * O(n) * cond(A)` where `ϵ` |
| 11 | +is the unit roundoff error in the standard precision. Iterative |
| 12 | +refinement with higher precision residuals can reduce this to |
| 13 | + `ϵ * O(n)`, as long as the matrix `A` is not very badly conditioned |
| 14 | +relative to `ϵ`. |
| 15 | + |
| 16 | +Why not do everything in high precision? The factorization step is |
| 17 | +typically *very* expensive (`O(n^3)`) in high precision, whereas the |
| 18 | +residual computation is relatively cheap (`O(n^2)`). Furthermore, IR |
| 19 | +schemes often provide useful error bounds. |
| 20 | + |
| 21 | +For typical use, one would have a basic working precision of `Float64` |
| 22 | +(`ϵ = 2.2e-16`), so that fast LAPACK/BLAS routines dominate the runtime. |
| 23 | +`rfldiv` will then (by default) use `BigFloat` for residuals. |
| 24 | +One might alternatively use `Double64` from |
| 25 | +[DoubleFloats.jl](https://github.com/JuliaMath/DoubleFloats.jl) |
| 26 | +or `Float128` from |
| 27 | +[Quadmath.jl](https://github.com/JuliaMath/Quadmath.jl) |
| 28 | + |
| 29 | +# Linear systems |
| 30 | + |
| 31 | +This package provides a function `rfldiv`, which |
| 32 | +handles linear matrix-vector problems of the form |
| 33 | + |
| 34 | +`A x = b`. |
| 35 | + |
| 36 | +## Basic Usage |
| 37 | +```julia |
| 38 | +julia> using LinearAlgebra, IterativeRefinement |
| 39 | +julia> x, bnorm, bcomp = rfldiv(A,b) |
| 40 | +``` |
| 41 | +This provides an accurate solution vector `x` and estimated bounds |
| 42 | +on norm-wise and component-wise relative error. By default `LU` decomposition |
| 43 | +is used. |
| 44 | + |
| 45 | +## Advanced Usage |
| 46 | +See the function docstring for details. |
| 47 | + |
| 48 | +If one has several right-hand-sides, one can equilibrate and factor |
| 49 | +`A` in advance; see the tests for an example. |
| 50 | + |
| 51 | +## Reference |
| 52 | +J.Demmel et al., "Error bounds from extra precise iterative refinement," |
| 53 | +LAPACK Working Note Nr. 165 (2005), also published as |
| 54 | +ACM TOMS, 32, 325 (2006). The work |
| 55 | +described therein eventually turned into a collection of subroutines |
| 56 | +included in some versions of LAPACK. This implementation is based on |
| 57 | +the paper; minor modifications were introduced based on experimentation. |
| 58 | +To be precise, this package implements Algorithm 3. |
| 59 | + |
| 60 | +# Eigensystems |
| 61 | + |
| 62 | +Additional methods (`rfeigen`) are provided for improving estimates of |
| 63 | +eigenvalue/subspace pairs of the form |
| 64 | + |
| 65 | +`A X = X λ`. |
| 66 | + |
| 67 | +For a simple eigenvalue, `X` is the corresponding eigenvector, and |
| 68 | +the user provides coarse estimates of both. In the case of |
| 69 | +multiple or defective eigenvalues, columns of `X` are generators for the |
| 70 | +corresponding invariant subspace, and the user provides a Schur decomposition |
| 71 | +with a list of indices for the cluster of interest. |
| 72 | + |
| 73 | +Problem-specific error bound estimates are not yet provided for eigensystems. |
| 74 | + |
| 75 | +## Basic Usage |
| 76 | +### isolated eigenvalue |
| 77 | +```julia |
| 78 | +julia> using LinearAlgebra, IterativeRefinement, Quadmath |
| 79 | +julia> E = eigen(A) |
| 80 | +julia> j = your_index_selection() |
| 81 | +julia> λrefined, xrefined = rfeigen(A, E.vectors[:,j], E.values[j], Float128) |
| 82 | +``` |
| 83 | + |
| 84 | +### eigenvalue cluster |
| 85 | +```julia |
| 86 | +julia> using LinearAlgebra, IterativeRefinement, Quadmath |
| 87 | +julia> S = schur(A) |
| 88 | +julia> S = LinearAlgebra.Schur{ComplexF64}(S) # if eltype of A is real |
| 89 | +julia> idx = findall(abs.(S.values .- target) .< 0.1) |
| 90 | +julia> λrefined, Vrefined = rfeigen(A, S, idx, Float128) |
| 91 | +``` |
| 92 | + |
| 93 | +## Reference |
| 94 | +J.J.Dongarra, C.B.Moler, and J.H.Wilkinson, "Improving the accuracy of computed |
| 95 | +eigenvalues and eigenvectors," SIAM J. Numer. Anal. 20, 23-45 (1983). |
| 96 | + |
0 commit comments