diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h index a4812dac6..cf6755bfc 100644 --- a/Eigen/src/Eigenvalues/ComplexSchur.h +++ b/Eigen/src/Eigenvalues/ComplexSchur.h @@ -192,6 +192,12 @@ template class ComplexSchur */ ComplexSchur& compute(const MatrixType& matrix, bool computeU = true); + /** \brief Maximum number of iterations. + * + * Maximum number of iterations allowed for an eigenvalue to converge. + */ + static const int m_maxIterations = 30; + protected: ComplexMatrixType m_matT, m_matU; HessenbergDecomposition m_hess; @@ -374,9 +380,9 @@ void ComplexSchur::reduceToTriangularForm(bool computeU) // if iu is zero then we are done; the whole matrix is triangularized if(iu==0) break; - // if we spent 30 iterations on the current element, we give up + // if we spent too many iterations on the current element, we give up iter++; - if(iter >= 30) break; + if(iter >= m_maxIterations) break; // find il, the top row of the active submatrix il = iu-1; @@ -406,7 +412,7 @@ void ComplexSchur::reduceToTriangularForm(bool computeU) } } - if(iter >= 30) + if(iter >= m_maxIterations) { // FIXME : what to do when iter==MAXITER ?? // std::cerr << "MAXITER" << std::endl; diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h index 41f74e530..584b1d05e 100644 --- a/Eigen/src/Eigenvalues/RealSchur.h +++ b/Eigen/src/Eigenvalues/RealSchur.h @@ -176,6 +176,12 @@ template class RealSchur */ RealSchur& compute(const MatrixType& matrix, bool computeU = true); + /** \brief Maximum number of iterations. + * + * Maximum number of iterations allowed for an eigenvalue to converge. + */ + static const int m_maxIterations = 40; + private: MatrixType m_matT; @@ -244,14 +250,19 @@ RealSchur& RealSchur::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; Index im; initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector); performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace); } } - m_isInitialized = true; - m_matUisUptodate = computeU; + if(iter < m_maxIterations) + { + m_isInitialized = true; + m_matUisUptodate = computeU; + } + return *this; } diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 20773a549..b08bbd7e2 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -360,6 +360,11 @@ template class SelfAdjointEigenSolver return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint(); } + /** \brief Maximum number of iterations. + * + * Maximum number of iterations allowed for an eigenvalue to converge. + */ + static const int m_maxIterations = 30; protected: MatrixType m_eivec; @@ -419,6 +424,8 @@ SelfAdjointEigenSolver& SelfAdjointEigenSolver::compute( Index end = n-1; Index start = 0; + Index iter = 0; // number of iterations we are working on one element + while (end>0) { for (Index i = start; i& SelfAdjointEigenSolver::compute( // find the largest unreduced block while (end>0 && m_subdiag[end-1]==0) + { + iter = 0; end--; + } if (end<=0) break; + + // if we spent too many iterations on the current element, we give up + iter++; + if(iter >= m_maxIterations) break; + start = end - 1; while (start>0 && m_subdiag[start-1]!=0) start--; @@ -437,6 +452,11 @@ SelfAdjointEigenSolver& SelfAdjointEigenSolver::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; + } + // Sort eigenvalues and corresponding vectors. // TODO make the sort optional ? // TODO use a better sort algorithm !!