bug #1557: fix RealSchur and EigenSolver for matrices with only zeros on the diagonal.

1 files changed, 66 insertions, 8 deletions

diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp
index e0e435151..086ecdf5e 100644
--- a/test/eigensolver_generic.cpp
+++ b/test/eigensolver_generic.cpp
@@ -12,6 +12,21 @@
 #include <limits>
 #include <Eigen/Eigenvalues>
 
+template<typename EigType,typename MatType>
+void check_eigensolver_for_given_mat(const EigType &eig, const MatType& a)
+{
+  typedef typename NumTraits<typename MatType::Scalar>::Real RealScalar;
+  typedef Matrix<RealScalar, MatType::RowsAtCompileTime, 1> RealVectorType;
+  typedef typename std::complex<RealScalar> Complex;
+  Index n = a.rows();
+  VERIFY_IS_EQUAL(eig.info(), Success);
+  VERIFY_IS_APPROX(a * eig.pseudoEigenvectors(), eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix());
+  VERIFY_IS_APPROX(a.template cast<Complex>() * eig.eigenvectors(),
+                   eig.eigenvectors() * eig.eigenvalues().asDiagonal());
+  VERIFY_IS_APPROX(eig.eigenvectors().colwise().norm(), RealVectorType::Ones(n).transpose());
+  VERIFY_IS_APPROX(a.eigenvalues(), eig.eigenvalues());
+}
+
 template<typename MatrixType> void eigensolver(const MatrixType& m)
 {
   /* this test covers the following files:
@@ -22,8 +37,7 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
 
   typedef typename MatrixType::Scalar Scalar;
   typedef typename NumTraits<Scalar>::Real RealScalar;
-  typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
-  typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
+  typedef typename std::complex<RealScalar> Complex;
 
   MatrixType a = MatrixType::Random(rows,cols);
   MatrixType a1 = MatrixType::Random(rows,cols);
@@ -36,12 +50,7 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
     (ei0.pseudoEigenvectors().template cast<Complex>()) * (ei0.eigenvalues().asDiagonal()));
 
   EigenSolver<MatrixType> ei1(a);
-  VERIFY_IS_EQUAL(ei1.info(), Success);
-  VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix());
-  VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(),
-                   ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
-  VERIFY_IS_APPROX(ei1.eigenvectors().colwise().norm(), RealVectorType::Ones(rows).transpose());
-  VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues());
+  CALL_SUBTEST( check_eigensolver_for_given_mat(ei1,a) );
 
   EigenSolver<MatrixType> ei2;
   ei2.setMaxIterations(RealSchur<MatrixType>::m_maxIterationsPerRow * rows).compute(a);
@@ -100,6 +109,19 @@ template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m
   VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors());
 }
 
+
+template<typename CoeffType>
+Matrix<typename CoeffType::Scalar,Dynamic,Dynamic>
+make_companion(const CoeffType& coeffs)
+{
+  Index n = coeffs.size()-1;
+  Matrix<typename CoeffType::Scalar,Dynamic,Dynamic> res(n,n);
+  res.setZero();
+	res.row(0) = -coeffs.tail(n) / coeffs(0);
+	res.diagonal(-1).setOnes();
+  return res;
+}
+
 template<int>
 void eigensolver_generic_extra()
 {
@@ -126,6 +148,42 @@ void eigensolver_generic_extra()
     VERIFY_IS_APPROX((a * eig.eigenvectors()).norm()+1., 1.);
     VERIFY_IS_APPROX((eig.eigenvectors() * eig.eigenvalues().asDiagonal()).norm()+1., 1.);
   }
+
+  // regression test for bug 933
+  {
+    {
+      VectorXd coeffs(5); coeffs << 1, -3, -175, -225, 2250;
+      MatrixXd C = make_companion(coeffs);
+      EigenSolver<MatrixXd> eig(C);
+      CALL_SUBTEST( check_eigensolver_for_given_mat(eig,C) );
+    }
+    {
+      // this test is tricky because it requires high accuracy in smallest eigenvalues
+      VectorXd coeffs(5); coeffs << 6.154671e-15, -1.003870e-10, -9.819570e-01, 3.995715e+03, 2.211511e+08;
+      MatrixXd C = make_companion(coeffs);
+      EigenSolver<MatrixXd> eig(C);
+      CALL_SUBTEST( check_eigensolver_for_given_mat(eig,C) );
+      Index n = C.rows();
+      for(Index i=0;i<n;++i)
+      {
+        typedef std::complex<double> Complex;
+        MatrixXcd ac = C.cast<Complex>();
+        ac.diagonal().array() -= eig.eigenvalues()(i);
+        VectorXd sv = ac.jacobiSvd().singularValues();
+        // comparing to sv(0) is not enough here to catch the "bug",
+        // the hard-coded 1.0 is important!
+        VERIFY_IS_MUCH_SMALLER_THAN(sv(n-1), 1.0);
+      }
+    }
+  }
+  // regression test for bug 1557
+  {
+    // this test is interesting because it contains zeros on the diagonal.
+    MatrixXd A_bug1557(3,3);
+    A_bug1557 << 0, 0, 0, 1, 0, 0.5887907064808635127, 0, 1, 0;
+    EigenSolver<MatrixXd> eig(A_bug1557);
+    CALL_SUBTEST( check_eigensolver_for_given_mat(eig,A_bug1557) );
+  }
 }
 
 EIGEN_DECLARE_TEST(eigensolver_generic)