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

RealSchur: change parameter lists; minor rewrite of computeShift().

Browse Source
This commit is contained in:
Jitse Niesen 2010-04-07 17:29:12 +01:00
parent b6829e1d5b
commit 7dea3a33a5
1 changed files with 75 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -96,11 +96,11 @@ template<typename _MatrixType> class RealSchur
typedef Matrix<Scalar,3,1> Vector3s; typedef Matrix<Scalar,3,1> Vector3s;
Scalar computeNormOfT(); Scalar computeNormOfT();
int findSmallSubdiagEntry(int n, Scalar norm); int findSmallSubdiagEntry(int iu, Scalar norm);
void computeShift(Scalar& x, Scalar& y, Scalar& w, int iu, Scalar& exshift, int iter);
void findTwoSmallSubdiagEntries(Scalar x, Scalar y, Scalar w, int il, int& m, int iu, Vector3s& firstHouseholderVector);
void doFrancisStep(int il, int m, int iu, const Vector3s& firstHouseholderVector, Scalar* workspace);
void splitOffTwoRows(int iu, Scalar exshift); void splitOffTwoRows(int iu, Scalar exshift);
void computeShift(int iu, int iter, Scalar& exshift, Vector3s& shiftInfo);
void initFrancisQRStep(int il, int iu, const Vector3s& shiftInfo, int& im, Vector3s& firstHouseholderVector);
void performFrancisQRStep(int il, int im, int iu, const Vector3s& firstHouseholderVector, Scalar* workspace);
}; };
@ -125,10 +125,10 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix)
// Rows il,...,iu is the part we are working on (the active window). // Rows il,...,iu is the part we are working on (the active window).
// Rows iu+1,...,end are already brought in triangular form. // Rows iu+1,...,end are already brought in triangular form.
int iu = m_matU.cols() - 1; int iu = m_matU.cols() - 1;
Scalar exshift = 0.0; int iter = 0; // iteration count
Scalar exshift = 0.0; // sum of exceptional shifts
Scalar norm = computeNormOfT(); Scalar norm = computeNormOfT();
int iter = 0;
while (iu >= 0) while (iu >= 0)
{ {
int il = findSmallSubdiagEntry(iu, norm); int il = findSmallSubdiagEntry(iu, norm);
@ -149,33 +149,33 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix)
} }
else // No convergence yet else // No convergence yet
{ {
Scalar x, y, w; Vector3s firstHouseholderVector, shiftInfo;
Vector3s firstHouseholderVector; computeShift(iu, iter, exshift, shiftInfo);
computeShift(x, y, w, iu, exshift, iter);
iter = iter + 1; // (Could check iteration count here.) iter = iter + 1; // (Could check iteration count here.)
int m; int im;
findTwoSmallSubdiagEntries(x, y, w, il, m, iu, firstHouseholderVector); initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
doFrancisStep(il, m, iu, firstHouseholderVector, workspace); performFrancisQRStep(il, im, iu, firstHouseholderVector, workspace);
} // check convergence }
} // while (iu >= 0) }
m_isInitialized = true; m_isInitialized = true;
} }
// Compute matrix norm /** \internal Computes and returns vector L1 norm of T */
template<typename MatrixType> template<typename MatrixType>
inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT() inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
{ {
const int size = m_matU.cols(); const int size = m_matU.cols();
// FIXME to be efficient the following would requires a triangular reduxion code // FIXME to be efficient the following would requires a triangular reduxion code
// Scalar norm = m_matT.upper().cwiseAbs().sum() + m_matT.corner(BottomLeft,size-1,size-1).diagonal().cwiseAbs().sum(); // Scalar norm = m_matT.upper().cwiseAbs().sum()
// + m_matT.corner(BottomLeft,size-1,size-1).diagonal().cwiseAbs().sum();
Scalar norm = 0.0; Scalar norm = 0.0;
for (int j = 0; j < size; ++j) for (int j = 0; j < size; ++j)
norm += m_matT.row(j).segment(std::max(j-1,0), size-std::max(j-1,0)).cwiseAbs().sum(); norm += m_matT.row(j).segment(std::max(j-1,0), size-std::max(j-1,0)).cwiseAbs().sum();
return norm; return norm;
} }
// Look for single small sub-diagonal element /** \internal Look for single small sub-diagonal element and returns its index */
template<typename MatrixType> template<typename MatrixType>
inline int RealSchur<MatrixType>::findSmallSubdiagEntry(int iu, Scalar norm) inline int RealSchur<MatrixType>::findSmallSubdiagEntry(int iu, Scalar norm)
{ {
@ -192,133 +192,134 @@ inline int RealSchur<MatrixType>::findSmallSubdiagEntry(int iu, Scalar norm)
return res; return res;
} }
/** \internal Update T given that rows iu-1 and iu decouple from the rest. */
template<typename MatrixType> template<typename MatrixType>
inline void RealSchur<MatrixType>::splitOffTwoRows(int iu, Scalar exshift) inline void RealSchur<MatrixType>::splitOffTwoRows(int iu, Scalar exshift)
{ {
const int size = m_matU.cols(); const int size = m_matU.cols();
// The eigenvalues of the 2x2 matrix [a b; c d] are
// trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc
Scalar w = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); Scalar w = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
Scalar p = (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu)) * Scalar(0.5); Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
Scalar q = p * p + w; Scalar q = p * p + w; // q = tr^2 / 4 - det = discr/4
Scalar z = ei_sqrt(ei_abs(q)); Scalar z = ei_sqrt(ei_abs(q));
m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift; m_matT.coeffRef(iu,iu) += exshift;
m_matT.coeffRef(iu-1,iu-1) = m_matT.coeff(iu-1,iu-1) + exshift; m_matT.coeffRef(iu-1,iu-1) += exshift;
Scalar x = m_matT.coeff(iu,iu);
// Scalar pair if (q >= 0) // Two real eigenvalues
if (q >= 0)
{ {
if (p >= 0)
z = p + z;
else
z = p - z;
m_eivalues.coeffRef(iu-1) = ComplexScalar(x + z, 0.0);
m_eivalues.coeffRef(iu) = ComplexScalar(z!=0.0 ? x - w / z : m_eivalues.coeff(iu-1).real(), 0.0);
PlanarRotation<Scalar> rot; PlanarRotation<Scalar> rot;
rot.makeGivens(z, m_matT.coeff(iu, iu-1)); if (p >= 0)
rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
else
rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
m_matT.block(0, iu-1, size, size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint()); m_matT.block(0, iu-1, size, size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
m_matT.block(0, 0, iu+1, size).applyOnTheRight(iu-1, iu, rot); m_matT.block(0, 0, iu+1, size).applyOnTheRight(iu-1, iu, rot);
m_matU.applyOnTheRight(iu-1, iu, rot); m_matU.applyOnTheRight(iu-1, iu, rot);
m_eivalues.coeffRef(iu-1) = ComplexScalar(m_matT.coeff(iu-1, iu-1), 0.0);
m_eivalues.coeffRef(iu) = ComplexScalar(m_matT.coeff(iu, iu), 0.0);
} }
else // Complex pair else // // Pair of complex conjugate eigenvalues
{ {
m_eivalues.coeffRef(iu-1) = ComplexScalar(x + p, z); m_eivalues.coeffRef(iu-1) = ComplexScalar(m_matT.coeff(iu,iu) + p, z);
m_eivalues.coeffRef(iu) = ComplexScalar(x + p, -z); m_eivalues.coeffRef(iu) = ComplexScalar(m_matT.coeff(iu,iu) + p, -z);
} }
} }
// Form shift /** \internal Form shift in shiftInfo, and update exshift if an exceptional shift is performed. */
template<typename MatrixType> template<typename MatrixType>
inline void RealSchur<MatrixType>::computeShift(Scalar& x, Scalar& y, Scalar& w, int iu, Scalar& exshift, int iter) inline void RealSchur<MatrixType>::computeShift(int iu, int iter, Scalar& exshift, Vector3s& shiftInfo)
{ {
x = m_matT.coeff(iu,iu); shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
y = m_matT.coeff(iu-1,iu-1); shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
w = 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 // Wilkinson's original ad hoc shift
if (iter == 10) if (iter == 10)
{ {
exshift += x; exshift += shiftInfo.coeff(0);
for (int i = 0; i <= iu; ++i) for (int i = 0; i <= iu; ++i)
m_matT.coeffRef(i,i) -= x; m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
Scalar s = ei_abs(m_matT.coeff(iu,iu-1)) + ei_abs(m_matT.coeff(iu-1,iu-2)); Scalar s = ei_abs(m_matT.coeff(iu,iu-1)) + ei_abs(m_matT.coeff(iu-1,iu-2));
x = y = Scalar(0.75) * s; shiftInfo.coeffRef(0) = Scalar(0.75) * s;
w = Scalar(-0.4375) * s * s; shiftInfo.coeffRef(1) = Scalar(0.75) * s;
shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
} }
// MATLAB's new ad hoc shift // MATLAB's new ad hoc shift
if (iter == 30) if (iter == 30)
{ {
Scalar s = Scalar((y - x) / 2.0); Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
s = s * s + w; s = s * s + shiftInfo.coeff(2);
if (s > 0) if (s > 0)
{ {
s = ei_sqrt(s); s = ei_sqrt(s);
if (y < x) if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
s = -s; s = -s;
s = Scalar(x - w / ((y - x) / 2.0 + s)); s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
exshift += s;
for (int i = 0; i <= iu; ++i) for (int i = 0; i <= iu; ++i)
m_matT.coeffRef(i,i) -= s; m_matT.coeffRef(i,i) -= s;
exshift += s; shiftInfo.setConstant(Scalar(0.964));
x = y = w = Scalar(0.964);
} }
} }
} }
// Look for two consecutive small sub-diagonal elements /** \internal Compute index im at which Francis QR step starts and the first Householder vector. */
template<typename MatrixType> template<typename MatrixType>
inline void RealSchur<MatrixType>::findTwoSmallSubdiagEntries(Scalar x, Scalar y, Scalar w, int il, int& m, int iu, 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; Scalar p = 0, q = 0, r = 0;
m = iu-2; for (im = iu-2; im >= il; --im)
while (m >= il)
{ {
Scalar z = m_matT.coeff(m,m); Scalar z = m_matT.coeff(im,im);
r = x - z; r = shiftInfo.coeff(0) - z;
Scalar s = y - z; Scalar s = shiftInfo.coeff(1) - z;
p = (r * s - w) / m_matT.coeff(m+1,m) + m_matT.coeff(m,m+1); p = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
q = m_matT.coeff(m+1,m+1) - z - r - s; q = m_matT.coeff(im+1,im+1) - z - r - s;
r = m_matT.coeff(m+2,m+1); r = m_matT.coeff(im+2,im+1);
s = ei_abs(p) + ei_abs(q) + ei_abs(r); s = ei_abs(p) + ei_abs(q) + ei_abs(r);
p = p / s; p = p / s;
q = q / s; q = q / s;
r = r / s; r = r / s;
if (m == il) { if (im == il) {
break; break;
} }
if (ei_abs(m_matT.coeff(m,m-1)) * (ei_abs(q) + ei_abs(r)) < if (ei_abs(m_matT.coeff(im,im-1)) * (ei_abs(q) + ei_abs(r)) <
NumTraits<Scalar>::epsilon() * (ei_abs(p) * (ei_abs(m_matT.coeff(m-1,m-1)) + ei_abs(z) + NumTraits<Scalar>::epsilon() * (ei_abs(p) * (ei_abs(m_matT.coeff(im-1,im-1)) + ei_abs(z) +
ei_abs(m_matT.coeff(m+1,m+1))))) ei_abs(m_matT.coeff(im+1,im+1)))))
{ {
break; break;
} }
m--;
} }
for (int i = m+2; i <= iu; ++i) for (int i = im+2; i <= iu; ++i)
{ {
m_matT.coeffRef(i,i-2) = 0.0; m_matT.coeffRef(i,i-2) = 0.0;
if (i > m+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; firstHouseholderVector << p, q, r;
} }
// Double QR step involving rows il:iu and columns m:iu /** Perform a Francis QR step involving rows il:iu and columns im:iu. */
template<typename MatrixType> template<typename MatrixType>
inline void RealSchur<MatrixType>::doFrancisStep(int il, int m, int iu, const Vector3s& firstHouseholderVector, Scalar* workspace) inline void RealSchur<MatrixType>::performFrancisQRStep(int il, int im, int iu, const Vector3s& firstHouseholderVector, Scalar* workspace)
{ {
assert(m >= il); assert(im >= il);
assert(m <= iu-2); assert(im <= iu-2);
const int size = m_matU.cols(); const int size = m_matU.cols();
for (int k = m; k <= iu-2; ++k) for (int k = im; k <= iu-2; ++k)
{ {
bool firstIteration = (k == m); bool firstIteration = (k == im);
Vector3s v; Vector3s v;
if (firstIteration) if (firstIteration)