diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h index 970500cd3..5cef6587b 100644 --- a/Eigen/src/Eigenvalues/RealSchur.h +++ b/Eigen/src/Eigenvalues/RealSchur.h @@ -408,28 +408,29 @@ inline void RealSchur::computeShift(Index iu, Index iter, Scalar& ex shiftInfo.coeffRef(1) = m_matT.coeff(iu - 1, iu - 1); shiftInfo.coeffRef(2) = m_matT.coeff(iu, iu - 1) * m_matT.coeff(iu - 1, iu); - // Wilkinson's original ad hoc shift - if (iter == 10) { - exshift += shiftInfo.coeff(0); - for (Index i = 0; i <= iu; ++i) m_matT.coeffRef(i, i) -= shiftInfo.coeff(0); - Scalar s = abs(m_matT.coeff(iu, iu - 1)) + abs(m_matT.coeff(iu - 1, iu - 2)); - shiftInfo.coeffRef(0) = Scalar(0.75) * s; - shiftInfo.coeffRef(1) = Scalar(0.75) * s; - shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s; - } - - // MATLAB's new ad hoc shift - if (iter == 30) { - Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0); - s = s * s + shiftInfo.coeff(2); - if (s > Scalar(0)) { - s = sqrt(s); - if (shiftInfo.coeff(1) < shiftInfo.coeff(0)) s = -s; - s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0); - s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s; - exshift += s; - for (Index i = 0; i <= iu; ++i) m_matT.coeffRef(i, i) -= s; - shiftInfo.setConstant(Scalar(0.964)); + // Alternate exceptional shifting strategy every 16 iterations. + if (iter % 16 == 0) { + // Wilkinson's original ad hoc shift + if (iter % 32 != 0) { + exshift += shiftInfo.coeff(0); + for (Index i = 0; i <= iu; ++i) m_matT.coeffRef(i, i) -= shiftInfo.coeff(0); + Scalar s = abs(m_matT.coeff(iu, iu - 1)) + abs(m_matT.coeff(iu - 1, iu - 2)); + shiftInfo.coeffRef(0) = Scalar(0.75) * s; + shiftInfo.coeffRef(1) = Scalar(0.75) * s; + shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s; + } else { + // MATLAB's new ad hoc shift + Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0); + s = s * s + shiftInfo.coeff(2); + if (s > Scalar(0)) { + s = sqrt(s); + if (shiftInfo.coeff(1) < shiftInfo.coeff(0)) s = -s; + s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0); + s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s; + exshift += s; + for (Index i = 0; i <= iu; ++i) m_matT.coeffRef(i, i) -= s; + shiftInfo.setConstant(Scalar(0.964)); + } } } } diff --git a/test/schur_real.cpp b/test/schur_real.cpp index cd0be9251..4a9dd891f 100644 --- a/test/schur_real.cpp +++ b/test/schur_real.cpp @@ -97,6 +97,13 @@ void schur(int size = MatrixType::ColsAtCompileTime) { } } +void test_bug2633() { + Eigen::MatrixXd A(4, 4); + A << 0, 0, 0, -2, 1, 0, 0, -0, 0, 1, 0, 2, 0, 0, 2, -0; + RealSchur schur(A); + VERIFY(schur.info() == Eigen::Success); +} + EIGEN_DECLARE_TEST(schur_real) { CALL_SUBTEST_1((schur())); CALL_SUBTEST_2((schur(internal::random(1, EIGEN_TEST_MAX_SIZE / 4)))); @@ -105,4 +112,6 @@ EIGEN_DECLARE_TEST(schur_real) { // Test problem size constructors CALL_SUBTEST_5(RealSchur(10)); + + CALL_SUBTEST_6((test_bug2633())); } diff --git a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h index ec5bc54b2..8c0ce3b71 100644 --- a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h +++ b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h @@ -320,6 +320,7 @@ class PolynomialSolver : public PolynomialSolverBase { internal::companion companion(poly); companion.balance(); m_eigenSolver.compute(companion.denseMatrix()); + eigen_assert(m_eigenSolver.info() == Eigen::Success); m_roots = m_eigenSolver.eigenvalues(); // cleanup noise in imaginary part of real roots: // if the imaginary part is rather small compared to the real part