mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-03 06:50:57 +08:00
Code Issues Packages Projects Releases Wiki Activity

implement other variants

Browse Source
This commit is contained in:
Gael Guennebaud 2010-06-17 10:11:38 +02:00
parent db160f2e0b
commit 43086d12d2
2 changed files with 49 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

64
Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h
View File

@ -43,7 +43,7 @@
* This class solves the generalized eigenvalue problem
* \f$ Av = \lambda Bv \f$. In this case, the matrix \f$ A \f$ should be
* selfadjoint and the matrix \f$ B \f$ should be positive definite.
*
*
* Only the \b lower \b triangular \b part of the input matrix is referenced.
*
* Call the function compute() to compute the eigenvalues and eigenvectors of
@ -138,7 +138,7 @@ class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixT
* Only the lower triangular part of the matrix is referenced.
* \param[in] options A or-ed set of flags {ComputeEigenvectors,EigenvaluesOnly} | {Ax_lBx,ABx_lx,BAx_lx}.
* Default is ComputeEigenvectors|Ax_lBx.
*
*
* \returns Reference to \c *this
*
* If \p options contains Ax_lBx (the default), this function computes eigenvalues
@ -147,7 +147,7 @@ class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixT
* matrix \f$ A \f$ and \a matB the positive definite
* matrix \f$ B \f$. In addition, each eigenvector \f$ x \f$
* satisfies the property \f$ x^* B x = 1 \f$.
*
*
* In addition, the two following variants can be solved via \p options:
* - \c ABx_lx: \f$ ABx = \lambda x \f$
* - \c BAx_lx: \f$ BAx = \lambda x \f$
@ -185,28 +185,58 @@ compute(const MatrixType& matA, const MatrixType& matB, int options)
ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
ei_assert((options&~(EigVecMask|GenEigMask))==0
&& (options&EigVecMask)!=EigVecMask
&& ((options&GenEigMask)==Ax_lBx || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
&& ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
|| (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
&& "invalid option parameter");
ei_assert((options&GenEigMask)==Ax_lBx && "other variants are not implemented yet, sorry.");
// TODO implements other variants !!
bool computeEigVecs = (options&EigVecMask)==ComputeEigenvectors;
bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
// Compute the cholesky decomposition of matB = L L' = U'U
LLT<MatrixType> cholB(matB);
// compute C = inv(L) A inv(L')
MatrixType matC = matA.template selfadjointView<Lower>();
cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
cholB.matrixU().template solveInPlace<OnTheRight>(matC);
int type = (options&GenEigMask);
if(type==0)
type = Ax_lBx;
Base::compute(matC, options&EigVecMask);
if(type==Ax_lBx)
{
// compute C = inv(L) A inv(L')
MatrixType matC = matA.template selfadjointView<Lower>();
cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
cholB.matrixU().template solveInPlace<OnTheRight>(matC);
// transform back the eigen vectors: evecs = inv(U) * evecs
if(computeEigVecs)
cholB.matrixU().solveInPlace(Base::m_eivec);
Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
// transform back the eigen vectors: evecs = inv(U) * evecs
if(computeEigVecs)
cholB.matrixU().solveInPlace(Base::m_eivec);
}
else if(type==ABx_lx)
{
// compute C = L' A L
MatrixType matC = matA.template selfadjointView<Lower>();
matC = matC * cholB.matrixL();
matC = cholB.matrixU() * matC;
Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
// transform back the eigen vectors: evecs = inv(U) * evecs
if(computeEigVecs)
cholB.matrixU().solveInPlace(Base::m_eivec);
}
else if(type==BAx_lx)
{
// compute C = L' A L
MatrixType matC = matA.template selfadjointView<Lower>();
matC = matC * cholB.matrixL();
matC = cholB.matrixU() * matC;
Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
// transform back the eigen vectors: evecs = L * evecs
if(computeEigVecs)
Base::m_eivec = cholB.matrixL() * Base::m_eivec;
}
return *this;
}

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

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2008-2010 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// Eigen is free software; you can redistribute it and/or
@ -66,7 +66,7 @@
*
* The documentation for SelfAdjointEigenSolver(const MatrixType&, bool)
* contains an example of the typical use of this class.
*
*
* To solve the \em generalized eigenvalue problem \f$ Av = \lambda Bv \f$ and
* the like see the class GeneralizedSelfAdjointEigenSolver.
*