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 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 */
GeneralizedSelfAdjointEigenSolver()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 */
GeneralizedSelfAdjointEigenSolver(Index size)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
162 template<typename MatrixType>
163 GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
compute(const MatrixType & matA,const MatrixType & matB,int options)164 compute(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