From 41e5625f96b9d9642c1724ff3859e709d9cfe8cb Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 10 Jun 2010 09:34:49 +0200 Subject: [PATCH] clean general symm eigensolver --- .../src/Eigenvalues/SelfAdjointEigenSolver.h | 23 +++++-------------- 1 file changed, 6 insertions(+), 17 deletions(-) diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 0fef247e0..636aa090c 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -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 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().solveTriangularInPlace(matC); -// matC = matC.transpose(); - // FIXME: this should work: (currently it only does for small matrices) -// Transpose trMatC(matC); -// cholB.matrixL().conjugate().eval().template marked().solveTriangularInPlace(trMatC); + cholB.matrixL().template solveInPlace(matC); + cholB.matrixU().template solveInPlace(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; }