mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-24 14:45:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

Revert part of changeset 5b3a6f51d3

Browse Source 
to keep accuracy of smallest eigenvalues.
This commit is contained in:
Gael Guennebaud 2016-08-29 14:14:18 +02:00
parent 7e029d1d6e
commit 2915e1fc5d
1 changed files with 5 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
Eigen/src/Eigenvalues/RealSchur.h
View File

@ -236,7 +236,7 @@ template<typename _MatrixType> class RealSchur
typedef Matrix<Scalar,3,1> Vector3s; typedef Matrix<Scalar,3,1> Vector3s;
Scalar computeNormOfT(); Scalar computeNormOfT();
Index findSmallSubdiagEntry(Index iu, const Scalar& maxDiagEntry); Index findSmallSubdiagEntry(Index iu);
void splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift); void splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift);
void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo); void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo);
void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector); void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector);
@ -293,18 +293,14 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMa
if(norm!=0) if(norm!=0)
{ {
Scalar maxDiagEntry = m_matT.cwiseAbs().diagonal().maxCoeff();
while (iu >= 0) while (iu >= 0)
{ {
Index il = findSmallSubdiagEntry(iu,maxDiagEntry); Index il = findSmallSubdiagEntry(iu);
// Check for convergence // Check for convergence
if (il == iu) // One root found if (il == iu) // One root found
{ {
m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift; m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
// keep track of the largest diagonal coefficient
maxDiagEntry = numext::maxi<Scalar>(maxDiagEntry,abs(m_matT.coeffRef(iu,iu)));
if (iu > 0) if (iu > 0)
m_matT.coeffRef(iu, iu-1) = Scalar(0); m_matT.coeffRef(iu, iu-1) = Scalar(0);
iu--; iu--;
@ -313,8 +309,6 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMa
else if (il == iu-1) // Two roots found else if (il == iu-1) // Two roots found
{ {
splitOffTwoRows(iu, computeU, exshift); splitOffTwoRows(iu, computeU, exshift);
// keep track of the largest diagonal coefficient
maxDiagEntry = numext::maxi<Scalar>(maxDiagEntry,numext::maxi(abs(m_matT.coeff(iu,iu)), abs(m_matT.coeff(iu-1,iu-1))));
iu -= 2; iu -= 2;
iter = 0; iter = 0;
} }
@ -329,8 +323,6 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMa
Index im; Index im;
initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector); initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace); performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
// keep track of the largest diagonal coefficient
maxDiagEntry = numext::maxi(maxDiagEntry,m_matT.cwiseAbs().diagonal().segment(im,iu-im).maxCoeff());
} }
} }
} }
@ -360,13 +352,14 @@ inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
/** \internal Look for single small sub-diagonal element and returns its index */ /** \internal Look for single small sub-diagonal element and returns its index */
template<typename MatrixType> template<typename MatrixType>
inline Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, const Scalar& maxDiagEntry) inline Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu)
{ {
using std::abs; using std::abs;
Index res = iu; Index res = iu;
while (res > 0) while (res > 0)
{ {
if (abs(m_matT.coeff(res,res-1)) <= NumTraits<Scalar>::epsilon() * maxDiagEntry) Scalar s = abs(m_matT.coeff(res-1,res-1)) + abs(m_matT.coeff(res,res));
if (abs(m_matT.coeff(res,res-1)) <= NumTraits<Scalar>::epsilon() * s)
break; break;
res--; res--;
} }