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

Author: haampie

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