From 3c66bb136bf2adcb9d73d3d66850a8b907bc9264 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Mon, 14 Apr 2014 22:00:27 +0200
Subject: [PATCH] bug #793: detect NaN and INF in EigenSolver instead of
 aborting with an assert.

---
 Eigen/src/Eigenvalues/EigenSolver.h | 25 ++++++++++++++++++++++---
 1 file changed, 22 insertions(+), 3 deletions(-)

diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h
index 739466949..d2563d470 100644
--- a/Eigen/src/Eigenvalues/EigenSolver.h
+++ b/Eigen/src/Eigenvalues/EigenSolver.h
@@ -275,10 +275,11 @@ template<typename _MatrixType> class EigenSolver
       */
     EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
 
+    /** \returns NumericalIssue if the input contains INF or NaN values or overflow occured. Returns Success otherwise. */
     ComputationInfo info() const
     {
       eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
-      return m_realSchur.info();
+      return m_info;
     }
 
     /** \brief Sets the maximum number of iterations allowed. */
@@ -302,6 +303,7 @@ template<typename _MatrixType> class EigenSolver
     EigenvalueType m_eivalues;
     bool m_isInitialized;
     bool m_eigenvectorsOk;
+    ComputationInfo m_info;
     RealSchur<MatrixType> m_realSchur;
     MatrixType m_matT;
 
@@ -367,12 +369,15 @@ EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvect
   using std::sqrt;
   using std::abs;
   using std::max;
+  using numext::isfinite;
   eigen_assert(matrix.cols() == matrix.rows());
 
   // Reduce to real Schur form.
   m_realSchur.compute(matrix, computeEigenvectors);
+  
+  m_info = m_realSchur.info();
 
-  if (m_realSchur.info() == Success)
+  if (m_info == Success)
   {
     m_matT = m_realSchur.matrixT();
     if (computeEigenvectors)
@@ -386,6 +391,13 @@ EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvect
       if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) 
       {
         m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
+        if(!isfinite(m_eivalues.coeffRef(i)))
+        {
+          m_isInitialized = true;
+          m_eigenvectorsOk = false;
+          m_info = NumericalIssue;
+          return *this;
+        }
         ++i;
       }
       else
@@ -406,6 +418,13 @@ EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvect
         
         m_eivalues.coeffRef(i)   = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
         m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
+        if(!(isfinite(m_eivalues.coeffRef(i)) && isfinite(m_eivalues.coeffRef(i+1))))
+        {
+          m_isInitialized = true;
+          m_eigenvectorsOk = false;
+          m_info = NumericalIssue;
+          return *this;
+        }
         i += 2;
       }
     }
@@ -594,7 +613,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
     }
     else
     {
-      eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen
+      eigen_assert(0 && "Internal bug in EigenSolver (INF or NaN has not been detected)"); // this should not happen
     }
   }