Merge pull request #268 from haampie/fix-division-by-zero

Fix division by zero due to exact off-diagonal zeros in Hessenberg matrix

Author: haampie

Project.toml

@@ -1,6 +1,6 @@
 name = "IterativeSolvers"
 uuid = "42fd0dbc-a981-5370-80f2-aaf504508153"
-version = "0.8.0"
+version = "0.8.2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/gmres.jl

@@ -205,10 +205,14 @@ function gmres!(x, A, b;
 end
 
 function update_residual!(r::Residual, arnoldi::ArnoldiDecomp, k::Int)
-    # Cheaply computes the current residual
-    r.nullvec[k + 1] = -conj(dot(view(r.nullvec, 1 : k), view(arnoldi.H, 1 : k, k)) / arnoldi.H[k + 1, k])
-    r.accumulator += abs2(r.nullvec[k + 1])
-    r.current = r.β / √r.accumulator
+    if iszero(arnoldi.H[k + 1, k])
+        r.current = zero(r.current)
+    else
+        # Cheaply computes the current residual
+        r.nullvec[k + 1] = -conj(dot(view(r.nullvec, 1 : k), view(arnoldi.H, 1 : k, k)) / arnoldi.H[k + 1, k])
+        r.accumulator += abs2(r.nullvec[k + 1])
+        r.current = r.β / √r.accumulator
+    end
 end
 
 function init!(arnoldi::ArnoldiDecomp{T}, x, b, Pl, Ax; initially_zero::Bool = false) where {T}

test/gmres.jl

@@ -62,4 +62,11 @@ end
     x = gmres(A, b; tol=tol, maxiter=2000)
     @test norm(A * x - b) / norm(b) ≤ tol
 end
+
+@testset "Off-diagonal in hessenberg matrix exactly zero" begin
+    A = Matrix(1.0I, 2, 2)
+    b = [1.0, 2.2]
+    x = gmres(A, b)
+    @test all(x .== b)
+end
 end