aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/Eigenvalues/RealSchur.h15
-rw-r--r--test/eigensolver_generic.cpp74
2 files changed, 77 insertions, 12 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--;
}
diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp
index e0e435151..086ecdf5e 100644
--- a/test/eigensolver_generic.cpp
+++ b/test/eigensolver_generic.cpp
@@ -12,6 +12,21 @@
#include <limits>
#include <Eigen/Eigenvalues>
+template<typename EigType,typename MatType>
+void check_eigensolver_for_given_mat(const EigType &eig, const MatType& a)
+{
+ typedef typename NumTraits<typename MatType::Scalar>::Real RealScalar;
+ typedef Matrix<RealScalar, MatType::RowsAtCompileTime, 1> RealVectorType;
+ typedef typename std::complex<RealScalar> Complex;
+ Index n = a.rows();
+ VERIFY_IS_EQUAL(eig.info(), Success);
+ VERIFY_IS_APPROX(a * eig.pseudoEigenvectors(), eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix());
+ VERIFY_IS_APPROX(a.template cast<Complex>() * eig.eigenvectors(),
+ eig.eigenvectors() * eig.eigenvalues().asDiagonal());
+ VERIFY_IS_APPROX(eig.eigenvectors().colwise().norm(), RealVectorType::Ones(n).transpose());
+ VERIFY_IS_APPROX(a.eigenvalues(), eig.eigenvalues());
+}
+
template<typename MatrixType> void eigensolver(const MatrixType& m)
{
/* this test covers the following files:
@@ -22,8 +37,7 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
- typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
- typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
+ typedef typename std::complex<RealScalar> Complex;
MatrixType a = MatrixType::Random(rows,cols);
MatrixType a1 = MatrixType::Random(rows,cols);
@@ -36,12 +50,7 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
(ei0.pseudoEigenvectors().template cast<Complex>()) * (ei0.eigenvalues().asDiagonal()));
EigenSolver<MatrixType> ei1(a);
- VERIFY_IS_EQUAL(ei1.info(), Success);
- VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix());
- VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(),
- ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
- VERIFY_IS_APPROX(ei1.eigenvectors().colwise().norm(), RealVectorType::Ones(rows).transpose());
- VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues());
+ CALL_SUBTEST( check_eigensolver_for_given_mat(ei1,a) );
EigenSolver<MatrixType> ei2;
ei2.setMaxIterations(RealSchur<MatrixType>::m_maxIterationsPerRow * rows).compute(a);
@@ -100,6 +109,19 @@ template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m
VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors());
}
+
+template<typename CoeffType>
+Matrix<typename CoeffType::Scalar,Dynamic,Dynamic>
+make_companion(const CoeffType& coeffs)
+{
+ Index n = coeffs.size()-1;
+ Matrix<typename CoeffType::Scalar,Dynamic,Dynamic> res(n,n);
+ res.setZero();
+ res.row(0) = -coeffs.tail(n) / coeffs(0);
+ res.diagonal(-1).setOnes();
+ return res;
+}
+
template<int>
void eigensolver_generic_extra()
{
@@ -126,6 +148,42 @@ void eigensolver_generic_extra()
VERIFY_IS_APPROX((a * eig.eigenvectors()).norm()+1., 1.);
VERIFY_IS_APPROX((eig.eigenvectors() * eig.eigenvalues().asDiagonal()).norm()+1., 1.);
}
+
+ // regression test for bug 933
+ {
+ {
+ VectorXd coeffs(5); coeffs << 1, -3, -175, -225, 2250;
+ MatrixXd C = make_companion(coeffs);
+ EigenSolver<MatrixXd> eig(C);
+ CALL_SUBTEST( check_eigensolver_for_given_mat(eig,C) );
+ }
+ {
+ // this test is tricky because it requires high accuracy in smallest eigenvalues
+ VectorXd coeffs(5); coeffs << 6.154671e-15, -1.003870e-10, -9.819570e-01, 3.995715e+03, 2.211511e+08;
+ MatrixXd C = make_companion(coeffs);
+ EigenSolver<MatrixXd> eig(C);
+ CALL_SUBTEST( check_eigensolver_for_given_mat(eig,C) );
+ Index n = C.rows();
+ for(Index i=0;i<n;++i)
+ {
+ typedef std::complex<double> Complex;
+ MatrixXcd ac = C.cast<Complex>();
+ ac.diagonal().array() -= eig.eigenvalues()(i);
+ VectorXd sv = ac.jacobiSvd().singularValues();
+ // comparing to sv(0) is not enough here to catch the "bug",
+ // the hard-coded 1.0 is important!
+ VERIFY_IS_MUCH_SMALLER_THAN(sv(n-1), 1.0);
+ }
+ }
+ }
+ // regression test for bug 1557
+ {
+ // this test is interesting because it contains zeros on the diagonal.
+ MatrixXd A_bug1557(3,3);
+ A_bug1557 << 0, 0, 0, 1, 0, 0.5887907064808635127, 0, 1, 0;
+ EigenSolver<MatrixXd> eig(A_bug1557);
+ CALL_SUBTEST( check_eigensolver_for_given_mat(eig,A_bug1557) );
+ }
}
EIGEN_DECLARE_TEST(eigensolver_generic)
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-04 15:27:57 +0000