1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
5// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
6//
7// This Source Code Form is subject to the terms of the Mozilla
8// Public License v. 2.0. If a copy of the MPL was not distributed
9// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11#ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
12#define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
13
14#include "./Tridiagonalization.h"
15
16namespace Eigen {
17
18/** \eigenvalues_module \ingroup Eigenvalues_Module
19  *
20  *
21  * \class GeneralizedSelfAdjointEigenSolver
22  *
23  * \brief Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem
24  *
25  * \tparam _MatrixType the type of the matrix of which we are computing the
26  * eigendecomposition; this is expected to be an instantiation of the Matrix
27  * class template.
28  *
29  * This class solves the generalized eigenvalue problem
30  * \f$ Av = \lambda Bv \f$. In this case, the matrix \f$ A \f$ should be
31  * selfadjoint and the matrix \f$ B \f$ should be positive definite.
32  *
33  * Only the \b lower \b triangular \b part of the input matrix is referenced.
34  *
35  * Call the function compute() to compute the eigenvalues and eigenvectors of
36  * a given matrix. Alternatively, you can use the
37  * GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
38  * constructor which computes the eigenvalues and eigenvectors at construction time.
39  * Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues()
40  * and eigenvectors() functions.
41  *
42  * The documentation for GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
43  * contains an example of the typical use of this class.
44  *
45  * \sa class SelfAdjointEigenSolver, class EigenSolver, class ComplexEigenSolver
46  */
47template<typename _MatrixType>
48class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType>
49{
50    typedef SelfAdjointEigenSolver<_MatrixType> Base;
51  public:
52
53    typedef typename Base::Index Index;
54    typedef _MatrixType MatrixType;
55
56    /** \brief Default constructor for fixed-size matrices.
57      *
58      * The default constructor is useful in cases in which the user intends to
59      * perform decompositions via compute(). This constructor
60      * can only be used if \p _MatrixType is a fixed-size matrix; use
61      * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices.
62      */
63    GeneralizedSelfAdjointEigenSolver() : Base() {}
64
65    /** \brief Constructor, pre-allocates memory for dynamic-size matrices.
66      *
67      * \param [in]  size  Positive integer, size of the matrix whose
68      * eigenvalues and eigenvectors will be computed.
69      *
70      * This constructor is useful for dynamic-size matrices, when the user
71      * intends to perform decompositions via compute(). The \p size
72      * parameter is only used as a hint. It is not an error to give a wrong
73      * \p size, but it may impair performance.
74      *
75      * \sa compute() for an example
76      */
77    GeneralizedSelfAdjointEigenSolver(Index size)
78        : Base(size)
79    {}
80
81    /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil.
82      *
83      * \param[in]  matA  Selfadjoint matrix in matrix pencil.
84      *                   Only the lower triangular part of the matrix is referenced.
85      * \param[in]  matB  Positive-definite matrix in matrix pencil.
86      *                   Only the lower triangular part of the matrix is referenced.
87      * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
88      *                     Default is #ComputeEigenvectors|#Ax_lBx.
89      *
90      * This constructor calls compute(const MatrixType&, const MatrixType&, int)
91      * to compute the eigenvalues and (if requested) the eigenvectors of the
92      * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the
93      * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix
94      * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property
95      * \f$ x^* B x = 1 \f$. The eigenvectors are computed if
96      * \a options contains ComputeEigenvectors.
97      *
98      * In addition, the two following variants can be solved via \p options:
99      * - \c ABx_lx: \f$ ABx = \lambda x \f$
100      * - \c BAx_lx: \f$ BAx = \lambda x \f$
101      *
102      * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp
103      * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out
104      *
105      * \sa compute(const MatrixType&, const MatrixType&, int)
106      */
107    GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
108                                      int options = ComputeEigenvectors|Ax_lBx)
109      : Base(matA.cols())
110    {
111      compute(matA, matB, options);
112    }
113
114    /** \brief Computes generalized eigendecomposition of given matrix pencil.
115      *
116      * \param[in]  matA  Selfadjoint matrix in matrix pencil.
117      *                   Only the lower triangular part of the matrix is referenced.
118      * \param[in]  matB  Positive-definite matrix in matrix pencil.
119      *                   Only the lower triangular part of the matrix is referenced.
120      * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
121      *                     Default is #ComputeEigenvectors|#Ax_lBx.
122      *
123      * \returns    Reference to \c *this
124      *
125      * Accoring to \p options, this function computes eigenvalues and (if requested)
126      * the eigenvectors of one of the following three generalized eigenproblems:
127      * - \c Ax_lBx: \f$ Ax = \lambda B x \f$
128      * - \c ABx_lx: \f$ ABx = \lambda x \f$
129      * - \c BAx_lx: \f$ BAx = \lambda x \f$
130      * with \a matA the selfadjoint matrix \f$ A \f$ and \a matB the positive definite
131      * matrix \f$ B \f$.
132      * In addition, each eigenvector \f$ x \f$ satisfies the property \f$ x^* B x = 1 \f$.
133      *
134      * The eigenvalues() function can be used to retrieve
135      * the eigenvalues. If \p options contains ComputeEigenvectors, then the
136      * eigenvectors are also computed and can be retrieved by calling
137      * eigenvectors().
138      *
139      * The implementation uses LLT to compute the Cholesky decomposition
140      * \f$ B = LL^* \f$ and computes the classical eigendecomposition
141      * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx
142      * and of \f$ L^{*} A L \f$ otherwise. This solves the
143      * generalized eigenproblem, because any solution of the generalized
144      * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution
145      * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the
146      * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. Similar statements
147      * can be made for the two other variants.
148      *
149      * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp
150      * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out
151      *
152      * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
153      */
154    GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB,
155                                               int options = ComputeEigenvectors|Ax_lBx);
156
157  protected:
158
159};
160
161
162template<typename MatrixType>
163GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
164compute(const MatrixType& matA, const MatrixType& matB, int options)
165{
166  eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
167  eigen_assert((options&~(EigVecMask|GenEigMask))==0
168          && (options&EigVecMask)!=EigVecMask
169          && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
170           || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
171          && "invalid option parameter");
172
173  bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
174
175  // Compute the cholesky decomposition of matB = L L' = U'U
176  LLT<MatrixType> cholB(matB);
177
178  int type = (options&GenEigMask);
179  if(type==0)
180    type = Ax_lBx;
181
182  if(type==Ax_lBx)
183  {
184    // compute C = inv(L) A inv(L')
185    MatrixType matC = matA.template selfadjointView<Lower>();
186    cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
187    cholB.matrixU().template solveInPlace<OnTheRight>(matC);
188
189    Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
190
191    // transform back the eigen vectors: evecs = inv(U) * evecs
192    if(computeEigVecs)
193      cholB.matrixU().solveInPlace(Base::m_eivec);
194  }
195  else if(type==ABx_lx)
196  {
197    // compute C = L' A L
198    MatrixType matC = matA.template selfadjointView<Lower>();
199    matC = matC * cholB.matrixL();
200    matC = cholB.matrixU() * matC;
201
202    Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
203
204    // transform back the eigen vectors: evecs = inv(U) * evecs
205    if(computeEigVecs)
206      cholB.matrixU().solveInPlace(Base::m_eivec);
207  }
208  else if(type==BAx_lx)
209  {
210    // compute C = L' A L
211    MatrixType matC = matA.template selfadjointView<Lower>();
212    matC = matC * cholB.matrixL();
213    matC = cholB.matrixU() * matC;
214
215    Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
216
217    // transform back the eigen vectors: evecs = L * evecs
218    if(computeEigVecs)
219      Base::m_eivec = cholB.matrixL() * Base::m_eivec;
220  }
221
222  return *this;
223}
224
225} // end namespace Eigen
226
227#endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
228