Michael Ernst Peter / Eigen

Eigen libary for mbed

src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h@1:3b8049da21b8, 2019-09-24 (annotated)

Committer:
jsoh91
Date:
Tue Sep 24 00:18:23 2019 +0000
Revision:
1:3b8049da21b8
Parent:
0:13a5d365ba16 
ignore and revise some of error parts

Who changed what in which revision?

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