Add info() method which can be queried to check whether iteration converged.

Eigen/src/Eigenvalues/ComplexEigenSolver.h

@@ -27,6 +27,9 @@
 #ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
 #define EIGEN_COMPLEX_EIGEN_SOLVER_H
 
+#include "./EigenvaluesCommon.h"
+#include "./ComplexSchur.h"
+
 /** \eigenvalues_module \ingroup Eigenvalues_Module
   * \nonstableyet
   *
@@ -220,6 +223,16 @@ template<typename _MatrixType> class ComplexEigenSolver
       */
     ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
 
+    /** \brief Reports whether previous computation was successful.
+      *
+      * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
+      */
+    ComputationInfo info() const
+    {
+      ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
+      return m_schur.info();
+    }
+
   protected:
     EigenvectorType m_eivec;
     EigenvalueType m_eivalues;
@@ -243,11 +256,14 @@ ComplexEigenSolver<MatrixType>& ComplexEigenSolver<MatrixType>::compute(const Ma
   // Do a complex Schur decomposition, A = U T U^*
   // The eigenvalues are on the diagonal of T.
   m_schur.compute(matrix, computeEigenvectors);
-  m_eivalues = m_schur.matrixT().diagonal();
 
-  if(computeEigenvectors)
-    doComputeEigenvectors(matrix.norm());
-  sortEigenvalues(computeEigenvectors);
+  if(m_schur.info() == Success) 
+  {
+    m_eivalues = m_schur.matrixT().diagonal();
+    if(computeEigenvectors)
+      doComputeEigenvectors(matrix.norm());
+    sortEigenvalues(computeEigenvectors);
+  }
 
   m_isInitialized = true;
   m_eigenvectorsOk = computeEigenvectors;

Eigen/src/Eigenvalues/ComplexSchur.h

@@ -27,6 +27,9 @@
 #ifndef EIGEN_COMPLEX_SCHUR_H
 #define EIGEN_COMPLEX_SCHUR_H
 
+#include "./EigenvaluesCommon.h"
+#include "./HessenbergDecomposition.h"
+
 template<typename MatrixType, bool IsComplex> struct ei_complex_schur_reduce_to_hessenberg;
 
 /** \eigenvalues_module \ingroup Eigenvalues_Module
@@ -192,6 +195,16 @@ template<typename _MatrixType> class ComplexSchur
       */
     ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
 
+    /** \brief Reports whether previous computation was successful.
+      *
+      * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
+      */
+    ComputationInfo info() const
+    {
+      ei_assert(m_isInitialized && "RealSchur is not initialized.");
+      return m_info;
+    }
+
     /** \brief Maximum number of iterations.
       *
       * Maximum number of iterations allowed for an eigenvalue to converge. 
@@ -201,6 +214,7 @@ template<typename _MatrixType> class ComplexSchur
   protected:
     ComplexMatrixType m_matT, m_matU;
     HessenbergDecomposition<MatrixType> m_hess;
+    ComputationInfo m_info;
     bool m_isInitialized;
     bool m_matUisUptodate;
 
@@ -312,6 +326,7 @@ ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& ma
   {
     m_matT = matrix.template cast<ComplexScalar>();
     if(computeU)  m_matU = ComplexMatrixType::Identity(1,1);
+    m_info = Success;
     m_isInitialized = true;
     m_matUisUptodate = computeU;
     return *this;
@@ -382,7 +397,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
 
     // if we spent too many iterations on the current element, we give up
     iter++;
-    if(iter >= m_maxIterations) break;
+    if(iter > m_maxIterations) break;
 
     // find il, the top row of the active submatrix
     il = iu-1;
@@ -412,12 +427,10 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
     }
   }
 
-  if(iter >= m_maxIterations) 
-  {
-    // FIXME : what to do when iter==MAXITER ??
-    // std::cerr << "MAXITER" << std::endl;
-    return;
-  }
+  if(iter <= m_maxIterations) 
+    m_info = Success;
+  else
+    m_info = NoConvergence;
 
   m_isInitialized = true;
   m_matUisUptodate = computeU;

Eigen/src/Eigenvalues/EigenSolver.h

@@ -26,6 +26,7 @@
 #ifndef EIGEN_EIGENSOLVER_H
 #define EIGEN_EIGENSOLVER_H
 
+#include "./EigenvaluesCommon.h"
 #include "./RealSchur.h"
 
 /** \eigenvalues_module \ingroup Eigenvalues_Module
@@ -286,6 +287,12 @@ template<typename _MatrixType> class EigenSolver
       */
     EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
 
+    ComputationInfo info() const
+    {
+      ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
+      return m_realSchur.info();
+    }
+
   private:
     void doComputeEigenvectors();
 
@@ -358,33 +365,36 @@ EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matr
 
   // Reduce to real Schur form.
   m_realSchur.compute(matrix, computeEigenvectors);
-  m_matT = m_realSchur.matrixT();
-  if (computeEigenvectors)
-    m_eivec = m_realSchur.matrixU();
-
-  // Compute eigenvalues from matT
-  m_eivalues.resize(matrix.cols());
-  Index i = 0;
-  while (i < matrix.cols()) 
+  if (m_realSchur.info() == Success)
   {
-    if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) 
-    {
-      m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
-      ++i;
-    }
-    else
-    {
-      Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
-      Scalar z = ei_sqrt(ei_abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
-      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);
-      i += 2;
-    }
-  }
+    m_matT = m_realSchur.matrixT();
+    if (computeEigenvectors)
+      m_eivec = m_realSchur.matrixU();
   
