abstol/reltol for BiCGStab(l)

Author: ranocha

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/chebyshev.jl

@@ -163,7 +163,7 @@ function chebyshev!(x, A, b, λmin::Real, λmax::Real;
 
     # TODO: Deprecations introduced in v0.8
     if tol !== nothing
-        Base.depwarn("The keyword argument `tol` is deprecated, use `reltol` instead.", :cg!)
+        Base.depwarn("The keyword argument `tol` is deprecated, use `reltol` instead.", :chebyshev!)
         reltol = tol
     end
 

src/gmres.jl

@@ -200,7 +200,7 @@ function gmres!(x, A, b;
 
     # TODO: Deprecations introduced in v0.8
     if tol !== nothing
-        Base.depwarn("The keyword argument `tol` is deprecated, use `reltol` instead.", :cg!)
+        Base.depwarn("The keyword argument `tol` is deprecated, use `reltol` instead.", :gmres!)
         reltol = tol
     end
 

src/idrs.jl

@@ -59,7 +59,7 @@ function idrs!(x, A, b;
 
     # TODO: Deprecations introduced in v0.8
     if tol !== nothing
-        Base.depwarn("The keyword argument `tol` is deprecated, use `reltol` instead.", :cg!)
+        Base.depwarn("The keyword argument `tol` is deprecated, use `reltol` instead.", :idrs!)
         reltol = tol
     end
 

test/bicgstabl.jl

@@ -13,20 +13,20 @@ n = 20
     x = ones(T, n)
     b = A * x
 
-    tol = √eps(real(T))
+    reltol = √eps(real(T))
 
     @testset "BiCGStab($l)" for l = (2, 4)
         # Solve without preconditioner
-        x1, his1 = bicgstabl(A, b, l, max_mv_products = 100, log = true, tol = tol)
+        x1, his1 = bicgstabl(A, b, l, max_mv_products = 100, log = true, reltol = reltol)
         @test isa(his1, ConvergenceHistory)
-        @test norm(A * x1 - b) / norm(b) ≤ tol
+        @test norm(A * x1 - b) / norm(b) ≤ reltol
 
         # With an initial guess
         x_guess = rand(T, n)
-        x2, his2 = bicgstabl!(x_guess, A, b, l, max_mv_products = 100, log = true, tol = tol)
+        x2, his2 = bicgstabl!(x_guess, A, b, l, max_mv_products = 100, log = true, reltol = reltol)
         @test isa(his2, ConvergenceHistory)
         @test x2 == x_guess
-        @test norm(A * x2 - b) / norm(b) ≤ tol
+        @test norm(A * x2 - b) / norm(b) ≤ reltol
 
         # The following tests fails CI on Windows and Ubuntu due to a
         # `SingularException(4)`
@@ -35,8 +35,34 @@ n = 20
         end
         # Do an exact LU decomp of a nearby matrix
         F = lu(A + rand(T, n, n))
-        x3, his3 = bicgstabl(A, b, Pl = F, l, max_mv_products = 100, log = true, tol = tol)
-        @test norm(A * x3 - b) / norm(b) ≤ tol
+        x3, his3 = bicgstabl(A, b, Pl = F, l, max_mv_products = 100, log = true, reltol = reltol)
+        @test norm(A * x3 - b) / norm(b) ≤ reltol
+    end
+end
+
+@testset "Termination criterion" begin
+    for T in (Float32, Float64, ComplexF32, ComplexF64)
+        A = T[ 2 -1  0
+              -1  2 -1
+               0 -1  2]
+        n = size(A, 2)
+        b = ones(T, n)
+        x0 = A \ b
+        perturbation = T[(-1)^i for i in 1:n]
+
+        # If the initial residual is small and a small relative tolerance is used,
+        # many iterations are necessary
+        x = x0 + sqrt(eps(real(T))) * perturbation
+        initial_residual = norm(A * x - b)
+        x, ch = bicgstabl!(x, A, b, log=true)
+        @test 1 ≤ niters(ch) ≤ n÷2 # BiCGStab(2) makes more than one RHS evaluations per iteration
+
+        # If the initial residual is small and a large absolute tolerance is used,
+        # no iterations are necessary
+        x = x0 + 10*sqrt(eps(real(T))) * perturbation
+        initial_residual = norm(A * x - b)
+        x, ch = bicgstabl!(x, A, b, abstol=2*initial_residual, reltol=zero(real(T)), log=true)
+        @test niters(ch) == 0
     end
 end
 end