From dce8c1c2869178cef26ac581e52bc6d7406ad19e Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sun, 30 Apr 2023 19:53:56 +0200 Subject: [PATCH] Improve accuracy of 2x2 symmetric/Hermitian eigendecomposition (#1158) * Improve accuracy of 2x2 symmetric/Hermitian eigendecomposition * bump version * Update src/eigen.jl Co-authored-by: Yuto Horikawa --------- Co-authored-by: Yuto Horikawa --- Project.toml | 2 +- src/eigen.jl | 11 ++--------- test/eigen.jl | 5 +++++ 3 files changed, 8 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index 8273252c..73efa9ee 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "StaticArrays" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.5.23" +version = "1.5.24" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/eigen.jl b/src/eigen.jl index a50304cf..58920db2 100644 --- a/src/eigen.jl +++ b/src/eigen.jl @@ -159,14 +159,10 @@ end @inline function _eig(::Size{(2,2)}, A::LinearAlgebra.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real} 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. + tmp = norm(SVector((a[1] - a[4])/2, a[3]')) vals = SVector(t_half - tmp, t_half + tmp) v11 = vals[1] - a[4] @@ -187,10 +183,7 @@ 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. + tmp = norm(SVector((a[1] - a[4])/2, a[2])) vals = SVector(t_half - tmp, t_half + tmp) v11 = vals[1] - a[4] diff --git a/test/eigen.jl b/test/eigen.jl index 1a79f59d..a2b01a48 100644 --- a/test/eigen.jl +++ b/test/eigen.jl @@ -93,6 +93,11 @@ using StaticArrays, Test, LinearAlgebra @test (@inferred SVector{2,ComplexF64} eigvals(Symmetric(m1), Symmetric(m2))) ≈ eigvals(Symmetric(m1_a), Symmetric(m2_a)) end + @testset "2×2, difficult case" begin + m1 = SA[1.6590891025248637 -2.7087777909606814e-7; -2.7087777909606814e-7 1.659089317183428] + @test norm(m1 - Array(eigen(m1))) < 1e-15 + end + @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])