From fde53e4167a3e9b05f708827b705d5ba7e8bb980 Mon Sep 17 00:00:00 2001
From: Alexander Plavin <alexander@plav.in>
Date: Sat, 13 Nov 2021 13:18:40 +0300
Subject: [PATCH] use Symmetric() to get rid of complex eigen() results

---
 src/bounds/ellipsoid.jl | 4 +++-
 test/bounds/bounds.jl   | 7 +++++++
 2 files changed, 10 insertions(+), 1 deletion(-)

diff --git a/src/bounds/ellipsoid.jl b/src/bounds/ellipsoid.jl
index 6cfac229..3ff7d13f 100644
--- a/src/bounds/ellipsoid.jl
+++ b/src/bounds/ellipsoid.jl
@@ -38,7 +38,9 @@ volume(ell::Ellipsoid) = ell.volume
 # Returns the principal axes
 axes(ell::Ellipsoid) = ell.axes
 
-function decompose(A::AbstractMatrix)
+decompose(A::AbstractMatrix) = decompose(Symmetric(A))  # ensure that eigen() always returns real values
+
+function decompose(A::Symmetric)
     E = eigen(A)
     axlens = @. 1 / sqrt(E.values)
     axes = E.vectors * Diagonal(axlens)
diff --git a/test/bounds/bounds.jl b/test/bounds/bounds.jl
index 5d625a9f..d6bb2bdd 100644
--- a/test/bounds/bounds.jl
+++ b/test/bounds/bounds.jl
@@ -12,6 +12,13 @@ const BOUNDST = [
     Bounds.MultiEllipsoid
 ]
 
+@testset "pathological cases" begin
+    A_almost_symmetric = [2.6081830533175096e8 -5.4107420917559285e6 -1.9314298704966028e9 -2.360066561768968e9; -5.410742091755895e6 379882.440454782 6.715028007245775e7 2.0195280814040575e7; -1.931429870496611e9 6.715028007245693e7 9.811342987452753e10 -4.6579127705367036e7; -2.3600665617689605e9 2.0195280814042665e7 -4.6579127705418006e7 9.80946804720486e10]
+    # shouldn't fail:
+    ell = Bounds.Ellipsoid(zeros(4), A_almost_symmetric)
+    Bounds.volume(ell)
+end
+
 @testset "interface - $B, $T, D=$D" for B in BOUNDST, T in [Float32, Float64], D in 1:20
     # creation, inspection
     bound = B(T, D)