Improve 2x2 eigen

src/eigen.jl

@@ -156,59 +156,40 @@ end
     a = A.data
     TA = eltype(A)
 
-    @inbounds if A.uplo == 'U'
-        if !iszero(a[3]) # A is not diagonal
-            t_half = real(a[1] + a[4]) / 2
-            d = real(a[1] * a[4] - a[3]' * a[3]) # Should be real
-
-            tmp2 = t_half * t_half - d
-            tmp = tmp2 < 0 ? zero(tmp2) : sqrt(tmp2) # Numerically stable for identity matrices, etc.
-            vals = SVector(t_half - tmp, t_half + tmp)
-
-            v11 = vals[1] - a[4]
-            n1 = sqrt(v11' * v11 + a[3]' * a[3])
+    A11 = real(a[1])
+    A22 = real(a[4])
+    if A.uplo == 'U'
+        @inbounds A21 = a[3]'
+    else
+        @inbounds A21 = a[2]
+    end
+    @inbounds if !iszero(A21) # A is not diagonal
+        t_half = (A11 + A22) / 2
+        d_half = (A11 - A22) / 2
+
+        tmp = hypot(d_half, A21)
+        vals = SVector(t_half - tmp, t_half + tmp)
+
+        v11 = d_half - tmp
+        n1 = hypot(v11, A21) # always > 0
+        if n1 < floatmin(T) # n1 subnormal
+            scale = inv(floatmin(T))
+            n1 *= scale
+            v11 = (v11 * scale) / n1
+            v12 = (A21 * scale) / n1
+        else
             v11 = v11 / n1
-            v12 = a[3]' / n1
-
-            v21 = vals[2] - a[4]
-            n2 = sqrt(v21' * v21 + a[3]' * a[3])
-            v21 = v21 / n2
-            v22 = a[3]' / n2
-
-            vecs = @SMatrix [ v11  v21 ;
-                              v12  v22 ]
-
-            return Eigen(vals, vecs)
+            v12 = A21 / n1
         end
-    else # A.uplo == 'L'
-        if !iszero(a[2]) # A is not diagonal
-            t_half = real(a[1] + a[4]) / 2
-            d = real(a[1] * a[4] - a[2]' * a[2]) # Should be real
-
-            tmp2 = t_half * t_half - d
-            tmp = tmp2 < 0 ? zero(tmp2) : sqrt(tmp2) # Numerically stable for identity matrices, etc.
-            vals = SVector(t_half - tmp, t_half + tmp)
-
-            v11 = vals[1] - a[4]
-            n1 = sqrt(v11' * v11 + a[2]' * a[2])
-            v11 = v11 / n1
-            v12 = a[2] / n1
 
-            v21 = vals[2] - a[4]
-            n2 = sqrt(v21' * v21 + a[2]' * a[2])
-            v21 = v21 / n2
-            v22 = a[2] / n2
+        vecs = @SMatrix [ v11  -v12' ;
+                          v12  v11' ]
 
-            vecs = @SMatrix [ v11  v21 ;
-                              v12  v22 ]
-
-            return Eigen(vals,vecs)
-        end
+        return Eigen(vals, vecs)
     end
 
-    # A must be diagonal if we reached this point; treatment of uplo 'L' and 'U' is then identical
-    A11 = real(a[1])
-    A22 = real(a[4])
+    # A must be diagonal if we reached this point;
+    # treatment of uplo 'L' and 'U' is then identical
     if A11 < A22
         vals = SVector(A11, A22)
         vecs = @SMatrix [convert(TA, 1) convert(TA, 0);

test/eigen.jl

@@ -75,24 +75,31 @@ using StaticArrays, Test, LinearAlgebra
         @test vals::SVector ≈ sort(m_d)
         @test eigvals(m_c) ≈ sort(m_d)
         @test eigvals(Hermitian(m_c)) ≈ sort(m_d)
+    end
 
-        # issue #523
-        for (i, j) in ((1, 2), (2, 1)), uplo in (:U, :L)
-            A = SMatrix{2,2,Float64}((i, 0, 0, j))
-            E = eigen(Symmetric(A, uplo))
-            @test eigvecs(E) * SDiagonal(eigvals(E)) * eigvecs(E)' ≈ A
-        end
-
-        m1_a = randn(2,2)
-        m1_a = m1_a*m1_a'
-        m1 = SMatrix{2,2}(m1_a)
-        m2_a = randn(2,2)
-        m2_a = m2_a*m2_a'
-        m2 = SMatrix{2,2}(m2_a)
-        @test (@inferred_maybe_allow SVector{2,ComplexF64} eigvals(m1, m2)) ≈ eigvals(m1_a, m2_a)
-        @test (@inferred_maybe_allow SVector{2,ComplexF64} eigvals(Symmetric(m1), Symmetric(m2))) ≈ eigvals(Symmetric(m1_a), Symmetric(m2_a))
+    # issue #523, #694
+    zero = 0.0
+    smallest_non_zero = nextfloat(zero)
+    smallest_normal = floatmin(zero)
+    largest_subnormal = prevfloat(smallest_normal)
+    epsilon = eps(1.0)
+    one_p_epsilon = nextfloat(1.0)
+    degenerate = (zero, -1, 1, smallest_non_zero, smallest_normal, largest_subnormal, epsilon, one_p_epsilon, -one_p_epsilon)
+    @testset "2×2 degenerate cases" for (i, j, k) in Iterators.product(degenerate,degenerate,degenerate), uplo in (:U, :L)
+        A = SMatrix{2,2,Float64}((i, k, k, j))
+        E = eigen(Symmetric(A, uplo))
+        @test eigvecs(E) * SDiagonal(eigvals(E)) * eigvecs(E)' ≈ A
     end
 
+    m1_a = randn(2,2)
+    m1_a = m1_a*m1_a'
+    m1 = SMatrix{2,2}(m1_a)
+    m2_a = randn(2,2)
+    m2_a = m2_a*m2_a'
+    m2 = SMatrix{2,2}(m2_a)
+    @test (@inferred_maybe_allow SVector{2,ComplexF64} eigvals(m1, m2)) ≈ eigvals(m1_a, m2_a)
+    @test (@inferred_maybe_allow SVector{2,ComplexF64} eigvals(Symmetric(m1), Symmetric(m2))) ≈ eigvals(Symmetric(m1_a), Symmetric(m2_a))
+
     @test_throws DimensionMismatch eigvals(SA[1 2 3; 4 5 6], SA[1 2 3; 4 5 5])
     @test_throws DimensionMismatch eigvals(SA[1 2; 4 5], SA[1 2 3; 4 5 5; 3 4 5])