diff options
author | Gael Guennebaud <g.gael@free.fr> | 2018-12-12 17:30:08 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2018-12-12 17:30:08 +0100 |
commit | 2de8da70fd0b35849845dc76b2741bb0689f0643 (patch) | |
tree | 5f322f6dde3e97799c62c4f766279967b769516b /Eigen/src/Eigenvalues | |
parent | 72c0bbe2bd1c49c75b6efdb81d0558f8b62578d1 (diff) |
bug #1557: fix RealSchur and EigenSolver for matrices with only zeros on the diagonal.
Diffstat (limited to 'Eigen/src/Eigenvalues')
-rw-r--r-- | Eigen/src/Eigenvalues/RealSchur.h | 15 |
1 files changed, 11 insertions, 4 deletions
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h index aca8a8279..8dbd9e314 100644 --- a/Eigen/src/Eigenvalues/RealSchur.h +++ b/Eigen/src/Eigenvalues/RealSchur.h @@ -236,7 +236,7 @@ template<typename _MatrixType> class RealSchur typedef Matrix<Scalar,3,1> Vector3s; Scalar computeNormOfT(); - Index findSmallSubdiagEntry(Index iu); + Index findSmallSubdiagEntry(Index iu, const Scalar& considerAsZero); void splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift); void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo); void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector); @@ -307,12 +307,16 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMa Index totalIter = 0; // iteration count for whole matrix Scalar exshift(0); // sum of exceptional shifts Scalar norm = computeNormOfT(); + // sub-diagonal entries smaller than considerAsZero will be treated as zero. + // We use eps^2 to enable more precision in small eigenvalues. + Scalar considerAsZero = numext::maxi( norm * numext::abs2(NumTraits<Scalar>::epsilon()), + (std::numeric_limits<Scalar>::min)() ); if(norm!=Scalar(0)) { while (iu >= 0) { - Index il = findSmallSubdiagEntry(iu); + Index il = findSmallSubdiagEntry(iu,considerAsZero); // Check for convergence if (il == iu) // One root found @@ -369,14 +373,17 @@ inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT() /** \internal Look for single small sub-diagonal element and returns its index */ template<typename MatrixType> -inline Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu) +inline Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, const Scalar& considerAsZero) { using std::abs; Index res = iu; while (res > 0) { 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) + + s = numext::maxi(s * NumTraits<Scalar>::epsilon(), considerAsZero); + + if (abs(m_matT.coeff(res,res-1)) <= s) break; res--; } |