Fix real schur and polynomial solver.
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h
index 970500c..5cef658 100644
--- a/Eigen/src/Eigenvalues/RealSchur.h
+++ b/Eigen/src/Eigenvalues/RealSchur.h
@@ -408,28 +408,29 @@
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 cd0be92..4a9dd89 100644
--- a/test/schur_real.cpp
+++ b/test/schur_real.cpp
@@ -97,6 +97,13 @@
}
}
+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<Eigen::MatrixXd> schur(A);
+ VERIFY(schur.info() == Eigen::Success);
+}
+
EIGEN_DECLARE_TEST(schur_real) {
CALL_SUBTEST_1((schur<Matrix4f>()));
CALL_SUBTEST_2((schur<MatrixXd>(internal::random<int>(1, EIGEN_TEST_MAX_SIZE / 4))));
@@ -105,4 +112,6 @@
// Test problem size constructors
CALL_SUBTEST_5(RealSchur<MatrixXf>(10));
+
+ CALL_SUBTEST_6((test_bug2633()));
}
diff --git a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h
index ec5bc54..8c0ce3b 100644
--- a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h
+++ b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h
@@ -320,6 +320,7 @@
internal::companion<Scalar, Deg_> 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
