GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-24 02:29:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

clean general symm eigensolver

Browse Source
This commit is contained in:
Gael Guennebaud 2010-06-10 09:34:49 +02:00
parent 3f388282ae
commit 41e5625f96
1 changed files with 6 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@ -499,31 +499,20 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors
{
ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
// Compute the cholesky decomposition of matB = L L'
// Compute the cholesky decomposition of matB = L L' = U'U
LLT<MatrixType> cholB(matB);
// compute C = inv(L) A inv(L')
MatrixType matC = matA;
cholB.matrixL().solveInPlace(matC);
// FIXME since we currently do not support A * inv(L'), let's do (inv(L) A')' :
matC.adjointInPlace();
cholB.matrixL().solveInPlace(matC);
matC.adjointInPlace();
// this version works too:
// matC = matC.transpose();
// cholB.matrixL().conjugate().template marked<Lower>().solveTriangularInPlace(matC);
// matC = matC.transpose();
// FIXME: this should work: (currently it only does for small matrices)
// Transpose<MatrixType> trMatC(matC);
// cholB.matrixL().conjugate().eval().template marked<Lower>().solveTriangularInPlace(trMatC);
cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
cholB.matrixU().template solveInPlace<OnTheRight>(matC);
compute(matC, computeEigenvectors);
if (computeEigenvectors)
{
// transform back the eigen vectors: evecs = inv(U) * evecs
// transform back the eigen vectors: evecs = inv(U) * evecs
if(computeEigenvectors)
cholB.matrixU().solveInPlace(m_eivec);
}
return *this;
}