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

Take advantage that T is already diagonal in the extraction of generalized complex eigenvalues.

Browse Source
This commit is contained in:
Gael Guennebaud 2016-06-09 17:12:03 +02:00
parent c1f9ca9254
commit 2bd59b0e0d
1 changed files with 7 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
Eigen/src/Eigenvalues/GeneralizedEigenSolver.h
View File

@ -327,24 +327,13 @@ GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixTyp
}
else
{
// We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a triangular 2x2 block T
// From the eigen decomposition of T = U * E * U^-1,
// we can extract the eigenvalues of (U^-1 * S * U) / E
// Here, we can take advantage that E = diag(T), and U = [ 1 T_01 ; 0 T_11-T_00], and U^-1 = [1 -T_11/(T_11-T_00) ; 0 1/(T_11-T_00)].
// Then taking beta=T_00*T_11*(T_11-T_00), we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00) * (T_11-T_00):
// We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T
// Then taking beta=T_00*T_11, we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00):
// T = [a b ; 0 c]
// S = [e f ; g h]
RealScalar a = m_realQZ.matrixT().coeff(i, i), b = m_realQZ.matrixT().coeff(i, i+1), c = m_realQZ.matrixT().coeff(i+1, i+1);
RealScalar e = m_matS.coeff(i, i), f = m_matS.coeff(i, i+1), g = m_matS.coeff(i+1, i), h = m_matS.coeff(i+1, i+1);
RealScalar d = c-a;
RealScalar gb = g*b;
Matrix<RealScalar,2,2> S2;
S2 << (e*d-gb)*c, ((e*b+f*d-h*b)*d-gb*b)*a,
g*c , (gb+h*d)*a;
// NOTE, we could also compute the SVD of T's block during the QZ factorization so that the respective T block is guaranteed to be diagonal,
// and then we could directly apply the formula below (while taking care of scaling S columns by T11,T00):
// T = [a 0]
// [0 b]
RealScalar a = m_realQZ.matrixT().coeff(i, i), b = m_realQZ.matrixT().coeff(i+1, i+1);
Matrix<RealScalar,2,2> S2 = m_matS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal();
Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1));
Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1)));
@ -352,7 +341,7 @@ GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixTyp
m_alphas.coeffRef(i+1) = ComplexScalar(S2.coeff(1,1) + p, -z);
m_betas.coeffRef(i) =
m_betas.coeffRef(i+1) = a*c*d;
m_betas.coeffRef(i+1) = a*b;
i += 2;
}