Allow user to specify max number of iterations (bug #479).

Eigen/src/Eigenvalues/ComplexEigenSolver.h

@@ -3,7 +3,7 @@
 //
 // Copyright (C) 2009 Claire Maurice
 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -208,7 +208,27 @@ template<typename _MatrixType> class ComplexEigenSolver
       * Example: \include ComplexEigenSolver_compute.cpp
       * Output: \verbinclude ComplexEigenSolver_compute.out
       */
-    ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
+    ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true)
+    {
+      return computeInternal(matrix, computeEigenvectors, false, 0);
+    }
+
+    /** \brief Computes eigendecomposition with specified maximum number of iterations. 
+      * 
+      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
+      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
+      *    eigenvalues are computed; if false, only the eigenvalues are
+      *    computed. 
+      * \param[in]  maxIter  Maximum number of iterations.
+      * \returns    Reference to \c *this
+      *
+      * This method provides the same functionality as compute(const MatrixType&, bool) but it also allows the
+      * user to specify the maximum number of iterations to be used when computing the Schur decomposition.
+      */
+    ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors, Index maxIter)
+    {
+      return computeInternal(matrix, computeEigenvectors, true, maxIter);
+    }
 
     /** \brief Reports whether previous computation was successful.
       *
@@ -231,18 +251,26 @@ template<typename _MatrixType> class ComplexEigenSolver
   private:
     void doComputeEigenvectors(RealScalar matrixnorm);
     void sortEigenvalues(bool computeEigenvectors);
+    ComplexEigenSolver& computeInternal(const MatrixType& matrix, bool computeEigenvectors, 
+                                        bool maxIterSpecified, Index maxIter);
+
 };
 
 
 template<typename MatrixType>
-ComplexEigenSolver<MatrixType>& ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
+ComplexEigenSolver<MatrixType>& 
+ComplexEigenSolver<MatrixType>::computeInternal(const MatrixType& matrix, bool computeEigenvectors, 
+                                                bool maxIterSpecified, Index maxIter)
 {
   // this code is inspired from Jampack
   assert(matrix.cols() == matrix.rows());
 
   // Do a complex Schur decomposition, A = U T U^*
   // The eigenvalues are on the diagonal of T.
-  m_schur.compute(matrix, computeEigenvectors);
+  if (maxIterSpecified)
+    m_schur.compute(matrix, computeEigenvectors, maxIter);
+  else
+    m_schur.compute(matrix, computeEigenvectors);
 
   if(m_schur.info() == Success)
   {

Eigen/src/Eigenvalues/ComplexSchur.h

@@ -3,7 +3,7 @@
 //
 // Copyright (C) 2009 Claire Maurice
 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -166,6 +166,7 @@ template<typename _MatrixType> class ComplexSchur
       * 
       * \param[in]  matrix  Square matrix whose Schur decomposition is to be computed.
       * \param[in]  computeU  If true, both T and U are computed; if false, only T is computed.
+
       * \returns    Reference to \c *this
       *
       * The Schur decomposition is computed by first reducing the
@@ -180,8 +181,27 @@ template<typename _MatrixType> class ComplexSchur
       *
       * Example: \include ComplexSchur_compute.cpp
       * Output: \verbinclude ComplexSchur_compute.out
+      *
+      * \sa compute(const MatrixType&, bool, Index)
       */
