change Chebyshev to be in accordance with docs (reltol wrt initial residual instead of RHS)

Author: ranocha

src/chebyshev.jl

@@ -88,10 +88,7 @@ function chebyshev_iterable!(x, A, b, λmin::Real, λmax::Real;
         r .-= c
     end
     resnorm = norm(r)
-    # TODO: According to the docs, the code below should use the initial residual
-    #       instead of the norm of the RHS `b` to set the relative tolerance.
-    # tolerance = max(reltol * resnorm, abstol)
-    tolerance = max(reltol * norm(b), abstol)
+    tolerance = max(reltol * resnorm, abstol)
 
     ChebyshevIterable(Pl, A, x, r, u, c,
         zero(real(T)),

test/chebyshev.jl

@@ -25,6 +25,7 @@ Random.seed!(1234321)
     A = randSPD(T, n)
     b = rand(T, n)
     reltol = √(eps(real(T)))
+    abstol = reltol
 
     @testset "Without preconditioner" begin
         λ_min, λ_max = approx_eigenvalue_bounds(A)
@@ -38,11 +39,19 @@ Random.seed!(1234321)
     @testset "With an initial guess" begin
         λ_min, λ_max = approx_eigenvalue_bounds(A)
         x0 = rand(T, n)
+        initial_residual = norm(A * x0 - b)
         x, history = chebyshev!(x0, A, b, λ_min, λ_max, reltol=reltol, maxiter=10n, log=true)
         @test isa(history, ConvergenceHistory)
         @test history.isconverged
         @test x == x0
-        @test norm(A * x - b) / norm(b) ≤ reltol
+        @test norm(A * x - b) ≤ reltol * initial_residual
+
+        x0 = rand(T, n)
+        x, history = chebyshev!(x0, A, b, λ_min, λ_max, abstol=abstol, reltol=zero(real(T)), maxiter=10n, log=true)
+        @test isa(history, ConvergenceHistory)
+        @test history.isconverged
+        @test x == x0
+        @test norm(A * x - b) ≤ 2*abstol
     end
 
     @testset "With a preconditioner" begin
@@ -54,4 +63,31 @@ Random.seed!(1234321)
         @test norm(A * x - 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]
+        λ_min, λ_max = approx_eigenvalue_bounds(A)
+
+        # 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 = chebyshev!(x, A, b, λ_min, λ_max, 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 = chebyshev!(x, A, b, λ_min, λ_max, abstol=2*initial_residual, reltol=zero(real(T)), log=true)
+        @test niters(ch) == 0
+    end
+end
 end