mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-07-31 09:12:02 +08:00
215 lines
9.6 KiB
C++
215 lines
9.6 KiB
C++
// This file is part of Eigen, a lightweight C++ template library
|
|
// for linear algebra.
|
|
//
|
|
// 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
|
|
// modify it under the terms of the GNU Lesser General Public
|
|
// License as published by the Free Software Foundation; either
|
|
// version 3 of the License, or (at your option) any later version.
|
|
//
|
|
// Alternatively, you can redistribute it and/or
|
|
// modify it under the terms of the GNU General Public License as
|
|
// published by the Free Software Foundation; either version 2 of
|
|
// the License, or (at your option) any later version.
|
|
//
|
|
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
|
|
// GNU General Public License for more details.
|
|
//
|
|
// You should have received a copy of the GNU Lesser General Public
|
|
// License and a copy of the GNU General Public License along with
|
|
// Eigen. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
#ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
|
|
#define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
|
|
|
|
#include "./EigenvaluesCommon.h"
|
|
#include "./Tridiagonalization.h"
|
|
|
|
/** \eigenvalues_module \ingroup Eigenvalues_Module
|
|
* \nonstableyet
|
|
*
|
|
* \class GeneralizedSelfAdjointEigenSolver
|
|
*
|
|
* \brief Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem
|
|
*
|
|
* \tparam _MatrixType the type of the matrix of which we are computing the
|
|
* eigendecomposition; this is expected to be an instantiation of the Matrix
|
|
* class template.
|
|
*
|
|
* 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
|
|
* a given matrix. Alternatively, you can use the
|
|
* GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
|
|
* constructor which computes the eigenvalues and eigenvectors at construction time.
|
|
* Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues()
|
|
* and eigenvectors() functions.
|
|
*
|
|
* The documentation for GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
|
|
* contains an example of the typical use of this class.
|
|
*
|
|
* \sa class SelfAdjointEigenSolver, class EigenSolver, class ComplexEigenSolver
|
|
*/
|
|
template<typename _MatrixType>
|
|
class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType>
|
|
{
|
|
typedef SelfAdjointEigenSolver<_MatrixType> Base;
|
|
public:
|
|
|
|
typedef typename Base::Index Index;
|
|
typedef _MatrixType MatrixType;
|
|
|
|
/** \brief Default constructor for fixed-size matrices.
|
|
*
|
|
* The default constructor is useful in cases in which the user intends to
|
|
* perform decompositions via compute(const MatrixType&, bool) or
|
|
* compute(const MatrixType&, const MatrixType&, bool). This constructor
|
|
* can only be used if \p _MatrixType is a fixed-size matrix; use
|
|
* SelfAdjointEigenSolver(Index) for dynamic-size matrices.
|
|
*
|
|
* Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp
|
|
* Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out
|
|
*/
|
|
GeneralizedSelfAdjointEigenSolver() : Base() {}
|
|
|
|
/** \brief Constructor, pre-allocates memory for dynamic-size matrices.
|
|
*
|
|
* \param [in] size Positive integer, size of the matrix whose
|
|
* eigenvalues and eigenvectors will be computed.
|
|
*
|
|
* This constructor is useful for dynamic-size matrices, when the user
|
|
* intends to perform decompositions via compute(const MatrixType&, bool)
|
|
* or compute(const MatrixType&, const MatrixType&, bool). The \p size
|
|
* parameter is only used as a hint. It is not an error to give a wrong
|
|
* \p size, but it may impair performance.
|
|
*
|
|
* \sa compute(const MatrixType&, bool) for an example
|
|
*/
|
|
GeneralizedSelfAdjointEigenSolver(Index size)
|
|
: Base(size)
|
|
{}
|
|
|
|
/** \brief Constructor; computes generalized eigendecomposition of given matrix pencil.
|
|
*
|
|
* \param[in] matA Selfadjoint matrix in matrix pencil.
|
|
* Only the lower triangular part of the matrix is referenced.
|
|
* \param[in] matB Positive-definite matrix in matrix pencil.
|
|
* 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.
|
|
*
|
|
* This constructor calls compute(const MatrixType&, const MatrixType&, int)
|
|
* to compute the eigenvalues and (if requested) the eigenvectors of the
|
|
* generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the
|
|
* selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix
|
|
* \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property
|
|
* \f$ x^* B x = 1 \f$. The eigenvectors are computed if
|
|
* \a options contains ComputeEigenvectors.
|
|
*
|
|
* 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$
|
|
*
|
|
* Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp
|
|
* Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out
|
|
*
|
|
* \sa compute(const MatrixType&, const MatrixType&, int)
|
|
*/
|
|
GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
|
|
int options = ComputeEigenvectors|Ax_lBx)
|
|
: Base(matA.cols())
|
|
{
|
|
compute(matA, matB, options);
|
|
}
|
|
|
|
/** \brief Computes generalized eigendecomposition of given matrix pencil.
|
|
*
|
|
* \param[in] matA Selfadjoint matrix in matrix pencil.
|
|
* Only the lower triangular part of the matrix is referenced.
|
|
* \param[in] matB Positive-definite matrix in matrix pencil.
|
|
* 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
|
|
* and (if requested) the eigenvectors of the generalized eigenproblem
|
|
* \f$ Ax = \lambda B x \f$ with \a matA the selfadjoint
|
|
* 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$
|
|
*
|
|
* The eigenvalues() function can be used to retrieve
|
|
* the eigenvalues. If \p options contains ComputeEigenvectors, then the
|
|
* eigenvectors are also computed and can be retrieved by calling
|
|
* eigenvectors().
|
|
*
|
|
* The implementation uses LLT to compute the Cholesky decomposition
|
|
* \f$ B = LL^* \f$ and calls compute(const MatrixType&, bool) to compute
|
|
* the eigendecomposition \f$ L^{-1} A (L^*)^{-1} \f$. This solves the
|
|
* generalized eigenproblem, because any solution of the generalized
|
|
* eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution
|
|
* \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the
|
|
* eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$.
|
|
*
|
|
* Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp
|
|
* Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out
|
|
*
|
|
* \sa SelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
|
|
*/
|
|
GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB,
|
|
int options = ComputeEigenvectors|Ax_lBx);
|
|
|
|
protected:
|
|
|
|
};
|
|
|
|
|
|
template<typename MatrixType>
|
|
GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
|
|
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)
|
|
&& "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;
|
|
|
|
// 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);
|
|
|
|
Base::compute(matC, options&EigVecMask);
|
|
|
|
// transform back the eigen vectors: evecs = inv(U) * evecs
|
|
if(computeEigVecs)
|
|
cholB.matrixU().solveInPlace(Base::m_eivec);
|
|
|
|
return *this;
|
|
}
|
|
|
|
#endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
|