-    ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
+    ComplexSchur& compute(const MatrixType& matrix, bool computeU = true)
+    {
+      return compute(matrix, computeU, m_maxIterations * matrix.rows());
+    }
+
+    /** \brief Computes Schur decomposition with specified maximum number of iterations.
+      *
+      * \param[in]  matrix    Square matrix whose Schur decomposition is to be computed.
+      * \param[in]  computeU  If true, both T and U are computed; if false, only T is computed.
+      * \param[in]  maxIter   Maximum number of iterations. 
+      *
+      * \returns    Reference to \c *this
+      *
+      * This method provides the same functionality as compute(const MatrixType&, bool) but it also allows the
+      * user to specify the maximum number of QR iterations to be used. The maximum number of iterations for
+      * compute(const MatrixType&, bool) is #m_maxIterations times the size of the matrix.
+      */
+    ComplexSchur& compute(const MatrixType& matrix, bool computeU, Index maxIter);
 
     /** \brief Reports whether previous computation was successful.
       *
@@ -189,13 +209,14 @@ template<typename _MatrixType> class ComplexSchur
       */
     ComputationInfo info() const
     {
-      eigen_assert(m_isInitialized && "RealSchur is not initialized.");
+      eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
       return m_info;
     }
 
     /** \brief Maximum number of iterations.
       *
-      * Maximum number of iterations allowed for an eigenvalue to converge. 
+      * If not otherwise specified, the maximum number of iterations is this number times the size of the
+      * matrix. It is currently set to 30.
       */
     static const int m_maxIterations = 30;
 
@@ -209,7 +230,7 @@ template<typename _MatrixType> class ComplexSchur
   private:  
     bool subdiagonalEntryIsNeglegible(Index i);
     ComplexScalar computeShift(Index iu, Index iter);
-    void reduceToTriangularForm(bool computeU);
+    void reduceToTriangularForm(bool computeU, Index maxIter);
     friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
 };
 
@@ -268,7 +289,7 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu
 
 
 template<typename MatrixType>
-ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
+ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU, Index maxIter)
 {
   m_matUisUptodate = false;
   eigen_assert(matrix.cols() == matrix.rows());
@@ -284,7 +305,7 @@ ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& ma
   }
 
   internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU);
-  reduceToTriangularForm(computeU);
+  reduceToTriangularForm(computeU, maxIter);
   return *this;
 }
 
@@ -327,7 +348,7 @@ struct complex_schur_reduce_to_hessenberg<MatrixType, false>
 
 // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
 template<typename MatrixType>
-void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
+void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU, Index maxIter)
 {  
   // The matrix m_matT is divided in three parts. 
   // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 
@@ -354,7 +375,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
     // if we spent too many iterations, we give up
     iter++;
     totalIter++;
-    if(totalIter > m_maxIterations * m_matT.cols()) break;
+    if(totalIter > maxIter) break;
 
     // find il, the top row of the active submatrix
     il = iu-1;
@@ -384,7 +405,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
     }
   }
 
-  if(totalIter <= m_maxIterations * m_matT.cols())
+  if(totalIter <= maxIter)
     m_info = Success;
   else
     m_info = NoConvergence;

Eigen/src/Eigenvalues/EigenSolver.h

@@ -2,7 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -273,7 +273,27 @@ template<typename _MatrixType> class EigenSolver
       * Example: \include EigenSolver_compute.cpp
       * Output: \verbinclude EigenSolver_compute.out
       */
-    EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
+    EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true)
+    {
+      return computeInternal(matrix, computeEigenvectors, false, 0);
+    }
+
+    /** \brief Computes eigendecomposition with specified maximum number of iterations. 
+      * 
+      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
+      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
+      *    eigenvalues are computed; if false, only the eigenvalues are
+      *    computed. 
+      * \param[in]  maxIter  Maximum number of iterations.
+      * \returns    Reference to \c *this
+      *
+      * This method provides the same functionality as compute(const MatrixType&, bool) but it also allows the
+      * user to specify the maximum number of iterations to be used when computing the Schur decomposition.
+      */
+    EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors, Index maxIter)
+    {
+      return computeInternal(matrix, computeEigenvectors, true, maxIter);
+    }
 
     ComputationInfo info() const
     {
@@ -283,6 +303,8 @@ template<typename _MatrixType> class EigenSolver
 
   private:
     void doComputeEigenvectors();
+    EigenSolver& computeInternal(const MatrixType& matrix, bool computeEigenvectors, 
+                                 bool maxIterSpecified, Index maxIter);
 
   protected:
     MatrixType m_eivec;
@@ -348,12 +370,18 @@ typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eige
 }
 
 template<typename MatrixType>
-EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
+EigenSolver<MatrixType>& 
+EigenSolver<MatrixType>::computeInternal(const MatrixType& matrix, bool computeEigenvectors, 
+                                         bool maxIterSpecified, Index maxIter)
 {
   assert(matrix.cols() == matrix.rows());
 
   // Reduce to real Schur form.
-  m_realSchur.compute(matrix, computeEigenvectors);
+  if (maxIterSpecified)
+    m_realSchur.compute(matrix, computeEigenvectors, maxIter);
+  else
+    m_realSchur.compute(matrix, computeEigenvectors);
+
   if (m_realSchur.info() == Success)
   {
     m_matT = m_realSchur.matrixT();

Eigen/src/Eigenvalues/RealSchur.h

@@ -2,7 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -160,8 +160,27 @@ template<typename _MatrixType> class RealSchur
       *
       * Example: \include RealSchur_compute.cpp
       * Output: \verbinclude RealSchur_compute.out
+      *
+      * \sa compute(const MatrixType&, bool, Index)
       */
-    RealSchur& compute(const MatrixType& matrix, bool computeU = true);
+    RealSchur& compute(const MatrixType& matrix, bool computeU = true)
+    {
+      return compute(matrix, computeU, m_maxIterations * matrix.rows());
+    }
+
+    /** \brief Computes Schur decomposition with specified maximum number of iterations.
+      *
+      * \param[in]  matrix    Square matrix whose Schur decomposition is to be computed.
+      * \param[in]  computeU  If true, both T and U are computed; if false, only T is computed.
+      * \param[in]  maxIter   Maximum number of iterations. 
+      *
+      * \returns    Reference to \c *this
+      *
+      * This method provides the same functionality as compute(const MatrixType&, bool) but it also allows the
+      * user to specify the maximum number of QR iterations to be used. The maximum number of iterations for
+      * compute(const MatrixType&, bool) is #m_maxIterations times the size of the matrix.
+      */
+    RealSchur& compute(const MatrixType& matrix, bool computeU, Index maxIter);
 
     /** \brief Reports whether previous computation was successful.
       *
@@ -175,7 +194,8 @@ template<typename _MatrixType> class RealSchur
 
     /** \brief Maximum number of iterations.
       *
-      * Maximum number of iterations allowed for an eigenvalue to converge. 
+      * If not otherwise specified, the maximum number of iterations is this number times the size of the
+      * matrix. It is currently set to 40.
       */
     static const int m_maxIterations = 40;
 
@@ -201,7 +221,7 @@ template<typename _MatrixType> class RealSchur
 
 
 template<typename MatrixType>
-RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
+RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU, Index maxIter)
 {
   assert(matrix.cols() == matrix.rows());
 
@@ -253,14 +273,14 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix,
         computeShift(iu, iter, exshift, shiftInfo);
         iter = iter + 1;
         totalIter = totalIter + 1;
-        if (totalIter > m_maxIterations * matrix.cols()) break;
+        if (totalIter > maxIter) break;
         Index im;
         initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
         performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
       }
     }
   }
-  if(totalIter <= m_maxIterations * matrix.cols()) 
+  if(totalIter <= maxIter)
     m_info = Success;
   else
     m_info = NoConvergence;

test/eigensolver_complex.cpp

@@ -59,6 +59,16 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
   // another algorithm so results may differ slightly
   verify_is_approx_upto_permutation(a.eigenvalues(), ei1.eigenvalues());
 
+  ComplexEigenSolver<MatrixType> ei2;
+  ei2.compute(a, true, ComplexSchur<MatrixType>::m_maxIterations * rows);
+  VERIFY_IS_EQUAL(ei2.info(), Success);
+  VERIFY_IS_EQUAL(ei2.eigenvectors(), ei1.eigenvectors());
+  VERIFY_IS_EQUAL(ei2.eigenvalues(), ei1.eigenvalues());
+  if (rows > 2) {
+    ei2.compute(a, true, 1);
+    VERIFY_IS_EQUAL(ei2.info(), NoConvergence);
+  }
+
   ComplexEigenSolver<MatrixType> eiNoEivecs(a, false);
   VERIFY_IS_EQUAL(eiNoEivecs.info(), Success);
   VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());

