From 3ee06ec52f7356614e58bf4b0868bbbebc7d80f6 Mon Sep 17 00:00:00 2001
From: Antonio Sanchez <cantonios@google.com>
Date: Fri, 16 Feb 2024 13:11:54 -0800
Subject: [PATCH] Fix real schur and polynomial solver.

For certain inputs, the real schur decomposition might get stuck in a cycle.
Exceptional shifts are supposed to knock us out of that - but previously
they were only ever applied at iteration 10 and 30, which doesn't help if
the cycle starts after cycle 30.  Modified to apply a shift every 16 iterations
(for reference, LAPACK seems to do it every 6 iterations).

Also added an assert in polynomial solver to verify that the schur decomposition
was successful.

Fixes #2633.
---
 Eigen/src/Eigenvalues/RealSchur.h             | 53 +++++++++----------
 test/schur_real.cpp                           | 12 +++++
 .../Eigen/src/Polynomials/PolynomialSolver.h  |  1 +
 3 files changed, 39 insertions(+), 27 deletions(-)

diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h
index 7304ef344..37394e191 100644
--- a/Eigen/src/Eigenvalues/RealSchur.h
+++ b/Eigen/src/Eigenvalues/RealSchur.h
@@ -435,34 +435,33 @@ inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& ex
   shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
   shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
 
-  // Wilkinson's original ad hoc shift
-  if (iter == 10)
-  {
-    exshift += shiftInfo.coeff(0);
-    for (Index i = 0; i <= iu; ++i)
-      m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
-    Scalar s = abs(m_matT.coeff(iu,iu-1)) + abs(m_matT.coeff(iu-1,iu-2));
-    shiftInfo.coeffRef(0) = Scalar(0.75) * s;
-    shiftInfo.coeffRef(1) = Scalar(0.75) * s;
-    shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
-  }
-
-  // MATLAB's new ad hoc shift
-  if (iter == 30)
-  {
-    Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
-    s = s * s + shiftInfo.coeff(2);
-    if (s > Scalar(0))
-    {
-      s = sqrt(s);
-      if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
-        s = -s;
-      s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
-      s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
-      exshift += s;
+  // Alternate exceptional shifting strategy every 16 iterations.
+  if (iter % 16 == 0) {
+    // Wilkinson's original ad hoc shift
+    if (iter % 32 != 0) {
+      exshift += shiftInfo.coeff(0);
       for (Index i = 0; i <= iu; ++i)
-        m_matT.coeffRef(i,i) -= s;
-      shiftInfo.setConstant(Scalar(0.964));
+        m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
+      Scalar s = abs(m_matT.coeff(iu,iu-1)) + abs(m_matT.coeff(iu-1,iu-2));
+      shiftInfo.coeffRef(0) = Scalar(0.75) * s;
+      shiftInfo.coeffRef(1) = Scalar(0.75) * s;
+      shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
+    } else {
+      // MATLAB's new ad hoc shift
+      Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
+      s = s * s + shiftInfo.coeff(2);
+      if (s > Scalar(0))
+      {
+        s = sqrt(s);
+        if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
+          s = -s;
+        s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
+        s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
+        exshift += s;
+        for (Index i = 0; i <= iu; ++i)
+          m_matT.coeffRef(i,i) -= s;
+        shiftInfo.setConstant(Scalar(0.964));
+      }
     }
   }
 }
diff --git a/test/schur_real.cpp b/test/schur_real.cpp
index 945461027..7097fbae3 100644
--- a/test/schur_real.cpp
+++ b/test/schur_real.cpp
@@ -98,6 +98,16 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
   }
 }
 
+void test_bug2633() {
+  Eigen::MatrixXd A(4, 4);
+  A << 0,  0,  0, -2,
+       1,  0,  0, -0,
+       0,  1,  0,  2,
+       0,  0,  2, -0;
+  RealSchur<Eigen::MatrixXd> schur(A);
+  VERIFY(schur.info() == Eigen::Success);
+}
+
 EIGEN_DECLARE_TEST(schur_real)
 {
   CALL_SUBTEST_1(( schur<Matrix4f>() ));
@@ -107,4 +117,6 @@ EIGEN_DECLARE_TEST(schur_real)
 
   // Test problem size constructors
   CALL_SUBTEST_5(RealSchur<MatrixXf>(10));
+
+  CALL_SUBTEST_6(( test_bug2633() ));
 }
diff --git a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h
index 5e0ecbb43..874a6f337 100644
--- a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h
+++ b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h
@@ -354,6 +354,7 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
         internal::companion<Scalar,_Deg> companion( poly );
         companion.balance();
         m_eigenSolver.compute( companion.denseMatrix() );
+        eigen_assert(m_eigenSolver.info() == Eigen::Success);
         m_roots = m_eigenSolver.eigenvalues();
         // cleanup noise in imaginary part of real roots:
         // if the imaginary part is rather small compared to the real part