Use EISPACK's strategy re max number of iters in Schur decomposition (bug #479).

This commit is contained in:
Jitse Niesen 2012-07-22 22:02:50 +01:00
parent fd5749f51c
commit b7ac053b9c
2 changed files with 11 additions and 7 deletions

View File

@ -336,6 +336,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
Index iu = m_matT.cols() - 1;
Index il;
Index iter = 0; // number of iterations we are working on the (iu,iu) element
Index totalIter = 0; // number of iterations for whole matrix
while(true)
{
@ -350,9 +351,10 @@ 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 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;
totalIter++;
if(totalIter > m_maxIterations * m_matT.cols()) break;
// find il, the top row of the active submatrix
il = iu-1;
@ -382,7 +384,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
}
}
if(iter <= m_maxIterations)
if(totalIter <= m_maxIterations * m_matT.cols())
m_info = Success;
else
m_info = NoConvergence;

View File

@ -220,8 +220,9 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix,
// Rows il,...,iu is the part we are working on (the active window).
// Rows iu+1,...,end are already brought in triangular form.
Index iu = m_matT.cols() - 1;
Index iter = 0; // iteration count
Scalar exshift(0); // sum of exceptional shifts
Index iter = 0; // iteration count for current eigenvalue
Index totalIter = 0; // iteration count for whole matrix
Scalar exshift(0); // sum of exceptional shifts
Scalar norm = computeNormOfT();
if(norm!=0)
@ -251,14 +252,15 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix,
Vector3s firstHouseholderVector(0,0,0), shiftInfo;
computeShift(iu, iter, exshift, shiftInfo);
iter = iter + 1;
if (iter > m_maxIterations) break;
totalIter = totalIter + 1;
if (totalIter > m_maxIterations * matrix.cols()) break;
Index im;
initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
}
}
}
if(iter <= m_maxIterations)
if(totalIter <= m_maxIterations * matrix.cols())
m_info = Success;
else
m_info = NoConvergence;