abstol/reltol for GMRES

Author: ranocha

src/gmres.jl

@@ -42,18 +42,21 @@ mutable struct GMRESIterable{preclT, precrT, solT, rhsT, vecT, arnoldiT <: Arnol
     restart::Int
     k::Int
     maxiter::Int
-    reltol::resT
+    tol::resT
     β::resT
 end
 
-converged(g::GMRESIterable) = g.residual.current ≤ g.reltol
+converged(g::GMRESIterable) = g.residual.current ≤ g.tol
 
 start(::GMRESIterable) = 0
 
 done(g::GMRESIterable, iteration::Int) = iteration ≥ g.maxiter || converged(g)
 
 function iterate(g::GMRESIterable, iteration::Int=start(g))
-    if done(g, iteration) return nothing end
+    # Check for termination first
+    if done(g, iteration)
+        return nothing
+    end
 
     # Arnoldi step: expand
     expand!(g.arnoldi, g.Pl, g.Pr, g.k, g.Ax)
@@ -100,15 +103,22 @@ function iterate(g::GMRESIterable, iteration::Int=start(g))
 end
 
 function gmres_iterable!(x, A, b;
-    Pl = Identity(),
-    Pr = Identity(),
-    tol = sqrt(eps(real(eltype(b)))),
-    restart::Int = min(20, size(A, 2)),
-    maxiter::Int = size(A, 2),
-    initially_zero::Bool = false
-)
+                         Pl = Identity(),
+                         Pr = Identity(),
+                         abstol::Real = zero(real(eltype(b))),
+                         reltol::Real = sqrt(eps(real(eltype(b)))),
+                         tol = nothing, # TODO: Deprecations introduced in v0.8
+                         restart::Int = min(20, size(A, 2)),
+                         maxiter::Int = size(A, 2),
+                         initially_zero::Bool = false)
     T = eltype(x)
 
+    # TODO: Deprecations introduced in v0.8
+    if tol !== nothing
+        Base.depwarn("The keyword argument `tol` is deprecated, use `reltol` instead.", :gmres_iterable!)
+        reltol = tol
+    end
+
     # Approximate solution
     arnoldi = ArnoldiDecomp(A, restart, T)
     residual = Residual(restart, T)