test/eigensolver_generic.cpp

@@ -2,7 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -45,6 +45,16 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
   VERIFY_IS_APPROX(ei1.eigenvectors().colwise().norm(), RealVectorType::Ones(rows).transpose());
   VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues());
 
+  EigenSolver<MatrixType> ei2;
+  ei2.compute(a, true, RealSchur<MatrixType>::m_maxIterations * rows);
+  VERIFY_IS_EQUAL(ei2.info(), Success);
+  VERIFY_IS_EQUAL(ei2.eigenvectors(), ei1.eigenvectors());
+  VERIFY_IS_EQUAL(ei2.eigenvalues(), ei1.eigenvalues());
+  if (rows > 2) {
+    ei2.compute(a, true, 1);
+    VERIFY_IS_EQUAL(ei2.info(), NoConvergence);
+  }
+
   EigenSolver<MatrixType> eiNoEivecs(a, false);
   VERIFY_IS_EQUAL(eiNoEivecs.info(), Success);
   VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());

test/schur_complex.cpp

@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -47,6 +47,22 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
   VERIFY_IS_EQUAL(cs1.matrixT(), cs2.matrixT());
   VERIFY_IS_EQUAL(cs1.matrixU(), cs2.matrixU());
 
+  // Test maximum number of iterations
+  ComplexSchur<MatrixType> cs3;
+  cs3.compute(A, true, ComplexSchur<MatrixType>::m_maxIterations * size);
+  VERIFY_IS_EQUAL(cs3.info(), Success);
+  VERIFY_IS_EQUAL(cs3.matrixT(), cs1.matrixT());
+  VERIFY_IS_EQUAL(cs3.matrixU(), cs1.matrixU());
+  cs3.compute(A, true, 1);
+  VERIFY_IS_EQUAL(cs3.info(), size > 1 ? NoConvergence : Success);
+
+  MatrixType Atriangular = A;
+  Atriangular.template triangularView<StrictlyLower>().setZero(); 
+  cs3.compute(Atriangular, true, 1); // triangular matrices do not need any iterations
+  VERIFY_IS_EQUAL(cs3.info(), Success);
+  VERIFY_IS_EQUAL(cs3.matrixT(), Atriangular.template cast<ComplexScalar>());
+  VERIFY_IS_EQUAL(cs3.matrixU(), ComplexMatrixType::Identity(size, size));
+
   // Test computation of only T, not U
   ComplexSchur<MatrixType> csOnlyT(A, false);
   VERIFY_IS_EQUAL(csOnlyT.info(), Success);

test/schur_real.cpp

@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -66,6 +66,24 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
   VERIFY_IS_EQUAL(rs1.matrixT(), rs2.matrixT());
   VERIFY_IS_EQUAL(rs1.matrixU(), rs2.matrixU());
 
+  // Test maximum number of iterations
+  RealSchur<MatrixType> rs3;
+  rs3.compute(A, true, RealSchur<MatrixType>::m_maxIterations * size);
+  VERIFY_IS_EQUAL(rs3.info(), Success);
+  VERIFY_IS_EQUAL(rs3.matrixT(), rs1.matrixT());
+  VERIFY_IS_EQUAL(rs3.matrixU(), rs1.matrixU());
+  if (size > 2) {
+    rs3.compute(A, true, 1);
+    VERIFY_IS_EQUAL(rs3.info(), NoConvergence);
+  }
+
+  MatrixType Atriangular = A;
+  Atriangular.template triangularView<StrictlyLower>().setZero(); 
+  rs3.compute(Atriangular, true, 1); // triangular matrices do not need any iterations
+  VERIFY_IS_EQUAL(rs3.info(), Success);
+  VERIFY_IS_EQUAL(rs3.matrixT(), Atriangular);
+  VERIFY_IS_EQUAL(rs3.matrixU(), MatrixType::Identity(size, size));
+
   // Test computation of only T, not U
   RealSchur<MatrixType> rsOnlyT(A, false);
   VERIFY_IS_EQUAL(rsOnlyT.info(), Success);