diff options
author | Gael Guennebaud <g.gael@free.fr> | 2010-06-16 23:48:16 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2010-06-16 23:48:16 +0200 |
commit | 74006a9fe917584974f29295f40e81b2c64fc6f4 (patch) | |
tree | d88547f79eaeaa0b782605b362d4ff6ad2ebf540 /Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h | |
parent | 197ce96c0087e48201e5bda525ab67163dc9b661 (diff) |
* decouple the generalized selfadjoint eigenvalue problem to the standard one
* uses named values instead of bools
Diffstat (limited to 'Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h')
-rw-r--r-- | Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h | 210 |
1 files changed, 210 insertions, 0 deletions
diff --git a/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h new file mode 100644 index 000000000..66d32dee6 --- /dev/null +++ b/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h @@ -0,0 +1,210 @@ +// 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"); + + 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 |