mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-02-17 18:09:55 +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 57207239f3
commit a594ac3966
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

@ -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;

3
doc/examples/TutorialLinAlgSelfAdjointEigenSolver.cpp
View File

@ -10,8 +10,9 @@ int main()
A << 1, 2, 2, 3;
cout << "Here is the matrix A:\n" << A << endl;
SelfAdjointEigenSolver<Matrix2f> eigensolver(A);
if (eigensolver.info() != Success) abort();
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"
<< eigensolver.eigenvectors() << endl;
}