GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-07 21:51:47 +08:00
Code Issues Packages Projects Releases Wiki Activity

Allow for more iterations in SelfAdjointEigenSolver (bug #354).

Browse Source 
Add an assert to guard against using eigenvalues that have not converged.
Add call to info() in tutorial example to cover non-convergence.
This commit is contained in:
Jitse Niesen 2011-11-02 14:18:20 +00:00
parent 411b4a1b1d
commit 1ab1f7b125
2 changed files with 12 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@ -220,6 +220,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
const MatrixType& eigenvectors() const const MatrixType& eigenvectors() const
{ {
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); 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."); eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
return m_eivec; return m_eivec;
} }
@ -242,6 +243,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
const RealVectorType& eigenvalues() const const RealVectorType& eigenvalues() const
{ {
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
eigen_assert(info() == Success && "Eigenvalue computation did not converge.");
return m_eivalues; return m_eivalues;
} }
@ -266,6 +268,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
MatrixType operatorSqrt() const MatrixType operatorSqrt() const
{ {
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); 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."); eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint(); return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
} }
@ -291,6 +294,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
MatrixType operatorInverseSqrt() const MatrixType operatorInverseSqrt() const
{ {
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); 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."); 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(); return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
} }
@ -307,7 +311,8 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
/** \brief Maximum number of iterations. /** \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; static const int m_maxIterations = 30;
@ -407,7 +412,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
Index end = n-1; Index end = n-1;
Index start = 0; 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) while (end>0)
{ {
@ -418,15 +423,14 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
// find the largest unreduced block // find the largest unreduced block
while (end>0 && m_subdiag[end-1]==0) while (end>0 && m_subdiag[end-1]==0)
{ {
iter = 0;
end--; end--;
} }
if (end<=0) if (end<=0)
break; break;
// if we spent too many iterations on the current element, we give up // if we spent too many iterations, we give up
iter++; iter++;
if(iter > m_maxIterations) break; if(iter > m_maxIterations * n) break;
start = end - 1; start = end - 1;
while (start>0 && m_subdiag[start-1]!=0) while (start>0 && m_subdiag[start-1]!=0)
@ -435,7 +439,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); 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; m_info = Success;
else else
m_info = NoConvergence; m_info = NoConvergence;

3
doc/examples/TutorialLinAlgSelfAdjointEigenSolver.cpp
View File

@ -10,8 +10,9 @@ int main()
A << 1, 2, 2, 3; A << 1, 2, 2, 3;
cout << "Here is the matrix A:\n" << A << endl; cout << "Here is the matrix A:\n" << A << endl;
SelfAdjointEigenSolver<Matrix2f> eigensolver(A); SelfAdjointEigenSolver<Matrix2f> eigensolver(A);
if (eigensolver.info() != Success) abort();
cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << endl; cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << endl;
cout << "Here's a matrix whose columns are eigenvectors of A " cout << "Here's a matrix whose columns are eigenvectors of A \n"
<< "corresponding to these eigenvalues:\n" << "corresponding to these eigenvalues:\n"
<< eigensolver.eigenvectors() << endl; << eigensolver.eigenvectors() << endl;
} }