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

ComplexSchur: compute shift more stably, introduce exceptional shifts.

Browse Source 
Both the new computation of  the eigenvalues of a 2x2 block and the
exceptional shifts are taken from EISPACK routine COMQR.
This commit is contained in:
Jitse Niesen 2010-02-25 22:33:38 +00:00
parent 53bae6b3f8
commit 90e4a605ef
1 changed files with 23 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
Eigen/src/Eigenvalues/ComplexSchur.h
View File

@ -154,6 +154,14 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
m_matT = hess.matrixH(); m_matT = hess.matrixH();
if(!skipU) m_matU = hess.matrixQ(); if(!skipU) m_matU = hess.matrixQ();
// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
// 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 submatrix).
// Rows iu+1,...,end are already brought in triangular form.
int iu = m_matT.cols() - 1; int iu = m_matT.cols() - 1;
int il; int il;
RealScalar d,sd,sf; RealScalar d,sd,sf;
@ -164,7 +172,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
int iter = 0; int iter = 0;
while(true) while(true)
{ {
//locate the range in which to iterate // find iu, the bottom row of the active submatrix
while(iu > 0) while(iu > 0)
{ {
d = ei_norm1(m_matT.coeff(iu,iu)) + ei_norm1(m_matT.coeff(iu-1,iu-1)); d = ei_norm1(m_matT.coeff(iu,iu)) + ei_norm1(m_matT.coeff(iu-1,iu-1));
@ -187,6 +195,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
return; return;
} }
// find il, the top row of the active submatrix
il = iu-1; il = iu-1;
while(il > 0) while(il > 0)
{ {
@ -202,15 +211,16 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
if( il != 0 ) m_matT.coeffRef(il,il-1) = Complex(0); if( il != 0 ) m_matT.coeffRef(il,il-1) = Complex(0);
// compute the shift (the normalization by sf is to avoid under/overflow) // compute the shift kappa as one of the eigenvalues of the 2x2
// diagonal block on the bottom of the active submatrix
Matrix<Scalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1); Matrix<Scalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
sf = t.cwiseAbs().sum(); sf = t.cwiseAbs().sum();
t /= sf; t /= sf; // the normalization by sf is to avoid under/overflow
c = t.determinant(); b = t.coeff(0,0) + t.coeff(1,1);
b = t.diagonal().sum(); c = t.coeff(0,0) - t.coeff(1,1);
disc = ei_sqrt(c*c + RealScalar(4)*t.coeff(0,1)*t.coeff(1,0));
disc = ei_sqrt(b*b - RealScalar(4)*c);
r1 = (b+disc)/RealScalar(2); r1 = (b+disc)/RealScalar(2);
r2 = (b-disc)/RealScalar(2); r2 = (b-disc)/RealScalar(2);
@ -225,6 +235,12 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
else else
kappa = sf * r2; kappa = sf * r2;
if (iter == 10 || iter == 20)
{
// exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
kappa = ei_abs(ei_real(m_matT.coeff(iu,iu-1))) + ei_abs(ei_real(m_matT.coeff(iu-1,iu-2)));
}
// perform the QR step using Givens rotations // perform the QR step using Givens rotations
PlanarRotation<Complex> rot; PlanarRotation<Complex> rot;
rot.makeGivens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il)); rot.makeGivens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il));
@ -246,18 +262,6 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
} }
} }
// FIXME : is it necessary ?
/*
for(int i=0 ; i<n ; i++)
for(int j=0 ; j<n ; j++)
{
if(ei_abs(ei_real(m_matT.coeff(i,j))) < eps)
ei_real_ref(m_matT.coeffRef(i,j)) = 0;
if(ei_imag(ei_abs(m_matT.coeff(i,j))) < eps)
ei_imag_ref(m_matT.coeffRef(i,j)) = 0;
}
*/
m_isInitialized = true; m_isInitialized = true;
m_matUisUptodate = !skipU; m_matUisUptodate = !skipU;
} }