Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h - mirror - Git at Google

 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.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
   *
   *
   * \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
       *
       * Accoring to \p options, this function computes eigenvalues and (if requested)
       * the eigenvectors of one of the following three generalized eigenproblems:
       * - \c Ax_lBx: \f$ Ax = \lambda B x \f$
       * - \c ABx_lx: \f$ ABx = \lambda x \f$
       * - \c BAx_lx: \f$ BAx = \lambda 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$.
       *
       * 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 computes the classical eigendecomposition
       * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx
       * and of \f$ L^{*} A L \f$ otherwise. 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$. Similar statements
       * can be made for the two other variants.
       *
       * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp
       * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out
       *
       * \sa GeneralizedSelfAdjointEigenSolver(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)
 {
   eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
   eigen_assert((options&~(EigVecMask|GenEigMask))==0
           && (options&EigVecMask)!=EigVecMask
           && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
            || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
           && "invalid option parameter");

   bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);

   // Compute the cholesky decomposition of matB = L L' = U'U
   LLT<MatrixType> cholB(matB);

   int type = (options&GenEigMask);
   if(type==0)
     type = Ax_lBx;

   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);

     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;
 }

 #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
	// This file is part of Eigen, a lightweight C++ template library
	// for linear algebra.
	//
	// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.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
	*
	*
	* \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
	*
	* Accoring to \p options, this function computes eigenvalues and (if requested)
	* the eigenvectors of one of the following three generalized eigenproblems:
	* - \c Ax_lBx: \f$ Ax = \lambda B x \f$
	* - \c ABx_lx: \f$ ABx = \lambda x \f$
	* - \c BAx_lx: \f$ BAx = \lambda 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$.
	*
	* 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 computes the classical eigendecomposition
	* of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx
	* and of \f$ L^{*} A L \f$ otherwise. 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$. Similar statements
	* can be made for the two other variants.
	*
	* Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp
	* Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out
	*
	* \sa GeneralizedSelfAdjointEigenSolver(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)
	{
	eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
	eigen_assert((options&~(EigVecMask\|GenEigMask))==0
	&& (options&EigVecMask)!=EigVecMask
	&& ((options&GenEigMask)==0 \|\| (options&GenEigMask)==Ax_lBx
	\|\| (options&GenEigMask)==ABx_lx \|\| (options&GenEigMask)==BAx_lx)
	&& "invalid option parameter");

	bool computeEigVecs = ((options&EigVecMask)==0) \|\| ((options&EigVecMask)==ComputeEigenvectors);

	// Compute the cholesky decomposition of matB = L L' = U'U
	LLT<MatrixType> cholB(matB);

	int type = (options&GenEigMask);
	if(type==0)
	type = Ax_lBx;

	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);

	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;
	}

	#endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H