@@ -119,11 +129,11 @@ function gmres_iterable!(x, A, b;
     residual.current = init!(arnoldi, x, b, Pl, Ax, initially_zero = initially_zero)
     init_residual!(residual, residual.current)
 
-    reltol = tol * residual.current
+    tolerance = max(reltol * residual.current, abstol)
 
     GMRESIterable(Pl, Pr, x, b, Ax,
         arnoldi, residual,
-        mv_products, restart, 1, maxiter, reltol, residual.current
+        mv_products, restart, 1, maxiter, tolerance, residual.current
     )
 end
 
@@ -150,7 +160,10 @@ Solves the problem ``Ax = b`` with restarted GMRES.
 - `initially_zero::Bool`: If `true` assumes that `iszero(x)` so that one
   matrix-vector product can be saved when computing the initial
   residual vector;
-- `tol`: relative tolerance;
+- `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`
 - `restart::Int = min(20, size(A, 2))`: restarts GMRES after specified number of iterations;
 - `maxiter::Int = size(A, 2)`: maximum number of inner iterations of GMRES;
 - `Pl`: left preconditioner;
@@ -170,20 +183,29 @@ Solves the problem ``Ax = b`` with restarted GMRES.
 - `history`: convergence history.
 """
 function gmres!(x, A, b;
-  Pl = Identity(),
-  Pr = Identity(),
-  tol = sqrt(eps(real(eltype(b)))),
-  restart::Int = min(20, size(A, 2)),
-  maxiter::Int = size(A, 2),
-  log::Bool = false,
-  initially_zero::Bool = false,
-  verbose::Bool = false
-)
+                Pl = Identity(),
+                Pr = Identity(),
+                abstol::Real = zero(real(eltype(b))),
+                reltol::Real = sqrt(eps(real(eltype(b)))),
+                tol = nothing, # TODO: Deprecations introduced in v0.8
+                restart::Int = min(20, size(A, 2)),
+                maxiter::Int = size(A, 2),
+                log::Bool = false,
+                initially_zero::Bool = false,
+                verbose::Bool = false)
     history = ConvergenceHistory(partial = !log, restart = restart)
     history[:tol] = tol
     log && reserve!(history, :resnorm, maxiter)
 
-    iterable = gmres_iterable!(x, A, b; Pl = Pl, Pr = Pr, tol = tol, maxiter = maxiter, restart = restart, initially_zero = initially_zero)
+    # TODO: Deprecations introduced in v0.8
+    if tol !== nothing
+        Base.depwarn("The keyword argument `tol` is deprecated, use `reltol` instead.", :cg!)
+        reltol = tol
+    end
+
+    iterable = gmres_iterable!(x, A, b; Pl = Pl, Pr = Pr,
+                               abstol = abstol, reltol = reltol, maxiter = maxiter,
+                               restart = restart, initially_zero = initially_zero)
 
     verbose && @printf("=== gmres ===\n%4s\t%4s\t%7s\n","rest","iter","resnorm")
 

test/gmres.jl

@@ -15,29 +15,29 @@ n = 10
     A = rand(T, n, n)
     b = rand(T, n)
     F = lu(A)
-    tol = √eps(real(T))
+    reltol = √eps(real(T))
 
     # Test optimality condition: residual should be non-increasing
-    x, history = gmres(A, b, log = true, restart = 3, maxiter = 10, tol = tol);
+    x, history = gmres(A, b, log=true, restart=3, maxiter=10, reltol=reltol);
     @test isa(history, ConvergenceHistory)
     @test all(diff(history[:resnorm]) .<= 0.0)
 
     # Left exact preconditioner
-    x, history = gmres(A, b, Pl=F, maxiter=1, restart=1, tol=tol, log=true)
+    x, history = gmres(A, b, Pl=F, maxiter=1, restart=1, reltol=reltol, log=true)
     @test history.isconverged
-    @test norm(F \ (A * x - b)) / norm(b) ≤ tol
+    @test norm(F \ (A * x - b)) / norm(b) ≤ reltol
 
     # Right exact preconditioner
-    x, history = gmres(A, b, Pl=Identity(), Pr=F, maxiter=1, restart=1, tol=tol, log=true)
+    x, history = gmres(A, b, Pl=Identity(), Pr=F, maxiter=1, restart=1, reltol=reltol, log=true)
     @test history.isconverged
-    @test norm(A * x - b) / norm(b) ≤ tol
+    @test norm(A * x - b) / norm(b) ≤ reltol
 end
 
 @testset "SparseMatrixCSC{$T, $Ti}" for T in (Float64, ComplexF64), Ti in (Int64, Int32)
     A = sprand(T, n, n, 0.5) + I
     b = rand(T, n)
     F = lu(A)
-    tol = √eps(real(T))
+    reltol = √eps(real(T))
 
     # Test optimality condition: residual should be non-increasing
     x, history = gmres(A, b, log = true, restart = 3, maxiter = 10);
@@ -46,21 +46,21 @@ end
     # Left exact preconditioner
     x, history = gmres(A, b, Pl=F, maxiter=1, restart=1, log=true)
     @test history.isconverged
-    @test norm(F \ (A * x - b)) / norm(b) ≤ tol
+    @test norm(F \ (A * x - b)) / norm(b) ≤ reltol
 
     # Right exact preconditioner
     x, history = gmres(A, b, Pl = Identity(), Pr=F, maxiter=1, restart=1, log=true)
     @test history.isconverged
-    @test norm(A * x - b) / norm(b) ≤ tol
+    @test norm(A * x - b) / norm(b) ≤ reltol
 end
 
 @testset "Linear operator defined as a function" begin
     A = LinearMap(cumsum!, 100; ismutating=true)
     b = rand(100)
-    tol = 1e-5
+    reltol = 1e-5
 
-    x = gmres(A, b; tol=tol, maxiter=2000)
-    @test norm(A * x - b) / norm(b) ≤ tol
+    x = gmres(A, b; reltol=reltol, maxiter=2000)
+    @test norm(A * x - b) / norm(b) ≤ reltol
 end
 
 @testset "Off-diagonal in hessenberg matrix exactly zero" begin
@@ -69,4 +69,30 @@ end
     x = gmres(A, b)
     @test all(x .== b)
 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 = gmres!(x, A, b, log=true)
+        @test 2 ≤ niters(ch) ≤ n
+
+        # 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 = gmres!(x, A, b, abstol=2*initial_residual, reltol=zero(real(T)), log=true)
+        @test niters(ch) == 0
+    end
+end
 end