Merge pull request #280 from ranocha/hr/abstol_reltol

Absolute and relative tolerance for linear solvers

Author: haampie

src/bicgstabl.jl

@@ -13,7 +13,7 @@ mutable struct BiCGStabIterable{precT, matT, solT, vecT <: AbstractVector, small
 
     max_mv_products::Int
     mv_products::Int
-    reltol::realT
+    tol::realT
     residual::realT
 
     Pl::precT
@@ -25,11 +25,12 @@ mutable struct BiCGStabIterable{precT, matT, solT, vecT <: AbstractVector, small
 end
 
 function bicgstabl_iterator!(x, A, b, l::Int = 2;
-    Pl = Identity(),
-    max_mv_products = size(A, 2),
-    initial_zero = false,
-    tol = sqrt(eps(real(eltype(b))))
-)
+                             Pl = Identity(),
+                             max_mv_products = size(A, 2),
+                             abstol::Real = zero(real(eltype(b))),
+                             reltol::Real = sqrt(eps(real(eltype(b)))),
+                             tol = nothing, # TODO: Deprecations introduced in v0.8
+                             initial_zero = false)
     T = eltype(x)
     n = size(A, 1)
     mv_products = 0
@@ -41,6 +42,12 @@ function bicgstabl_iterator!(x, A, b, l::Int = 2;
 
     residual = view(rs, :, 1)
 
+    # TODO: Deprecations introduced in v0.8
+    if tol !== nothing
+        Base.depwarn("The keyword argument `tol` is deprecated, use `reltol` instead.", :bicgstabl_iterator!)
+        reltol = tol
+    end
+
     # Compute the initial residual rs[:, 1] = b - A * x
     # Avoid computing A * 0.
     if initial_zero
@@ -62,17 +69,17 @@ function bicgstabl_iterator!(x, A, b, l::Int = 2;
     # For the least-squares problem
     M = zeros(T, l + 1, l + 1)
 
-    # Stopping condition based on relative tolerance.
-    reltol = nrm * tol
+    # Stopping condition based on absolute and relative tolerance.
+    tolerance = max(reltol * nrm, abstol)
 
     BiCGStabIterable(A, l, x, r_shadow, rs, us,
-        max_mv_products, mv_products, reltol, nrm,
+        max_mv_products, mv_products, tolerance, nrm,
         Pl,
         γ, ω, σ, M
     )
 end
 
-@inline converged(it::BiCGStabIterable) = it.residual ≤ it.reltol
+@inline converged(it::BiCGStabIterable) = it.residual ≤ it.tol
 @inline start(::BiCGStabIterable) = 0
 @inline done(it::BiCGStabIterable, iteration::Int) = it.mv_products ≥ it.max_mv_products || converged(it)
 
@@ -157,10 +164,16 @@ bicgstabl(A, b, l = 2; kwargs...) = bicgstabl!(zerox(A, b), A, b, l; initial_zer
 - `max_mv_products::Int = size(A, 2)`: maximum number of matrix vector products.
 For BiCGStab(l) this is a less dubious term than "number of iterations";
 - `Pl = Identity()`: left preconditioner of the method;
-- `tol::Real = sqrt(eps(real(eltype(b))))`: tolerance for stopping condition `|r_k| / |r_0| ≤ tol`.
-   Note that (1) the true residual norm is never computed during the iterations,
-   only an approximation; and (2) if a preconditioner is given, the stopping condition is based on the
-   *preconditioned residual*.
+- `abstol::Real = zero(real(eltype(b)))`,
+  `reltol::Real = sqrt(eps(real(eltype(b))))`: absolute and relative
+  tolerance for the stopping condition
+  `|r_k| / |r_0| ≤ max(reltol * resnorm, abstol)`, where `r_k = A * x_k - b`
+  is the residual in the `k`th iteration;
+  !!! note
+      1. The true residual norm is never computed during the iterations,
+         only an approximation;
+      2. If a left preconditioner is given, the stopping condition is based on the
+         *preconditioned residual*.
 
 # Return values
 
@@ -174,21 +187,31 @@ For BiCGStab(l) this is a less dubious term than "number of iterations";
 - `history`: convergence history.
 """
 function bicgstabl!(x, A, b, l = 2;
-    tol = sqrt(eps(real(eltype(b)))),
-    max_mv_products::Int = size(A, 2),
-    log::Bool = false,
-    verbose::Bool = false,
-    Pl = Identity(),
-    kwargs...
-)
+                    abstol::Real = zero(real(eltype(b))),
+                    reltol::Real = sqrt(eps(real(eltype(b)))),
+                    tol = nothing, # TODO: Deprecations introduced in v0.8
+                    max_mv_products::Int = size(A, 2),
+                    log::Bool = false,
+                    verbose::Bool = false,
+                    Pl = Identity(),
+                    kwargs...)
     history = ConvergenceHistory(partial = !log)
-    history[:tol] = tol
+    history[:abstol] = abstol
+    history[:reltol] = reltol
 
     # This doesn't yet make sense: the number of iters is smaller.
     log && reserve!(history, :resnorm, max_mv_products)
 
-    # Actually perform CG
-    iterable = bicgstabl_iterator!(x, A, b, l; Pl = Pl, tol = tol, max_mv_products = max_mv_products, kwargs...)
+    # TODO: Deprecations introduced in v0.8
+    if tol !== nothing
+        Base.depwarn("The keyword argument `tol` is deprecated, use `reltol` instead.", :bicgstabl!)
+        reltol = tol
+    end
+
+    # Actually perform iterative solve
+    iterable = bicgstabl_iterator!(x, A, b, l; Pl = Pl,
+                                   abstol = abstol, reltol = reltol,
+                                   max_mv_products = max_mv_products, kwargs...)
 
     if log
         history.mvps = iterable.mv_products

src/cg.jl

@@ -8,7 +8,7 @@ mutable struct CGIterable{matT, solT, vecT, numT <: Real}
     r::vecT
     c::vecT
     u::vecT
-    reltol::numT
+    tol::numT
     residual::numT
     prev_residual::numT
     maxiter::Int
@@ -22,14 +22,14 @@ mutable struct PCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Numb
     r::vecT
     c::vecT
     u::vecT
-    reltol::numT
+    tol::numT
     residual::numT
     ρ::paramT
     maxiter::Int
     mv_products::Int
 end
 
-@inline converged(it::Union{CGIterable, PCGIterable}) = it.residual ≤ it.reltol
+@inline converged(it::Union{CGIterable, PCGIterable}) = it.residual ≤ it.tol
 
 @inline start(it::Union{CGIterable, PCGIterable}) = 0
 
@@ -41,7 +41,10 @@ end
 ###############
 
 function iterate(it::CGIterable, iteration::Int=start(it))
-    if done(it, iteration) return nothing end
+    # Check for termination first
+    if done(it, iteration)
+        return nothing
+    end
 
     # u := r + βu (almost an axpy)
     β = it.residual^2 / it.prev_residual^2
@@ -72,6 +75,7 @@ function iterate(it::PCGIterable, iteration::Int=start(it))
         return nothing
     end
 
+    # Apply left preconditioner
     ldiv!(it.c, it.Pl, it.r)
 
     ρ_prev = it.ρ
@@ -114,40 +118,44 @@ struct CGStateVariables{T,Tx<:AbstractArray{T}}
 end
 
 function cg_iterator!(x, A, b, Pl = Identity();
-    tol = sqrt(eps(real(eltype(b)))),
-    maxiter::Int = size(A, 2),
-    statevars::CGStateVariables = CGStateVariables(zero(x), similar(x), similar(x)),
-    initially_zero::Bool = false
-)
+                      abstol::Real = zero(real(eltype(b))),
+                      reltol::Real = sqrt(eps(real(eltype(b)))),
+                      tol = nothing, # TODO: Deprecations introduced in v0.8
+                      maxiter::Int = size(A, 2),
+                      statevars::CGStateVariables = CGStateVariables(zero(x), similar(x), similar(x)),
+                      initially_zero::Bool = false)
     u = statevars.u
     r = statevars.r
     c = statevars.c
     u .= zero(eltype(x))
     copyto!(r, b)
 
+    # TODO: Deprecations introduced in v0.8
+    if tol !== nothing
+        Base.depwarn("The keyword argument `tol` is deprecated, use `reltol` instead.", :cg_iterator!)
+        reltol = tol
+    end
+
     # Compute r with an MV-product or not.
     if initially_zero
         mv_products = 0
-        c = similar(x)
-        residual = norm(b)
-        reltol = residual * tol # Save one dot product
     else
         mv_products = 1
         mul!(c, A, x)
         r .-= c
-        residual = norm(r)
-        reltol = norm(b) * tol
     end
+    residual = norm(r)
+    tolerance = max(reltol * residual, abstol)
 
     # Return the iterable
     if isa(Pl, Identity)
         return CGIterable(A, x, r, c, u,
-            reltol, residual, one(residual),
+            tolerance, residual, one(residual),
             maxiter, mv_products
         )
     else
         return PCGIterable(Pl, A, x, r, c, u,
-            reltol, residual, one(eltype(x)),
+            tolerance, residual, one(eltype(x)),
             maxiter, mv_products
         )
     end
@@ -177,7 +185,11 @@ cg(A, b; kwargs...) = cg!(zerox(A, b), A, b; initially_zero = true, kwargs...)
   residual vector;
 - `Pl = Identity()`: left preconditioner of the method. Should be symmetric,
   positive-definite like `A`;
-- `tol::Real = sqrt(eps(real(eltype(b))))`: tolerance for stopping condition `|r_k| / |r_0| ≤ tol`;
+- `abstol::Real = zero(real(eltype(b)))`,
+  `reltol::Real = sqrt(eps(real(eltype(b))))`: absolute and relative
+  tolerance for the stopping condition
+  `|r_k| / |r_0| ≤ max(reltol * resnorm, abstol)`, where `r_k = A * x_k - b`
+  is the residual in the `k`th iteration;
 - `maxiter::Int = size(A,2)`: maximum number of iterations;
 - `verbose::Bool = false`: print method information;
 - `log::Bool = false`: keep track of the residual norm in each iteration.
@@ -199,20 +211,29 @@ cg(A, b; kwargs...) = cg!(zerox(A, b), A, b; initially_zero = true, kwargs...)
 - `:resnom` => `::Vector`: residual norm at each iteration.
 """
 function cg!(x, A, b;
-    tol = sqrt(eps(real(eltype(b)))),
-    maxiter::Int = size(A, 2),
-    log::Bool = false,
-    statevars::CGStateVariables = CGStateVariables(zero(x), similar(x), similar(x)),
-    verbose::Bool = false,
-    Pl = Identity(),
-    kwargs...
-)
+             abstol::Real = zero(real(eltype(b))),
+             reltol::Real = sqrt(eps(real(eltype(b)))),
+             tol = nothing, # TODO: Deprecations introduced in v0.8
+             maxiter::Int = size(A, 2),
+             log::Bool = false,
+             statevars::CGStateVariables = CGStateVariables(zero(x), similar(x), similar(x)),
+             verbose::Bool = false,
+             Pl = Identity(),
+             kwargs...)
     history = ConvergenceHistory(partial = !log)
-    history[:tol] = tol
+    history[:abstol] = abstol
+    history[:reltol] = reltol
     log && reserve!(history, :resnorm, maxiter + 1)
 
+    # TODO: Deprecations introduced in v0.8
+    if tol !== nothing
+        Base.depwarn("The keyword argument `tol` is deprecated, use `reltol` instead.", :cg!)
+        reltol = tol
+    end
+
     # Actually perform CG
-    iterable = cg_iterator!(x, A, b, Pl; tol = tol, maxiter = maxiter, statevars = statevars, kwargs...)
+    iterable = cg_iterator!(x, A, b, Pl; abstol = abstol, reltol = reltol, maxiter = maxiter,
+                            statevars = statevars, kwargs...)
     if log
         history.mvps = iterable.mv_products
     end