From 6975534cb2e3070ea8dc403ee12a6e874eb07482 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Tue, 31 Jan 2017 14:22:42 +0100
Subject: [PATCH] bug #478: fix regression in the eigen decomposition of zero
 matrices. (grafted from 53026d29d41e81065b28631445e8eb5c4044c187 )

---
 Eigen/src/Eigenvalues/ComplexEigenSolver.h |  6 ++++--
 Eigen/src/Eigenvalues/RealSchur.h          | 12 ++++++++++++
 test/eigensolver_complex.cpp               |  9 +++++++++
 test/eigensolver_generic.cpp               |  9 +++++++++
 test/eigensolver_selfadjoint.cpp           |  9 +++++++++
 5 files changed, 43 insertions(+), 2 deletions(-)

diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
index ec3b1633e..dc5fae06a 100644
--- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h
+++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -250,7 +250,7 @@ template<typename _MatrixType> class ComplexEigenSolver
     EigenvectorType m_matX;
 
   private:
-    void doComputeEigenvectors(const RealScalar& matrixnorm);
+    void doComputeEigenvectors(RealScalar matrixnorm);
     void sortEigenvalues(bool computeEigenvectors);
 };
 
@@ -284,10 +284,12 @@ ComplexEigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool
 
 
 template<typename MatrixType>
-void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(const RealScalar& matrixnorm)
+void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm)
 {
   const Index n = m_eivalues.size();
 
+  matrixnorm = numext::maxi(matrixnorm,(std::numeric_limits<RealScalar>::min)());
+
   // Compute X such that T = X D X^(-1), where D is the diagonal of T.
   // The matrix X is unit triangular.
   m_matX = EigenvectorType::Zero(n, n);
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h
index d6a339f07..f5c86041d 100644
--- a/Eigen/src/Eigenvalues/RealSchur.h
+++ b/Eigen/src/Eigenvalues/RealSchur.h
@@ -248,12 +248,24 @@ template<typename MatrixType>
 template<typename InputType>
 RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU)
 {
+  const Scalar considerAsZero = (std::numeric_limits<Scalar>::min)();
+
   eigen_assert(matrix.cols() == matrix.rows());
   Index maxIters = m_maxIters;
   if (maxIters == -1)
     maxIters = m_maxIterationsPerRow * matrix.rows();
 
   Scalar scale = matrix.derived().cwiseAbs().maxCoeff();
+  if(scale<considerAsZero)
+  {
+    m_matT.setZero(matrix.rows(),matrix.cols());
+    if(computeU)
+      m_matU.setIdentity(matrix.rows(),matrix.cols());
+    m_info = Success;
+    m_isInitialized = true;
+    m_matUisUptodate = computeU;
+    return *this;
+  }
 
   // Step 1. Reduce to Hessenberg form
   m_hess.compute(matrix.derived()/scale);
diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp
index 8e2bb9ef0..293b1b265 100644
--- a/test/eigensolver_complex.cpp
+++ b/test/eigensolver_complex.cpp
@@ -131,6 +131,15 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
     ComplexEigenSolver<MatrixType> eig(a.adjoint() * a);
     eig.compute(a.adjoint() * a);
   }
+
+  // regression test for bug 478
+  {
+    a.setZero();
+    ComplexEigenSolver<MatrixType> ei3(a);
+    VERIFY_IS_EQUAL(ei3.info(), Success);
+    VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(),RealScalar(1));
+    VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity());
+  }
 }
 
 template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)
diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp
index e18fbf687..d0e644d4b 100644
--- a/test/eigensolver_generic.cpp
+++ b/test/eigensolver_generic.cpp
@@ -76,6 +76,15 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
     EigenSolver<MatrixType> eig(a.adjoint() * a);
     eig.compute(a.adjoint() * a);
   }
+
+  // regression test for bug 478
+  {
+    a.setZero();
+    EigenSolver<MatrixType> ei3(a);
+    VERIFY_IS_EQUAL(ei3.info(), Success);
+    VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(),RealScalar(1));
+    VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity());
+  }
 }
 
 template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)
diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp
index 4ed126116..39ad4130e 100644
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -180,6 +180,15 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
     SelfAdjointEigenSolver<MatrixType> eig(a.adjoint() * a);
     eig.compute(a.adjoint() * a);
   }
+
+  // regression test for bug 478
+  {
+    a.setZero();
+    SelfAdjointEigenSolver<MatrixType> ei3(a);
+    VERIFY_IS_EQUAL(ei3.info(), Success);
+    VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(),RealScalar(1));
+    VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity());
+  }
 }
 
 template<int>