Add field m_maxIterations; break loop when this limit is exceeded.

This commit is contained in:
Jitse Niesen 2010-06-02 17:31:02 +01:00
parent 9ff0d67156
commit 38d8352b7b
3 changed files with 42 additions and 5 deletions

View File

@ -192,6 +192,12 @@ template<typename _MatrixType> 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<MatrixType> m_hess;
@ -374,9 +380,9 @@ void ComplexSchur<MatrixType>::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<MatrixType>::reduceToTriangularForm(bool computeU)
}
}
if(iter >= 30)
if(iter >= m_maxIterations)
{
// FIXME : what to do when iter==MAXITER ??
// std::cerr << "MAXITER" << std::endl;

View File

@ -176,6 +176,12 @@ template<typename _MatrixType> 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<MatrixType>& RealSchur<MatrixType>::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;
}

View File

@ -360,6 +360,11 @@ template<typename _MatrixType> 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<MatrixType>& SelfAdjointEigenSolver<MatrixType>::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<end; ++i)
@ -427,9 +434,17 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::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<MatrixType>& SelfAdjointEigenSolver<MatrixType>::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 !!