-  // Compute eigenvectors.
-  if (computeEigenvectors)
-    doComputeEigenvectors();
+    // Compute eigenvalues from matT
+    m_eivalues.resize(matrix.cols());
+    Index i = 0;
+    while (i < matrix.cols()) 
+    {
+      if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) 
+      {
+        m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
+        ++i;
+      }
+      else
+      {
+        Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
+        Scalar z = ei_sqrt(ei_abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
+        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);
+        i += 2;
+      }
+    }
+    
+    // Compute eigenvectors.
+    if (computeEigenvectors)
+      doComputeEigenvectors();
+  }
 
   m_isInitialized = true;
   m_eigenvectorsOk = computeEigenvectors;

Eigen/src/Eigenvalues/EigenvaluesCommon.h

@@ -0,0 +1,39 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_EIGENVALUES_COMMON_H
+#define EIGEN_EIGENVALUES_COMMON_H
+
+/** \eigenvalues_module \ingroup Eigenvalues_Module 
+  * \nonstableyet
+  *
+  * \brief Enum for reporting the status of a computation. 
+  */
+enum ComputationInfo {
+  Success = 0,       /**< \brief Computation was successful. */
+  NoConvergence = 1  /**< \brief Iterative procedure did not converge. */
+};
+
+#endif // EIGEN_EIGENVALUES_COMMON_H
+

Eigen/src/Eigenvalues/RealSchur.h

@@ -26,6 +26,7 @@
 #ifndef EIGEN_REAL_SCHUR_H
 #define EIGEN_REAL_SCHUR_H
 
+#include "./EigenvaluesCommon.h"
 #include "./HessenbergDecomposition.h"
 
 /** \eigenvalues_module \ingroup Eigenvalues_Module
@@ -176,6 +177,16 @@ template<typename _MatrixType> class RealSchur
       */
     RealSchur& compute(const MatrixType& matrix, bool computeU = true);
 
