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