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: Rename l and n to il and iu.

Browse Source
This commit is contained in:
Jitse Niesen 2010-04-06 18:26:30 +01:00
parent 9fad1e392b
commit cc57df9bea
1 changed files with 62 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -95,10 +95,10 @@ template<typename _MatrixType> class RealSchur
Scalar computeNormOfT(); Scalar computeNormOfT();
int findSmallSubdiagEntry(int n, Scalar norm); int findSmallSubdiagEntry(int n, Scalar norm);
void computeShift(Scalar& x, Scalar& y, Scalar& w, int l, int n, Scalar& exshift, int iter); void computeShift(Scalar& x, Scalar& y, Scalar& w, int iu, Scalar& exshift, int iter);
void findTwoSmallSubdiagEntries(Scalar x, Scalar y, Scalar w, int l, int& m, int n, Scalar& p, Scalar& q, Scalar& r); void findTwoSmallSubdiagEntries(Scalar x, Scalar y, Scalar w, int l, int& m, int iu, Scalar& p, Scalar& q, Scalar& r);
void doFrancisStep(int l, int m, int n, Scalar p, Scalar q, Scalar r, Scalar x, Scalar* workspace); void doFrancisStep(int l, int m, int iu, Scalar p, Scalar q, Scalar r, Scalar x, Scalar* workspace);
void splitOffTwoRows(int n, Scalar exshift); void splitOffTwoRows(int iu, Scalar exshift);
}; };
@ -118,39 +118,43 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix)
ColumnVectorType workspaceVector(m_matU.cols()); ColumnVectorType workspaceVector(m_matU.cols());
Scalar* workspace = &workspaceVector.coeffRef(0); Scalar* workspace = &workspaceVector.coeffRef(0);
int n = m_matU.cols() - 1; // The matrix m_matT is divided in three parts.
// Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
// Rows il,...,iu is the part we are working on (the active window).
// Rows iu+1,...,end are already brought in triangular form.
int iu = m_matU.cols() - 1;
Scalar exshift = 0.0; Scalar exshift = 0.0;
Scalar norm = computeNormOfT(); Scalar norm = computeNormOfT();
int iter = 0; int iter = 0;
while (n >= 0) while (iu >= 0)
{ {
int l = findSmallSubdiagEntry(n, norm); int il = findSmallSubdiagEntry(iu, norm);
// Check for convergence // Check for convergence
if (l == n) // One root found if (il == iu) // One root found
{ {
m_matT.coeffRef(n,n) = m_matT.coeff(n,n) + exshift; m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
m_eivalues.coeffRef(n) = ComplexScalar(m_matT.coeff(n,n), 0.0); m_eivalues.coeffRef(iu) = ComplexScalar(m_matT.coeff(iu,iu), 0.0);
n--; iu--;
iter = 0; iter = 0;
} }
else if (l == n-1) // Two roots found else if (il == iu-1) // Two roots found
{ {
splitOffTwoRows(n, exshift); splitOffTwoRows(iu, exshift);
n = n - 2; iu -= 2;
iter = 0; iter = 0;
} }
else // No convergence yet else // No convergence yet
{ {
Scalar p = 0, q = 0, r = 0, x, y, w; Scalar p = 0, q = 0, r = 0, x, y, w;
computeShift(x, y, w, l, n, exshift, iter); 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 m;
findTwoSmallSubdiagEntries(x, y, w, l, m, n, p, q, r); findTwoSmallSubdiagEntries(x, y, w, il, m, iu, p, q, r);
doFrancisStep(l, m, n, p, q, r, x, workspace); doFrancisStep(il, m, iu, p, q, r, x, workspace);
} // check convergence } // check convergence
} // while (n >= 0) } // while (iu >= 0)
m_isInitialized = true; m_isInitialized = true;
} }
@ -170,32 +174,32 @@ inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
// Look for single small sub-diagonal element // Look for single small sub-diagonal element
template<typename MatrixType> template<typename MatrixType>
inline int RealSchur<MatrixType>::findSmallSubdiagEntry(int n, Scalar norm) inline int RealSchur<MatrixType>::findSmallSubdiagEntry(int iu, Scalar norm)
{ {
int l = n; int res = iu;
while (l > 0) while (res > 0)
{ {
Scalar s = ei_abs(m_matT.coeff(l-1,l-1)) + ei_abs(m_matT.coeff(l,l)); Scalar s = ei_abs(m_matT.coeff(res-1,res-1)) + ei_abs(m_matT.coeff(res,res));
if (s == 0.0) if (s == 0.0)
s = norm; s = norm;
if (ei_abs(m_matT.coeff(l,l-1)) < NumTraits<Scalar>::epsilon() * s) if (ei_abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
break; break;
l--; res--;
} }
return l; return res;
} }
template<typename MatrixType> template<typename MatrixType>
inline void RealSchur<MatrixType>::splitOffTwoRows(int n, Scalar exshift) inline void RealSchur<MatrixType>::splitOffTwoRows(int iu, Scalar exshift)
{ {
const int size = m_matU.cols(); const int size = m_matU.cols();
Scalar w = m_matT.coeff(n,n-1) * m_matT.coeff(n-1,n); Scalar w = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
Scalar p = (m_matT.coeff(n-1,n-1) - m_matT.coeff(n,n)) * Scalar(0.5); Scalar p = (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu)) * Scalar(0.5);
Scalar q = p * p + w; Scalar q = p * p + w;
Scalar z = ei_sqrt(ei_abs(q)); Scalar z = ei_sqrt(ei_abs(q));
m_matT.coeffRef(n,n) = m_matT.coeff(n,n) + exshift; m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
m_matT.coeffRef(n-1,n-1) = m_matT.coeff(n-1,n-1) + exshift; m_matT.coeffRef(iu-1,iu-1) = m_matT.coeff(iu-1,iu-1) + exshift;
Scalar x = m_matT.coeff(n,n); Scalar x = m_matT.coeff(iu,iu);
// Scalar pair // Scalar pair
if (q >= 0) if (q >= 0)
@ -205,42 +209,37 @@ inline void RealSchur<MatrixType>::splitOffTwoRows(int n, Scalar exshift)
else else
z = p - z; z = p - z;
m_eivalues.coeffRef(n-1) = ComplexScalar(x + z, 0.0); m_eivalues.coeffRef(iu-1) = ComplexScalar(x + z, 0.0);
m_eivalues.coeffRef(n) = ComplexScalar(z!=0.0 ? x - w / z : m_eivalues.coeff(n-1).real(), 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(n, n-1)); rot.makeGivens(z, m_matT.coeff(iu, iu-1));
m_matT.block(0, n-1, size, size-n+1).applyOnTheLeft(n-1, n, rot.adjoint()); m_matT.block(0, iu-1, size, size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
m_matT.block(0, 0, n+1, size).applyOnTheRight(n-1, n, rot); m_matT.block(0, 0, iu+1, size).applyOnTheRight(iu-1, iu, rot);
m_matU.applyOnTheRight(n-1, n, rot); m_matU.applyOnTheRight(iu-1, iu, rot);
} }
else // Complex pair else // Complex pair
{ {
m_eivalues.coeffRef(n-1) = ComplexScalar(x + p, z); m_eivalues.coeffRef(iu-1) = ComplexScalar(x + p, z);
m_eivalues.coeffRef(n) = ComplexScalar(x + p, -z); m_eivalues.coeffRef(iu) = ComplexScalar(x + p, -z);
} }
} }
// Form shift // Form shift
template<typename MatrixType> template<typename MatrixType>
inline void RealSchur<MatrixType>::computeShift(Scalar& x, Scalar& y, Scalar& w, int l, int n, Scalar& exshift, int iter) inline void RealSchur<MatrixType>::computeShift(Scalar& x, Scalar& y, Scalar& w, int iu, Scalar& exshift, int iter)
{ {
x = m_matT.coeff(n,n); x = m_matT.coeff(iu,iu);
y = 0.0; y = m_matT.coeff(iu-1,iu-1);
w = 0.0; w = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
if (l < n)
{
y = m_matT.coeff(n-1,n-1);
w = m_matT.coeff(n,n-1) * m_matT.coeff(n-1,n);
}
// Wilkinson's original ad hoc shift // Wilkinson's original ad hoc shift
if (iter == 10) if (iter == 10)
{ {
exshift += x; exshift += x;
for (int i = 0; i <= n; ++i) for (int i = 0; i <= iu; ++i)
m_matT.coeffRef(i,i) -= x; m_matT.coeffRef(i,i) -= x;
Scalar s = ei_abs(m_matT.coeff(n,n-1)) + ei_abs(m_matT.coeff(n-1,n-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; x = y = Scalar(0.75) * s;
w = Scalar(-0.4375) * s * s; w = Scalar(-0.4375) * s * s;
} }
@ -256,7 +255,7 @@ inline void RealSchur<MatrixType>::computeShift(Scalar& x, Scalar& y, Scalar& w,
if (y < x) if (y < x)
s = -s; s = -s;
s = Scalar(x - w / ((y - x) / 2.0 + s)); s = Scalar(x - w / ((y - x) / 2.0 + s));
for (int i = 0; i <= n; ++i) for (int i = 0; i <= iu; ++i)
m_matT.coeffRef(i,i) -= s; m_matT.coeffRef(i,i) -= s;
exshift += s; exshift += s;
x = y = w = Scalar(0.964); x = y = w = Scalar(0.964);
@ -266,10 +265,10 @@ inline void RealSchur<MatrixType>::computeShift(Scalar& x, Scalar& y, Scalar& w,
// Look for two consecutive small sub-diagonal elements // Look for two consecutive small sub-diagonal elements
template<typename MatrixType> template<typename MatrixType>
inline void RealSchur<MatrixType>::findTwoSmallSubdiagEntries(Scalar x, Scalar y, Scalar w, int l, int& m, int n, Scalar& p, Scalar& q, Scalar& r) inline void RealSchur<MatrixType>::findTwoSmallSubdiagEntries(Scalar x, Scalar y, Scalar w, int il, int& m, int iu, Scalar& p, Scalar& q, Scalar& r)
{ {
m = n-2; m = iu-2;
while (m >= l) while (m >= il)
{ {
Scalar z = m_matT.coeff(m,m); Scalar z = m_matT.coeff(m,m);
r = x - z; r = x - z;
@ -281,7 +280,7 @@ inline void RealSchur<MatrixType>::findTwoSmallSubdiagEntries(Scalar x, Scalar y
p = p / s; p = p / s;
q = q / s; q = q / s;
r = r / s; r = r / s;
if (m == l) { if (m == il) {
break; break;
} }
if (ei_abs(m_matT.coeff(m,m-1)) * (ei_abs(q) + ei_abs(r)) < if (ei_abs(m_matT.coeff(m,m-1)) * (ei_abs(q) + ei_abs(r)) <
@ -293,7 +292,7 @@ inline void RealSchur<MatrixType>::findTwoSmallSubdiagEntries(Scalar x, Scalar y
m--; m--;
} }
for (int i = m+2; i <= n; ++i) for (int i = m+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 > m+2)
@ -301,15 +300,15 @@ inline void RealSchur<MatrixType>::findTwoSmallSubdiagEntries(Scalar x, Scalar y
} }
} }
// Double QR step involving rows l:n and columns m:n // Double QR step involving rows il:iu and columns m:iu
template<typename MatrixType> template<typename MatrixType>
inline void RealSchur<MatrixType>::doFrancisStep(int l, int m, int n, Scalar p, Scalar q, Scalar r, Scalar x, Scalar* workspace) inline void RealSchur<MatrixType>::doFrancisStep(int il, int m, int iu, Scalar p, Scalar q, Scalar r, Scalar x, Scalar* workspace)
{ {
const int size = m_matU.cols(); const int size = m_matU.cols();
for (int k = m; k <= n-1; ++k) for (int k = m; k <= iu-1; ++k)
{ {
int notlast = (k != n-1); int notlast = (k != iu-1);
if (k != m) { if (k != m) {
p = m_matT.coeff(k,k-1); p = m_matT.coeff(k,k-1);
q = m_matT.coeff(k+1,k-1); q = m_matT.coeff(k+1,k-1);
@ -335,7 +334,7 @@ inline void RealSchur<MatrixType>::doFrancisStep(int l, int m, int n, Scalar p,
{ {
if (k != m) if (k != m)
m_matT.coeffRef(k,k-1) = -s * x; m_matT.coeffRef(k,k-1) = -s * x;
else if (l != m) else if (il != m)
m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1); m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
p = p + s; p = p + s;
@ -344,7 +343,7 @@ inline void RealSchur<MatrixType>::doFrancisStep(int l, int m, int n, Scalar p,
{ {
Matrix<Scalar, 2, 1> ess(q/p, r/p); Matrix<Scalar, 2, 1> ess(q/p, r/p);
m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, p/s, workspace); m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, p/s, workspace);
m_matT.block(0, k, std::min(n,k+3) + 1, 3).applyHouseholderOnTheRight(ess, p/s, workspace); m_matT.block(0, k, std::min(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, p/s, workspace);
m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, p/s, workspace); m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, p/s, workspace);
} }
else else
@ -352,7 +351,7 @@ inline void RealSchur<MatrixType>::doFrancisStep(int l, int m, int n, Scalar p,
Matrix<Scalar, 1, 1> ess; Matrix<Scalar, 1, 1> ess;
ess.coeffRef(0) = q/p; ess.coeffRef(0) = q/p;
m_matT.block(k, k, 2, size-k).applyHouseholderOnTheLeft(ess, p/s, workspace); m_matT.block(k, k, 2, size-k).applyHouseholderOnTheLeft(ess, p/s, workspace);
m_matT.block(0, k, std::min(n,k+3) + 1, 2).applyHouseholderOnTheRight(ess, p/s, workspace); m_matT.block(0, k, std::min(iu,k+3) + 1, 2).applyHouseholderOnTheRight(ess, p/s, workspace);
m_matU.block(0, k, size, 2).applyHouseholderOnTheRight(ess, p/s, workspace); m_matU.block(0, k, size, 2).applyHouseholderOnTheRight(ess, p/s, workspace);
} }
} // (s != 0) } // (s != 0)