diff options
-rw-r--r-- | Eigen/src/Eigenvalues/ComplexEigenSolver.h | 46 | ||||
-rw-r--r-- | Eigen/src/Eigenvalues/ComplexSchur.h | 71 | ||||
-rw-r--r-- | Eigen/src/Eigenvalues/EigenSolver.h | 41 | ||||
-rw-r--r-- | Eigen/src/Eigenvalues/RealSchur.h | 56 | ||||
-rw-r--r-- | test/eigensolver_complex.cpp | 5 | ||||
-rw-r--r-- | test/eigensolver_generic.cpp | 5 | ||||
-rw-r--r-- | test/schur_complex.cpp | 7 | ||||
-rw-r--r-- | test/schur_real.cpp | 7 |
8 files changed, 114 insertions, 124 deletions
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h index 39dfe2a0a..95c70aecb 100644 --- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h +++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h @@ -208,27 +208,7 @@ template<typename _MatrixType> class ComplexEigenSolver * Example: \include ComplexEigenSolver_compute.cpp * Output: \verbinclude ComplexEigenSolver_compute.out */ - 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); - } + ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); /** \brief Reports whether previous computation was successful. * @@ -240,6 +220,19 @@ template<typename _MatrixType> class ComplexEigenSolver return m_schur.info(); } + /** \brief Sets the maximum number of iterations allowed. */ + ComplexEigenSolver& setMaxIterations(Index maxIters) + { + m_schur.setMaxIterations(maxIters); + return *this; + } + + /** \brief Returns the maximum number of iterations. */ + Index getMaxIterations() + { + return m_schur.getMaxIterations(); + } + protected: EigenvectorType m_eivec; EigenvalueType m_eivalues; @@ -251,26 +244,19 @@ 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>::computeInternal(const MatrixType& matrix, bool computeEigenvectors, - bool maxIterSpecified, Index maxIter) +ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) { // 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. - if (maxIterSpecified) - m_schur.compute(matrix, computeEigenvectors, maxIter); - else - m_schur.compute(matrix, computeEigenvectors); + 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 710a8a276..62cbbb14f 100644 --- a/Eigen/src/Eigenvalues/ComplexSchur.h +++ b/Eigen/src/Eigenvalues/ComplexSchur.h @@ -96,7 +96,8 @@ template<typename _MatrixType> class ComplexSchur m_matU(size,size), m_hess(size), m_isInitialized(false), - m_matUisUptodate(false) + m_matUisUptodate(false), + m_maxIters(-1) {} /** \brief Constructor; computes Schur decomposition of given matrix. @@ -109,11 +110,12 @@ template<typename _MatrixType> class ComplexSchur * \sa matrixT() and matrixU() for examples. */ ComplexSchur(const MatrixType& matrix, bool computeU = true) - : m_matT(matrix.rows(),matrix.cols()), - m_matU(matrix.rows(),matrix.cols()), - m_hess(matrix.rows()), - m_isInitialized(false), - m_matUisUptodate(false) + : m_matT(matrix.rows(),matrix.cols()), + m_matU(matrix.rows(),matrix.cols()), + m_hess(matrix.rows()), + m_isInitialized(false), + m_matUisUptodate(false), + m_maxIters(-1) { compute(matrix, computeU); } @@ -184,24 +186,7 @@ template<typename _MatrixType> class ComplexSchur * * \sa compute(const MatrixType&, bool, Index) */ - 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); + ComplexSchur& compute(const MatrixType& matrix, bool computeU = true); /** \brief Reports whether previous computation was successful. * @@ -213,12 +198,29 @@ template<typename _MatrixType> class ComplexSchur return m_info; } - /** \brief Maximum number of iterations. + /** \brief Sets the maximum number of iterations allowed. + * + * If not specified by the user, the maximum number of iterations is m_maxIterationsPerRow times the size + * of the matrix. + */ + ComplexSchur& setMaxIterations(Index maxIters) + { + m_maxIters = maxIters; + return *this; + } + + /** \brief Returns the maximum number of iterations. */ + Index getMaxIterations() + { + return m_maxIters; + } + + /** \brief Maximum number of iterations per row. * * 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; + static const int m_maxIterationsPerRow = 30; protected: ComplexMatrixType m_matT, m_matU; @@ -226,11 +228,12 @@ template<typename _MatrixType> class ComplexSchur ComputationInfo m_info; bool m_isInitialized; bool m_matUisUptodate; + Index m_maxIters; private: bool subdiagonalEntryIsNeglegible(Index i); ComplexScalar computeShift(Index iu, Index iter); - void reduceToTriangularForm(bool computeU, Index maxIter); + void reduceToTriangularForm(bool computeU); friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>; }; @@ -289,7 +292,7 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu template<typename MatrixType> -ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU, Index maxIter) +ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU) { m_matUisUptodate = false; eigen_assert(matrix.cols() == matrix.rows()); @@ -305,7 +308,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, maxIter); + reduceToTriangularForm(computeU); return *this; } @@ -348,8 +351,12 @@ 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, Index maxIter) +void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU) { + Index maxIters = m_maxIters; + if (maxIters == -1) + maxIters = m_maxIterationsPerRow * m_matT.rows(); + // 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. // Rows il,...,iu is the part we are working on (the active submatrix). @@ -375,7 +382,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU, Index maxIt // if we spent too many iterations, we give up iter++; totalIter++; - if(totalIter > maxIter) break; + if(totalIter > maxIters) break; // find il, the top row of the active submatrix il = iu-1; @@ -405,7 +412,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU, Index maxIt } } - if(totalIter <= maxIter) + if(totalIter <= maxIters) m_info = Success; else m_info = NoConvergence; diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index 7f1c94162..9c3bba1e5 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h @@ -273,38 +273,29 @@ 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); + + ComputationInfo info() const { - return computeInternal(matrix, computeEigenvectors, false, 0); + eigen_assert(m_isInitialized && "EigenSolver is not initialized."); + return m_realSchur.info(); } - /** \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) + /** \brief Sets the maximum number of iterations allowed. */ + EigenSolver& setMaxIterations(Index maxIters) { - return computeInternal(matrix, computeEigenvectors, true, maxIter); + m_realSchur.setMaxIterations(maxIters); + return *this; } - ComputationInfo info() const + /** \brief Returns the maximum number of iterations. */ + Index getMaxIterations() { - eigen_assert(m_isInitialized && "EigenSolver is not initialized."); - return m_realSchur.info(); + return m_realSchur.getMaxIterations(); } private: void doComputeEigenvectors(); - EigenSolver& computeInternal(const MatrixType& matrix, bool computeEigenvectors, - bool maxIterSpecified, Index maxIter); protected: MatrixType m_eivec; @@ -371,16 +362,12 @@ typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eige template<typename MatrixType> EigenSolver<MatrixType>& -EigenSolver<MatrixType>::computeInternal(const MatrixType& matrix, bool computeEigenvectors, - bool maxIterSpecified, Index maxIter) +EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) { assert(matrix.cols() == matrix.rows()); // Reduce to real Schur form. - if (maxIterSpecified) - m_realSchur.compute(matrix, computeEigenvectors, maxIter); - else - m_realSchur.compute(matrix, computeEigenvectors); + m_realSchur.compute(matrix, computeEigenvectors); if (m_realSchur.info() == Success) { diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h index 4cdfb17c2..da069bc55 100644 --- a/Eigen/src/Eigenvalues/RealSchur.h +++ b/Eigen/src/Eigenvalues/RealSchur.h @@ -86,7 +86,8 @@ template<typename _MatrixType> class RealSchur m_workspaceVector(size), m_hess(size), m_isInitialized(false), - m_matUisUptodate(false) + m_matUisUptodate(false), + m_maxIters(-1) { } /** \brief Constructor; computes real Schur decomposition of given matrix. @@ -105,7 +106,8 @@ template<typename _MatrixType> class RealSchur m_workspaceVector(matrix.rows()), m_hess(matrix.rows()), m_isInitialized(false), - m_matUisUptodate(false) + m_matUisUptodate(false), + m_maxIters(-1) { compute(matrix, computeU); } @@ -163,24 +165,7 @@ template<typename _MatrixType> class RealSchur * * \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, Index maxIter); + RealSchur& compute(const MatrixType& matrix, bool computeU = true); /** \brief Reports whether previous computation was successful. * @@ -192,12 +177,29 @@ template<typename _MatrixType> class RealSchur return m_info; } - /** \brief Maximum number of iterations. + /** \brief Sets the maximum number of iterations allowed. + * + * If not specified by the user, the maximum number of iterations is m_maxIterationsPerRow times the size + * of the matrix. + */ + RealSchur& setMaxIterations(Index maxIters) + { + m_maxIters = maxIters; + return *this; + } + + /** \brief Returns the maximum number of iterations. */ + Index getMaxIterations() + { + return m_maxIters; + } + + /** \brief Maximum number of iterations per row. * * 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; + static const int m_maxIterationsPerRow = 40; private: @@ -208,6 +210,7 @@ template<typename _MatrixType> class RealSchur ComputationInfo m_info; bool m_isInitialized; bool m_matUisUptodate; + Index m_maxIters; typedef Matrix<Scalar,3,1> Vector3s; @@ -221,9 +224,12 @@ template<typename _MatrixType> class RealSchur template<typename MatrixType> -RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU, Index maxIter) +RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU) { assert(matrix.cols() == matrix.rows()); + Index maxIters = m_maxIters; + if (maxIters == -1) + maxIters = m_maxIterationsPerRow * matrix.rows(); // Step 1. Reduce to Hessenberg form m_hess.compute(matrix); @@ -273,14 +279,14 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, computeShift(iu, iter, exshift, shiftInfo); iter = iter + 1; totalIter = totalIter + 1; - if (totalIter > maxIter) break; + if (totalIter > maxIters) break; Index im; initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector); performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace); } } } - if(totalIter <= maxIter) + if(totalIter <= maxIters) m_info = Success; else m_info = NoConvergence; diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp index a569d8507..aef125739 100644 --- a/test/eigensolver_complex.cpp +++ b/test/eigensolver_complex.cpp @@ -60,13 +60,14 @@ template<typename MatrixType> void eigensolver(const MatrixType& m) verify_is_approx_upto_permutation(a.eigenvalues(), ei1.eigenvalues()); ComplexEigenSolver<MatrixType> ei2; - ei2.compute(a, true, ComplexSchur<MatrixType>::m_maxIterations * rows); + ei2.setMaxIterations(ComplexSchur<MatrixType>::m_maxIterationsPerRow * rows).compute(a); 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); + ei2.setMaxIterations(1).compute(a); VERIFY_IS_EQUAL(ei2.info(), NoConvergence); + VERIFY_IS_EQUAL(ei2.getMaxIterations(), 1); } ComplexEigenSolver<MatrixType> eiNoEivecs(a, false); diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp index 1b83afbc4..ef499a989 100644 --- a/test/eigensolver_generic.cpp +++ b/test/eigensolver_generic.cpp @@ -46,13 +46,14 @@ template<typename MatrixType> void eigensolver(const MatrixType& m) VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues()); EigenSolver<MatrixType> ei2; - ei2.compute(a, true, RealSchur<MatrixType>::m_maxIterations * rows); + ei2.setMaxIterations(RealSchur<MatrixType>::m_maxIterationsPerRow * rows).compute(a); 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); + ei2.setMaxIterations(1).compute(a); VERIFY_IS_EQUAL(ei2.info(), NoConvergence); + VERIFY_IS_EQUAL(ei2.getMaxIterations(), 1); } EigenSolver<MatrixType> eiNoEivecs(a, false); diff --git a/test/schur_complex.cpp b/test/schur_complex.cpp index 4d5cbbfdf..5e869790f 100644 --- a/test/schur_complex.cpp +++ b/test/schur_complex.cpp @@ -49,16 +49,17 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim // Test maximum number of iterations ComplexSchur<MatrixType> cs3; - cs3.compute(A, true, ComplexSchur<MatrixType>::m_maxIterations * size); + cs3.setMaxIterations(ComplexSchur<MatrixType>::m_maxIterationsPerRow * size).compute(A); 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); + cs3.setMaxIterations(1).compute(A); VERIFY_IS_EQUAL(cs3.info(), size > 1 ? NoConvergence : Success); + VERIFY_IS_EQUAL(cs3.getMaxIterations(), 1); MatrixType Atriangular = A; Atriangular.template triangularView<StrictlyLower>().setZero(); - cs3.compute(Atriangular, true, 1); // triangular matrices do not need any iterations + cs3.setMaxIterations(1).compute(Atriangular); // 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)); diff --git a/test/schur_real.cpp b/test/schur_real.cpp index e81a28d7d..36b9c24d1 100644 --- a/test/schur_real.cpp +++ b/test/schur_real.cpp @@ -68,18 +68,19 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim // Test maximum number of iterations RealSchur<MatrixType> rs3; - rs3.compute(A, true, RealSchur<MatrixType>::m_maxIterations * size); + rs3.setMaxIterations(RealSchur<MatrixType>::m_maxIterationsPerRow * size).compute(A); 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); + rs3.setMaxIterations(1).compute(A); VERIFY_IS_EQUAL(rs3.info(), NoConvergence); + VERIFY_IS_EQUAL(rs3.getMaxIterations(), 1); } MatrixType Atriangular = A; Atriangular.template triangularView<StrictlyLower>().setZero(); - rs3.compute(Atriangular, true, 1); // triangular matrices do not need any iterations + rs3.setMaxIterations(1).compute(Atriangular); // 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)); |