devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix real schur and polynomial solver.

Browse Source
This commit is contained in:
Antonio Sánchez
2024-02-17 15:22:11 +00:00
parent 8a4118746e
commit b14c5d0fa1
3 changed files with 33 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -408,28 +408,29 @@ inline void RealSchur<MatrixType>::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));
}
}
}
}