diff options
author | 2011-11-02 14:18:20 +0000 | |
---|---|---|
committer | 2011-11-02 14:18:20 +0000 | |
commit | a594ac3966541f46e88b785c56d82bfaf71fd302 (patch) | |
tree | 8e226b5ced4a50200656a83324ec935188d98036 /Eigen/src | |
parent | 57207239f330c2716442557a341f3a34e696b9ad (diff) |
Allow for more iterations in SelfAdjointEigenSolver (bug #354).
Add an assert to guard against using eigenvalues that have not converged.
Add call to info() in tutorial example to cover non-convergence.
Diffstat (limited to 'Eigen/src')
-rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 16 |
1 files changed, 10 insertions, 6 deletions
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 84be9bce1..e748f03b0 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -242,6 +242,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver const MatrixType& eigenvectors() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); + eigen_assert(info() == Success && "Eigenvalue computation did not converge."); eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); return m_eivec; } @@ -264,6 +265,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver const RealVectorType& eigenvalues() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); + eigen_assert(info() == Success && "Eigenvalue computation did not converge."); return m_eivalues; } @@ -288,6 +290,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver MatrixType operatorSqrt() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); + eigen_assert(info() == Success && "Eigenvalue computation did not converge."); eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint(); } @@ -313,6 +316,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver MatrixType operatorInverseSqrt() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); + eigen_assert(info() == Success && "Eigenvalue computation did not converge."); eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint(); } @@ -329,7 +333,8 @@ template<typename _MatrixType> class SelfAdjointEigenSolver /** \brief Maximum number of iterations. * - * Maximum number of iterations allowed for an eigenvalue to converge. + * The algorithm terminates if it does not converge within m_maxIterations * n iterations, where n + * denotes the size of the matrix. This value is currently set to 30 (copied from LAPACK). */ static const int m_maxIterations = 30; @@ -429,7 +434,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> Index end = n-1; Index start = 0; - Index iter = 0; // number of iterations we are working on one element + Index iter = 0; // total number of iterations while (end>0) { @@ -440,15 +445,14 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> // 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 + // if we spent too many iterations, we give up iter++; - if(iter > m_maxIterations) break; + if(iter > m_maxIterations * n) break; start = end - 1; while (start>0 && m_subdiag[start-1]!=0) @@ -457,7 +461,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n); } - if (iter <= m_maxIterations) + if (iter <= m_maxIterations * n) m_info = Success; else m_info = NoConvergence; |