• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
16 namespace 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   */
47 template<typename _MatrixType>
48 class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType>
49 {
50     typedef SelfAdjointEigenSolver<_MatrixType> Base;
51   public:
52 
53     typedef _MatrixType MatrixType;
54 
55     /** \brief Default constructor for fixed-size matrices.
56       *
57       * The default constructor is useful in cases in which the user intends to
58       * perform decompositions via compute(). This constructor
59       * can only be used if \p _MatrixType is a fixed-size matrix; use
60       * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices.
61       */
GeneralizedSelfAdjointEigenSolver()62     GeneralizedSelfAdjointEigenSolver() : Base() {}
63 
64     /** \brief Constructor, pre-allocates memory for dynamic-size matrices.
65       *
66       * \param [in]  size  Positive integer, size of the matrix whose
67       * eigenvalues and eigenvectors will be computed.
68       *
69       * This constructor is useful for dynamic-size matrices, when the user
70       * intends to perform decompositions via compute(). The \p size
71       * parameter is only used as a hint. It is not an error to give a wrong
72       * \p size, but it may impair performance.
73       *
74       * \sa compute() for an example
75       */
GeneralizedSelfAdjointEigenSolver(Index size)76     explicit GeneralizedSelfAdjointEigenSolver(Index size)
77         : Base(size)
78     {}
79 
80     /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil.
81       *
82       * \param[in]  matA  Selfadjoint matrix in matrix pencil.
83       *                   Only the lower triangular part of the matrix is referenced.
84       * \param[in]  matB  Positive-definite matrix in matrix pencil.
85       *                   Only the lower triangular part of the matrix is referenced.
86       * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
87       *                     Default is #ComputeEigenvectors|#Ax_lBx.
88       *
89       * This constructor calls compute(const MatrixType&, const MatrixType&, int)
90       * to compute the eigenvalues and (if requested) the eigenvectors of the
91       * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the
92       * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix
93       * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property
94       * \f$ x^* B x = 1 \f$. The eigenvectors are computed if
95       * \a options contains ComputeEigenvectors.
96       *
97       * In addition, the two following variants can be solved via \p options:
98       * - \c ABx_lx: \f$ ABx = \lambda x \f$
99       * - \c BAx_lx: \f$ BAx = \lambda x \f$
100       *
101       * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp
102       * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out
103       *
104       * \sa compute(const MatrixType&, const MatrixType&, int)
105       */
106     GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
107                                       int options = ComputeEigenvectors|Ax_lBx)
108       : Base(matA.cols())
109     {
110       compute(matA, matB, options);
111     }
112 
113     /** \brief Computes generalized eigendecomposition of given matrix pencil.
114       *
115       * \param[in]  matA  Selfadjoint matrix in matrix pencil.
116       *                   Only the lower triangular part of the matrix is referenced.
117       * \param[in]  matB  Positive-definite matrix in matrix pencil.
118       *                   Only the lower triangular part of the matrix is referenced.
119       * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
120       *                     Default is #ComputeEigenvectors|#Ax_lBx.
121       *
122       * \returns    Reference to \c *this
123       *
124       * Accoring to \p options, this function computes eigenvalues and (if requested)
125       * the eigenvectors of one of the following three generalized eigenproblems:
126       * - \c Ax_lBx: \f$ Ax = \lambda B x \f$
127       * - \c ABx_lx: \f$ ABx = \lambda x \f$
128       * - \c BAx_lx: \f$ BAx = \lambda x \f$
129       * with \a matA the selfadjoint matrix \f$ A \f$ and \a matB the positive definite
130       * matrix \f$ B \f$.
131       * In addition, each eigenvector \f$ x \f$ satisfies the property \f$ x^* B x = 1 \f$.
132       *
133       * The eigenvalues() function can be used to retrieve
134       * the eigenvalues. If \p options contains ComputeEigenvectors, then the
135       * eigenvectors are also computed and can be retrieved by calling
136       * eigenvectors().
137       *
138       * The implementation uses LLT to compute the Cholesky decomposition
139       * \f$ B = LL^* \f$ and computes the classical eigendecomposition
140       * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx
141       * and of \f$ L^{*} A L \f$ otherwise. This solves the
142       * generalized eigenproblem, because any solution of the generalized
143       * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution
144       * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the
145       * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. Similar statements
146       * can be made for the two other variants.
147       *
148       * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp
149       * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out
150       *
151       * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
152       */
153     GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB,
154                                                int options = ComputeEigenvectors|Ax_lBx);
155 
156   protected:
157 
158 };
159 
160 
161 template<typename MatrixType>
162 GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
compute(const MatrixType & matA,const MatrixType & matB,int options)163 compute(const MatrixType& matA, const MatrixType& matB, int options)
164 {
165   eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
166   eigen_assert((options&~(EigVecMask|GenEigMask))==0
167           && (options&EigVecMask)!=EigVecMask
168           && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
169            || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
170           && "invalid option parameter");
171 
172   bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
173 
174   // Compute the cholesky decomposition of matB = L L' = U'U
175   LLT<MatrixType> cholB(matB);
176 
177   int type = (options&GenEigMask);
178   if(type==0)
179     type = Ax_lBx;
180 
181   if(type==Ax_lBx)
182   {
183     // compute C = inv(L) A inv(L')
184     MatrixType matC = matA.template selfadjointView<Lower>();
185     cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
186     cholB.matrixU().template solveInPlace<OnTheRight>(matC);
187 
188     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
189 
190     // transform back the eigen vectors: evecs = inv(U) * evecs
191     if(computeEigVecs)
192       cholB.matrixU().solveInPlace(Base::m_eivec);
193   }
194   else if(type==ABx_lx)
195   {
196     // compute C = L' A L
197     MatrixType matC = matA.template selfadjointView<Lower>();
198     matC = matC * cholB.matrixL();
199     matC = cholB.matrixU() * matC;
200 
201     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
202 
203     // transform back the eigen vectors: evecs = inv(U) * evecs
204     if(computeEigVecs)
205       cholB.matrixU().solveInPlace(Base::m_eivec);
206   }
207   else if(type==BAx_lx)
208   {
209     // compute C = L' A L
210     MatrixType matC = matA.template selfadjointView<Lower>();
211     matC = matC * cholB.matrixL();
212     matC = cholB.matrixU() * matC;
213 
214     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
215 
216     // transform back the eigen vectors: evecs = L * evecs
217     if(computeEigVecs)
218       Base::m_eivec = cholB.matrixL() * Base::m_eivec;
219   }
220 
221   return *this;
222 }
223 
224 } // end namespace Eigen
225 
226 #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
227