GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-13 04:09:10 +08:00
Code Issues Packages Projects Releases Wiki Activity

RealSchur: Make sure zeros are really zero; simplify initFrancisQRStep().

Browse Source
This commit is contained in:
Jitse Niesen 2010-04-09 11:23:17 +01:00
parent 7dea3a33a5
commit 872df22ca4
1 changed files with 14 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -137,6 +137,7 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix)
if (il == iu) // One root found if (il == iu) // One root found
{ {
m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift; m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
m_matT.block(iu, std::max(0,iu-2), 1, iu - std::max(0,iu-2)).setZero();
m_eivalues.coeffRef(iu) = ComplexScalar(m_matT.coeff(iu,iu), 0.0); m_eivalues.coeffRef(iu) = ComplexScalar(m_matT.coeff(iu,iu), 0.0);
iu--; iu--;
iter = 0; iter = 0;
@ -227,6 +228,8 @@ inline void RealSchur<MatrixType>::splitOffTwoRows(int iu, Scalar exshift)
m_eivalues.coeffRef(iu-1) = ComplexScalar(m_matT.coeff(iu,iu) + p, z); m_eivalues.coeffRef(iu-1) = ComplexScalar(m_matT.coeff(iu,iu) + p, z);
m_eivalues.coeffRef(iu) = ComplexScalar(m_matT.coeff(iu,iu) + p, -z); m_eivalues.coeffRef(iu) = ComplexScalar(m_matT.coeff(iu,iu) + p, -z);
} }
m_matT.block(iu-1, std::max(0,iu-3), 2, iu-1 - std::max(0,iu-3)).setZero();
} }
/** \internal Form shift in shiftInfo, and update exshift if an exceptional shift is performed. */ /** \internal Form shift in shiftInfo, and update exshift if an exceptional shift is performed. */
@ -273,26 +276,22 @@ inline void RealSchur<MatrixType>::computeShift(int iu, int iter, Scalar& exshif
template<typename MatrixType> template<typename MatrixType>
inline void RealSchur<MatrixType>::initFrancisQRStep(int il, int iu, const Vector3s& shiftInfo, int& im, Vector3s& firstHouseholderVector) inline void RealSchur<MatrixType>::initFrancisQRStep(int il, int iu, const Vector3s& shiftInfo, int& im, Vector3s& firstHouseholderVector)
{ {
Scalar p = 0, q = 0, r = 0; Vector3s& v = firstHouseholderVector; // alias to save typing
for (im = iu-2; im >= il; --im) for (im = iu-2; im >= il; --im)
{ {
Scalar z = m_matT.coeff(im,im); const Scalar Tmm = m_matT.coeff(im,im);
r = shiftInfo.coeff(0) - z; const Scalar r = shiftInfo.coeff(0) - Tmm;
Scalar s = shiftInfo.coeff(1) - z; const Scalar s = shiftInfo.coeff(1) - Tmm;
p = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1); v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
q = m_matT.coeff(im+1,im+1) - z - r - s; v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s;
r = m_matT.coeff(im+2,im+1); v.coeffRef(2) = m_matT.coeff(im+2,im+1);
s = ei_abs(p) + ei_abs(q) + ei_abs(r);
p = p / s;
q = q / s;
r = r / s;
if (im == il) { if (im == il) {
break; break;
} }
if (ei_abs(m_matT.coeff(im,im-1)) * (ei_abs(q) + ei_abs(r)) < const Scalar lhs = m_matT.coeff(im,im-1) * (ei_abs(v.coeff(1)) + ei_abs(v.coeff(2)));
NumTraits<Scalar>::epsilon() * (ei_abs(p) * (ei_abs(m_matT.coeff(im-1,im-1)) + ei_abs(z) + const Scalar rhs = v.coeff(0) * (ei_abs(m_matT.coeff(im-1,im-1)) + ei_abs(Tmm) + ei_abs(m_matT.coeff(im+1,im+1)));
ei_abs(m_matT.coeff(im+1,im+1))))) if (ei_abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
{ {
break; break;
} }
@ -304,11 +303,9 @@ inline void RealSchur<MatrixType>::initFrancisQRStep(int il, int iu, const Vecto
if (i > im+2) if (i > im+2)
m_matT.coeffRef(i,i-3) = 0.0; m_matT.coeffRef(i,i-3) = 0.0;
} }
firstHouseholderVector << p, q, r;
} }
/** Perform a Francis QR step involving rows il:iu and columns im:iu. */ /** \internal Perform a Francis QR step involving rows il:iu and columns im:iu. */
template<typename MatrixType> template<typename MatrixType>
inline void RealSchur<MatrixType>::performFrancisQRStep(int il, int im, int iu, const Vector3s& firstHouseholderVector, Scalar* workspace) inline void RealSchur<MatrixType>::performFrancisQRStep(int il, int im, int iu, const Vector3s& firstHouseholderVector, Scalar* workspace)
{ {