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

* rename JacobiRotation => PlanarRotation

Browse Source 
* move the makeJacobi and make_givens_* to PlanarRotation
* rename applyJacobi* => apply*
This commit is contained in:
Gael Guennebaud
2009-09-02 15:04:10 +02:00
parent 496ea63972
commit b83654b5d0
8 changed files with 194 additions and 168 deletions
Download Patch File Download Diff File Expand all files Collapse all files

59
Eigen/src/QR/ComplexSchur.h
View File

@@ -80,34 +80,6 @@ template<typename _MatrixType> class ComplexSchur
bool m_isInitialized;
};
// computes the plane rotation G such that G' x |p| = | c s' |' |p| = |z|
// |q| |-s c' | |q| |0|
// and returns z if requested. Note that the returned c is real.
template<typename T> void ei_make_givens(const std::complex<T>& p, const std::complex<T>& q,
JacobiRotation<std::complex<T> >& rot, std::complex<T>* z=0)
{
typedef std::complex<T> Complex;
T scale, absx, absxy;
if(p==Complex(0))
{
// return identity
rot.c() = Complex(1,0);
rot.s() = Complex(0,0);
if(z) *z = p;
}
else
{
scale = cnorm1(p);
absx = scale * ei_sqrt(ei_abs2(p/scale));
scale = ei_abs(scale) + cnorm1(q);
absxy = scale * ei_sqrt((absx/scale)*(absx/scale) + ei_abs2(q/scale));
rot.c() = Complex(absx / absxy);
Complex np = p/absx;
rot.s() = -ei_conj(np) * q / absxy;
if(z) *z = np * absxy;
}
}
template<typename MatrixType>
void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
{
@@ -133,8 +105,8 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
//locate the range in which to iterate
while(iu > 0)
{
d = cnorm1(m_matT.coeffRef(iu,iu)) + cnorm1(m_matT.coeffRef(iu-1,iu-1));
sd = cnorm1(m_matT.coeffRef(iu,iu-1));
d = ei_norm1(m_matT.coeffRef(iu,iu)) + ei_norm1(m_matT.coeffRef(iu-1,iu-1));
sd = ei_norm1(m_matT.coeffRef(iu,iu-1));
if(sd >= eps * d) break; // FIXME : precision criterion ??
@@ -156,8 +128,8 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
while( il > 0 )
{
// check if the current 2x2 block on the diagonal is upper triangular
d = cnorm1(m_matT.coeffRef(il,il)) + cnorm1(m_matT.coeffRef(il-1,il-1));
sd = cnorm1(m_matT.coeffRef(il,il-1));
d = ei_norm1(m_matT.coeffRef(il,il)) + ei_norm1(m_matT.coeffRef(il-1,il-1));
sd = ei_norm1(m_matT.coeffRef(il,il-1));
if(sd < eps * d) break; // FIXME : precision criterion ??
@@ -179,32 +151,32 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
r1 = (b+disc)/RealScalar(2);
r2 = (b-disc)/RealScalar(2);
if(cnorm1(r1) > cnorm1(r2))
if(ei_norm1(r1) > ei_norm1(r2))
r2 = c/r1;
else
r1 = c/r2;
if(cnorm1(r1-t.coeff(1,1)) < cnorm1(r2-t.coeff(1,1)))
if(ei_norm1(r1-t.coeff(1,1)) < ei_norm1(r2-t.coeff(1,1)))
kappa = sf * r1;
else
kappa = sf * r2;
// perform the QR step using Givens rotations
JacobiRotation<Complex> rot;
ei_make_givens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il), rot);
PlanarRotation<Complex> rot;
rot.makeGivens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il));
for(int i=il ; i<iu ; i++)
{
m_matT.block(0,i,n,n-i).applyJacobiOnTheLeft(i, i+1, rot.adjoint());
m_matT.block(0,0,std::min(i+2,iu)+1,n).applyJacobiOnTheRight(i, i+1, rot);
m_matU.applyJacobiOnTheRight(i, i+1, rot);
m_matT.block(0,i,n,n-i).applyOnTheLeft(i, i+1, rot.adjoint());
m_matT.block(0,0,std::min(i+2,iu)+1,n).applyOnTheRight(i, i+1, rot);
m_matU.applyOnTheRight(i, i+1, rot);
if(i != iu-1)
{
int i1 = i+1;
int i2 = i+2;
ei_make_givens(m_matT.coeffRef(i1,i), m_matT.coeffRef(i2,i), rot, &m_matT.coeffRef(i1,i));
rot.makeGivens(m_matT.coeffRef(i1,i), m_matT.coeffRef(i2,i), &m_matT.coeffRef(i1,i));
m_matT.coeffRef(i2,i) = Complex(0);
}
}
@@ -223,13 +195,6 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
m_isInitialized = true;
}
// norm1 of complex numbers
template<typename T>
T cnorm1(const std::complex<T> &Z)
{
return(ei_abs(Z.real()) + ei_abs(Z.imag()));
}
/**
* Computes the principal value of the square root of the complex \a z.
*/

45
Eigen/src/QR/SelfAdjointEigenSolver.h
View File

@@ -135,28 +135,6 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
#ifndef EIGEN_HIDE_HEAVY_CODE
// from Golub's "Matrix Computations", algorithm 5.1.3
template<typename Scalar>
static void ei_givens_rotation(Scalar a, Scalar b, Scalar& c, Scalar& s)
{
if (b==0)
{
c = 1; s = 0;
}
else if (ei_abs(b)>ei_abs(a))
{
Scalar t = -a/b;
s = Scalar(1)/ei_sqrt(1+t*t);
c = s * t;
}
else
{
Scalar t = -b/a;
c = Scalar(1)/ei_sqrt(1+t*t);
s = c * t;
}
}
/** \internal
*
* \qr_module
@@ -353,34 +331,33 @@ static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int st
for (int k = start; k < end; ++k)
{
RealScalar c, s;
ei_givens_rotation(x, z, c, s);
PlanarRotation<RealScalar> rot;
rot.makeGivens(x, z);
// do T = G' T G
RealScalar sdk = s * diag[k] + c * subdiag[k];
RealScalar dkp1 = s * subdiag[k] + c * diag[k+1];
RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
diag[k] = c * (c * diag[k] - s * subdiag[k]) - s * (c * subdiag[k] - s * diag[k+1]);
diag[k+1] = s * sdk + c * dkp1;
subdiag[k] = c * sdk - s * dkp1;
diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
if (k > start)
subdiag[k - 1] = c * subdiag[k-1] - s * z;
subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
x = subdiag[k];
if (k < end - 1)
{
z = -s * subdiag[k+1];
subdiag[k + 1] = c * subdiag[k+1];
z = -rot.s() * subdiag[k+1];
subdiag[k + 1] = rot.c() * subdiag[k+1];
}
// apply the givens rotation to the unit matrix Q = Q * G
// G only modifies the two columns k and k+1
if (matrixQ)
{
Map<Matrix<Scalar,Dynamic,Dynamic> > q(matrixQ,n,n);
q.applyJacobiOnTheRight(k,k+1,JacobiRotation<RealScalar>(c,s));
q.applyOnTheRight(k,k+1,rot);
}
}
}