diff options
author | Jitse Niesen <jitse@maths.leeds.ac.uk> | 2012-07-24 15:17:59 +0100 |
---|---|---|
committer | Jitse Niesen <jitse@maths.leeds.ac.uk> | 2012-07-24 15:17:59 +0100 |
commit | ba5eecae53aa038374d1708573cf03a2df3f76f3 (patch) | |
tree | b1bb109190fc9b31cbbda8c227343170e2058307 /Eigen | |
parent | b7ac053b9c1196bf775ca0ec45765f6262648c43 (diff) |
Allow user to specify max number of iterations (bug #479).
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/src/Eigenvalues/ComplexEigenSolver.h | 36 | ||||
-rw-r--r-- | Eigen/src/Eigenvalues/ComplexSchur.h | 41 | ||||
-rw-r--r-- | Eigen/src/Eigenvalues/EigenSolver.h | 36 | ||||
-rw-r--r-- | Eigen/src/Eigenvalues/RealSchur.h | 32 |
4 files changed, 121 insertions, 24 deletions
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h index c4b8a308c..39dfe2a0a 100644 --- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h +++ b/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) { diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h index 55aeedb90..710a8a276 100644 --- a/Eigen/src/Eigenvalues/ComplexSchur.h +++ b/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; diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index c16ff2b74..315958421 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/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(); diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h index d1949b83c..da877dd04 100644 --- a/Eigen/src/Eigenvalues/RealSchur.h +++ b/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) + { + 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 = true); + 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; |