+    /** \brief Reports whether previous computation was successful.
+      *
+      * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
+      */
+    ComputationInfo info() const
+    {
+      ei_assert(m_isInitialized && "RealSchur is not initialized.");
+      return m_info;
+    }
+
     /** \brief Maximum number of iterations.
       *
       * Maximum number of iterations allowed for an eigenvalue to converge. 
@@ -188,6 +199,7 @@ template<typename _MatrixType> class RealSchur
     MatrixType m_matU;
     ColumnVectorType m_workspaceVector;
     HessenbergDecomposition<MatrixType> m_hess;
+    ComputationInfo m_info;
     bool m_isInitialized;
     bool m_matUisUptodate;
 
@@ -249,20 +261,21 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix,
     {
       Vector3s firstHouseholderVector, shiftInfo;
       computeShift(iu, iter, exshift, shiftInfo);
-      iter = iter + 1;   // (Could check iteration count here.)
-      if (iter >= m_maxIterations) break;
+      iter = iter + 1; 
+      if (iter > m_maxIterations) break;
       Index im;
       initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
       performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
     }
   } 
 
-  if(iter < m_maxIterations) 
-  {
-    m_isInitialized = true;
-    m_matUisUptodate = computeU;
-  }
+  if(iter <= m_maxIterations) 
+    m_info = Success;
+  else
+    m_info = NoConvergence;
 
+  m_isInitialized = true;
+  m_matUisUptodate = computeU;
   return *this;
 }
 

Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h

@@ -26,6 +26,9 @@
 #ifndef EIGEN_SELFADJOINTEIGENSOLVER_H
 #define EIGEN_SELFADJOINTEIGENSOLVER_H
 
+#include "./EigenvaluesCommon.h"
+#include "./Tridiagonalization.h"
+
 /** \eigenvalues_module \ingroup Eigenvalues_Module
   * \nonstableyet
   *
@@ -360,6 +363,16 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
       return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
     }
 
+    /** \brief Reports whether previous computation was successful.
+      *
+      * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
+      */
+    ComputationInfo info() const
+    {
+      ei_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
+      return m_info;
+    }
+
     /** \brief Maximum number of iterations.
       *
       * Maximum number of iterations allowed for an eigenvalue to converge. 
@@ -371,6 +384,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     RealVectorType m_eivalues;
     TridiagonalizationType m_tridiag;
     typename TridiagonalizationType::SubDiagonalType m_subdiag;
+    ComputationInfo m_info;
     bool m_isInitialized;
     bool m_eigenvectorsOk;
 };
@@ -410,6 +424,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(
     m_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0));
     if(computeEigenvectors)
       m_eivec.setOnes();
+    m_info = Success;
     m_isInitialized = true;
     m_eigenvectorsOk = computeEigenvectors;
     return *this;
@@ -443,7 +458,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(
 
     // if we spent too many iterations on the current element, we give up
     iter++;
-    if(iter >= m_maxIterations) break;
+    if(iter > m_maxIterations) break;
 
     start = end - 1;
     while (start>0 && m_subdiag[start-1]!=0)
@@ -452,23 +467,26 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(
     ei_tridiagonal_qr_step(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
   }
 
-  if(iter >= m_maxIterations) 
-  {
-    return *this;
-  }
+  if (iter <= m_maxIterations) 
+    m_info = Success;
+  else
+    m_info = NoConvergence;
 
   // Sort eigenvalues and corresponding vectors.
   // TODO make the sort optional ?
   // TODO use a better sort algorithm !!
-  for (Index i = 0; i < n-1; ++i)
+  if (m_info == Success)
   {
-    Index k;
-    m_eivalues.segment(i,n-i).minCoeff(&k);
-    if (k > 0)
+    for (Index i = 0; i < n-1; ++i)
     {
-      std::swap(m_eivalues[i], m_eivalues[k+i]);
-      if(computeEigenvectors)
-	m_eivec.col(i).swap(m_eivec.col(k+i));
+      Index k;
+      m_eivalues.segment(i,n-i).minCoeff(&k);
+      if (k > 0)
+      {
+        std::swap(m_eivalues[i], m_eivalues[k+i]);
+        if(computeEigenvectors)
+  	m_eivec.col(i).swap(m_eivec.col(k+i));
+      }
     }
   }