Use <= when testing for deflation for avoid issue with zeros in

diagonal. Fixes #52
JuliaLinearAlgebra · Apr 17, 2019 · 61faeb5 · 61faeb5
1 parent 31d2d80
commit 61faeb5
Show file tree

Hide file tree

Showing 2 changed files with 14 additions and 2 deletions.
diff --git a/src/eigenSelfAdjoint.jl b/src/eigenSelfAdjoint.jl
@@ -174,7 +174,7 @@ function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T), debug::Bool=false) wher
         while true
             # Check for zero off diagonal elements
             for be in blockstart + 1:n
-                if abs(e[be - 1]) < tol*sqrt(abs(d[be - 1]))*sqrt(abs(d[be]))
+                if abs(e[be - 1]) <= tol*sqrt(abs(d[be - 1]))*sqrt(abs(d[be]))
                     blockend = be - 1
                     break
                 end
@@ -231,7 +231,7 @@ function eigQL!(S::SymTridiagonal{T};
         while true
             # Check for zero off diagonal elements
             for be = blockstart+1:n
-                if abs(e[be-1]) < tol*sqrt(abs(d[be-1]))*sqrt(abs(d[be]))
+                if abs(e[be-1]) <= tol*sqrt(abs(d[be-1]))*sqrt(abs(d[be]))
                     blockend = be - 1
                     break
                 end

diff --git a/test/eigenselfadjoint.jl b/test/eigenselfadjoint.jl
@@ -90,4 +90,16 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0
         @test GenericLinearAlgebra._eigen!(SymTridiagonal([1.0], Float64[])).values == [1.0]
         @test GenericLinearAlgebra._eigen!(SymTridiagonal([1.0], Float64[])).vectors == fill(1.0, 1, 1)
     end
+
+    # Issue #52
+    @testset "convergence criterion when diagonal has zeros" begin
+        M1 = Hermitian(zeros(3, 3))
+        M1[1,1] = 0.01
+        M2 = Hermitian(zeros(3, 3))
+        M2[3, 3] = 0.01
+        @test eigvals(M1) == GenericLinearAlgebra._eigvals!(copy(M1))
+        @test eigvals(M2) == GenericLinearAlgebra._eigvals!(copy(M2))
+        @test eigen(M1).values == GenericLinearAlgebra._eigen!(copy(M1)).values
+        @test eigen(M2).values == GenericLinearAlgebra._eigen!(copy(M2)).values
+    end
 end