GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-15 01:21:48 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix real schur and polynomial solver.

Browse Source 
For certain inputs, the real schur decomposition might get stuck in a cycle.
Exceptional shifts are supposed to knock us out of that - but previously
they were only ever applied at iteration 10 and 30, which doesn't help if
the cycle starts after cycle 30.  Modified to apply a shift every 16 iterations
(for reference, LAPACK seems to do it every 6 iterations).

Also added an assert in polynomial solver to verify that the schur decomposition
was successful.

Fixes #2633.
This commit is contained in:
Antonio Sanchez 2024-02-16 13:11:54 -08:00
parent 287c801780
commit 3ee06ec52f
3 changed files with 39 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -435,34 +435,33 @@ inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& ex
shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1); 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); shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
// Wilkinson's original ad hoc shift // Alternate exceptional shifting strategy every 16 iterations.
if (iter == 10) if (iter % 16 == 0) {
{ // Wilkinson's original ad hoc shift
exshift += shiftInfo.coeff(0); if (iter % 32 != 0) {
for (Index i = 0; i <= iu; ++i) exshift += shiftInfo.coeff(0);
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) for (Index i = 0; i <= iu; ++i)
m_matT.coeffRef(i,i) -= s; m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
shiftInfo.setConstant(Scalar(0.964)); 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));
}
} }
} }
} }

12
test/schur_real.cpp
View File

@ -98,6 +98,16 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
} }
} }
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) EIGEN_DECLARE_TEST(schur_real)
{ {
CALL_SUBTEST_1(( schur<Matrix4f>() )); CALL_SUBTEST_1(( schur<Matrix4f>() ));
@ -107,4 +117,6 @@ EIGEN_DECLARE_TEST(schur_real)
// Test problem size constructors // Test problem size constructors
CALL_SUBTEST_5(RealSchur<MatrixXf>(10)); CALL_SUBTEST_5(RealSchur<MatrixXf>(10));
CALL_SUBTEST_6(( test_bug2633() ));
} }

1
unsupported/Eigen/src/Polynomials/PolynomialSolver.h
View File

@ -354,6 +354,7 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
internal::companion<Scalar,_Deg> companion( poly ); internal::companion<Scalar,_Deg> companion( poly );
companion.balance(); companion.balance();
m_eigenSolver.compute( companion.denseMatrix() ); m_eigenSolver.compute( companion.denseMatrix() );
eigen_assert(m_eigenSolver.info() == Eigen::Success);
m_roots = m_eigenSolver.eigenvalues(); m_roots = m_eigenSolver.eigenvalues();
// cleanup noise in imaginary part of real roots: // cleanup noise in imaginary part of real roots:
// if the imaginary part is rather small compared to the real part // if the imaginary part is rather small compared to the real part