Eigne Matrix Class Library
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
Eigen Matrix Class Library for mbed.
Finally, you can use Eigen on your mbed!!!
src/Eigenvalues/GeneralizedEigenSolver.h@0:13a5d365ba16, 2016-10-13 (annotated)
- Committer:
- ykuroda
- Date:
- Thu Oct 13 04:07:23 2016 +0000
- Revision:
- 0:13a5d365ba16
First commint, Eigne Matrix Class Library
Who changed what in which revision?
User | Revision | Line number | New 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) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> |
ykuroda | 0:13a5d365ba16 | 5 | // Copyright (C) 2010,2012 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_GENERALIZEDEIGENSOLVER_H |
ykuroda | 0:13a5d365ba16 | 12 | #define EIGEN_GENERALIZEDEIGENSOLVER_H |
ykuroda | 0:13a5d365ba16 | 13 | |
ykuroda | 0:13a5d365ba16 | 14 | #include "./RealQZ.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 GeneralizedEigenSolver |
ykuroda | 0:13a5d365ba16 | 22 | * |
ykuroda | 0:13a5d365ba16 | 23 | * \brief Computes the generalized eigenvalues and eigenvectors of a pair of general matrices |
ykuroda | 0:13a5d365ba16 | 24 | * |
ykuroda | 0:13a5d365ba16 | 25 | * \tparam _MatrixType the type of the matrices of which we are computing the |
ykuroda | 0:13a5d365ba16 | 26 | * eigen-decomposition; this is expected to be an instantiation of the Matrix |
ykuroda | 0:13a5d365ba16 | 27 | * class template. Currently, only real matrices are supported. |
ykuroda | 0:13a5d365ba16 | 28 | * |
ykuroda | 0:13a5d365ba16 | 29 | * The generalized eigenvalues and eigenvectors of a matrix pair \f$ A \f$ and \f$ B \f$ are scalars |
ykuroda | 0:13a5d365ba16 | 30 | * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda Bv \f$. If |
ykuroda | 0:13a5d365ba16 | 31 | * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and |
ykuroda | 0:13a5d365ba16 | 32 | * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V = |
ykuroda | 0:13a5d365ba16 | 33 | * B V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we |
ykuroda | 0:13a5d365ba16 | 34 | * have \f$ A = B V D V^{-1} \f$. This is called the generalized eigen-decomposition. |
ykuroda | 0:13a5d365ba16 | 35 | * |
ykuroda | 0:13a5d365ba16 | 36 | * The generalized eigenvalues and eigenvectors of a matrix pair may be complex, even when the |
ykuroda | 0:13a5d365ba16 | 37 | * matrices are real. Moreover, the generalized eigenvalue might be infinite if the matrix B is |
ykuroda | 0:13a5d365ba16 | 38 | * singular. To workaround this difficulty, the eigenvalues are provided as a pair of complex \f$ \alpha \f$ |
ykuroda | 0:13a5d365ba16 | 39 | * and real \f$ \beta \f$ such that: \f$ \lambda_i = \alpha_i / \beta_i \f$. If \f$ \beta_i \f$ is (nearly) zero, |
ykuroda | 0:13a5d365ba16 | 40 | * then one can consider the well defined left eigenvalue \f$ \mu = \beta_i / \alpha_i\f$ such that: |
ykuroda | 0:13a5d365ba16 | 41 | * \f$ \mu_i A v_i = B v_i \f$, or even \f$ \mu_i u_i^T A = u_i^T B \f$ where \f$ u_i \f$ is |
ykuroda | 0:13a5d365ba16 | 42 | * called the left eigenvector. |
ykuroda | 0:13a5d365ba16 | 43 | * |
ykuroda | 0:13a5d365ba16 | 44 | * Call the function compute() to compute the generalized eigenvalues and eigenvectors of |
ykuroda | 0:13a5d365ba16 | 45 | * a given matrix pair. Alternatively, you can use the |
ykuroda | 0:13a5d365ba16 | 46 | * GeneralizedEigenSolver(const MatrixType&, const MatrixType&, bool) constructor which computes the |
ykuroda | 0:13a5d365ba16 | 47 | * eigenvalues and eigenvectors at construction time. Once the eigenvalue and |
ykuroda | 0:13a5d365ba16 | 48 | * eigenvectors are computed, they can be retrieved with the eigenvalues() and |
ykuroda | 0:13a5d365ba16 | 49 | * eigenvectors() functions. |
ykuroda | 0:13a5d365ba16 | 50 | * |
ykuroda | 0:13a5d365ba16 | 51 | * Here is an usage example of this class: |
ykuroda | 0:13a5d365ba16 | 52 | * Example: \include GeneralizedEigenSolver.cpp |
ykuroda | 0:13a5d365ba16 | 53 | * Output: \verbinclude GeneralizedEigenSolver.out |
ykuroda | 0:13a5d365ba16 | 54 | * |
ykuroda | 0:13a5d365ba16 | 55 | * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver |
ykuroda | 0:13a5d365ba16 | 56 | */ |
ykuroda | 0:13a5d365ba16 | 57 | template<typename _MatrixType> class GeneralizedEigenSolver |
ykuroda | 0:13a5d365ba16 | 58 | { |
ykuroda | 0:13a5d365ba16 | 59 | public: |
ykuroda | 0:13a5d365ba16 | 60 | |
ykuroda | 0:13a5d365ba16 | 61 | /** \brief Synonym for the template parameter \p _MatrixType. */ |
ykuroda | 0:13a5d365ba16 | 62 | typedef _MatrixType MatrixType; |
ykuroda | 0:13a5d365ba16 | 63 | |
ykuroda | 0:13a5d365ba16 | 64 | enum { |
ykuroda | 0:13a5d365ba16 | 65 | RowsAtCompileTime = MatrixType::RowsAtCompileTime, |
ykuroda | 0:13a5d365ba16 | 66 | ColsAtCompileTime = MatrixType::ColsAtCompileTime, |
ykuroda | 0:13a5d365ba16 | 67 | Options = MatrixType::Options, |
ykuroda | 0:13a5d365ba16 | 68 | MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, |
ykuroda | 0:13a5d365ba16 | 69 | MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime |
ykuroda | 0:13a5d365ba16 | 70 | }; |
ykuroda | 0:13a5d365ba16 | 71 | |
ykuroda | 0:13a5d365ba16 | 72 | /** \brief Scalar type for matrices of type #MatrixType. */ |
ykuroda | 0:13a5d365ba16 | 73 | typedef typename MatrixType::Scalar Scalar; |
ykuroda | 0:13a5d365ba16 | 74 | typedef typename NumTraits<Scalar>::Real RealScalar; |
ykuroda | 0:13a5d365ba16 | 75 | typedef typename MatrixType::Index Index; |
ykuroda | 0:13a5d365ba16 | 76 | |
ykuroda | 0:13a5d365ba16 | 77 | /** \brief Complex scalar type for #MatrixType. |
ykuroda | 0:13a5d365ba16 | 78 | * |
ykuroda | 0:13a5d365ba16 | 79 | * This is \c std::complex<Scalar> if #Scalar is real (e.g., |
ykuroda | 0:13a5d365ba16 | 80 | * \c float or \c double) and just \c Scalar if #Scalar is |
ykuroda | 0:13a5d365ba16 | 81 | * complex. |
ykuroda | 0:13a5d365ba16 | 82 | */ |
ykuroda | 0:13a5d365ba16 | 83 | typedef std::complex<RealScalar> ComplexScalar; |
ykuroda | 0:13a5d365ba16 | 84 | |
ykuroda | 0:13a5d365ba16 | 85 | /** \brief Type for vector of real scalar values eigenvalues as returned by betas(). |
ykuroda | 0:13a5d365ba16 | 86 | * |
ykuroda | 0:13a5d365ba16 | 87 | * This is a column vector with entries of type #Scalar. |
ykuroda | 0:13a5d365ba16 | 88 | * The length of the vector is the size of #MatrixType. |
ykuroda | 0:13a5d365ba16 | 89 | */ |
ykuroda | 0:13a5d365ba16 | 90 | typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType; |
ykuroda | 0:13a5d365ba16 | 91 | |
ykuroda | 0:13a5d365ba16 | 92 | /** \brief Type for vector of complex scalar values eigenvalues as returned by betas(). |
ykuroda | 0:13a5d365ba16 | 93 | * |
ykuroda | 0:13a5d365ba16 | 94 | * This is a column vector with entries of type #ComplexScalar. |
ykuroda | 0:13a5d365ba16 | 95 | * The length of the vector is the size of #MatrixType. |
ykuroda | 0:13a5d365ba16 | 96 | */ |
ykuroda | 0:13a5d365ba16 | 97 | typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ComplexVectorType; |
ykuroda | 0:13a5d365ba16 | 98 | |
ykuroda | 0:13a5d365ba16 | 99 | /** \brief Expression type for the eigenvalues as returned by eigenvalues(). |
ykuroda | 0:13a5d365ba16 | 100 | */ |
ykuroda | 0:13a5d365ba16 | 101 | typedef CwiseBinaryOp<internal::scalar_quotient_op<ComplexScalar,Scalar>,ComplexVectorType,VectorType> EigenvalueType; |
ykuroda | 0:13a5d365ba16 | 102 | |
ykuroda | 0:13a5d365ba16 | 103 | /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). |
ykuroda | 0:13a5d365ba16 | 104 | * |
ykuroda | 0:13a5d365ba16 | 105 | * This is a square matrix with entries of type #ComplexScalar. |
ykuroda | 0:13a5d365ba16 | 106 | * The size is the same as the size of #MatrixType. |
ykuroda | 0:13a5d365ba16 | 107 | */ |
ykuroda | 0:13a5d365ba16 | 108 | typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType; |
ykuroda | 0:13a5d365ba16 | 109 | |
ykuroda | 0:13a5d365ba16 | 110 | /** \brief Default constructor. |
ykuroda | 0:13a5d365ba16 | 111 | * |
ykuroda | 0:13a5d365ba16 | 112 | * The default constructor is useful in cases in which the user intends to |
ykuroda | 0:13a5d365ba16 | 113 | * perform decompositions via EigenSolver::compute(const MatrixType&, bool). |
ykuroda | 0:13a5d365ba16 | 114 | * |
ykuroda | 0:13a5d365ba16 | 115 | * \sa compute() for an example. |
ykuroda | 0:13a5d365ba16 | 116 | */ |
ykuroda | 0:13a5d365ba16 | 117 | GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {} |
ykuroda | 0:13a5d365ba16 | 118 | |
ykuroda | 0:13a5d365ba16 | 119 | /** \brief Default constructor with memory preallocation |
ykuroda | 0:13a5d365ba16 | 120 | * |
ykuroda | 0:13a5d365ba16 | 121 | * Like the default constructor but with preallocation of the internal data |
ykuroda | 0:13a5d365ba16 | 122 | * according to the specified problem \a size. |
ykuroda | 0:13a5d365ba16 | 123 | * \sa GeneralizedEigenSolver() |
ykuroda | 0:13a5d365ba16 | 124 | */ |
ykuroda | 0:13a5d365ba16 | 125 | GeneralizedEigenSolver(Index size) |
ykuroda | 0:13a5d365ba16 | 126 | : m_eivec(size, size), |
ykuroda | 0:13a5d365ba16 | 127 | m_alphas(size), |
ykuroda | 0:13a5d365ba16 | 128 | m_betas(size), |
ykuroda | 0:13a5d365ba16 | 129 | m_isInitialized(false), |
ykuroda | 0:13a5d365ba16 | 130 | m_eigenvectorsOk(false), |
ykuroda | 0:13a5d365ba16 | 131 | m_realQZ(size), |
ykuroda | 0:13a5d365ba16 | 132 | m_matS(size, size), |
ykuroda | 0:13a5d365ba16 | 133 | m_tmp(size) |
ykuroda | 0:13a5d365ba16 | 134 | {} |
ykuroda | 0:13a5d365ba16 | 135 | |
ykuroda | 0:13a5d365ba16 | 136 | /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair. |
ykuroda | 0:13a5d365ba16 | 137 | * |
ykuroda | 0:13a5d365ba16 | 138 | * \param[in] A Square matrix whose eigendecomposition is to be computed. |
ykuroda | 0:13a5d365ba16 | 139 | * \param[in] B Square matrix whose eigendecomposition is to be computed. |
ykuroda | 0:13a5d365ba16 | 140 | * \param[in] computeEigenvectors If true, both the eigenvectors and the |
ykuroda | 0:13a5d365ba16 | 141 | * eigenvalues are computed; if false, only the eigenvalues are computed. |
ykuroda | 0:13a5d365ba16 | 142 | * |
ykuroda | 0:13a5d365ba16 | 143 | * This constructor calls compute() to compute the generalized eigenvalues |
ykuroda | 0:13a5d365ba16 | 144 | * and eigenvectors. |
ykuroda | 0:13a5d365ba16 | 145 | * |
ykuroda | 0:13a5d365ba16 | 146 | * \sa compute() |
ykuroda | 0:13a5d365ba16 | 147 | */ |
ykuroda | 0:13a5d365ba16 | 148 | GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true) |
ykuroda | 0:13a5d365ba16 | 149 | : m_eivec(A.rows(), A.cols()), |
ykuroda | 0:13a5d365ba16 | 150 | m_alphas(A.cols()), |
ykuroda | 0:13a5d365ba16 | 151 | m_betas(A.cols()), |
ykuroda | 0:13a5d365ba16 | 152 | m_isInitialized(false), |
ykuroda | 0:13a5d365ba16 | 153 | m_eigenvectorsOk(false), |
ykuroda | 0:13a5d365ba16 | 154 | m_realQZ(A.cols()), |
ykuroda | 0:13a5d365ba16 | 155 | m_matS(A.rows(), A.cols()), |
ykuroda | 0:13a5d365ba16 | 156 | m_tmp(A.cols()) |
ykuroda | 0:13a5d365ba16 | 157 | { |
ykuroda | 0:13a5d365ba16 | 158 | compute(A, B, computeEigenvectors); |
ykuroda | 0:13a5d365ba16 | 159 | } |
ykuroda | 0:13a5d365ba16 | 160 | |
ykuroda | 0:13a5d365ba16 | 161 | /* \brief Returns the computed generalized eigenvectors. |
ykuroda | 0:13a5d365ba16 | 162 | * |
ykuroda | 0:13a5d365ba16 | 163 | * \returns %Matrix whose columns are the (possibly complex) eigenvectors. |
ykuroda | 0:13a5d365ba16 | 164 | * |
ykuroda | 0:13a5d365ba16 | 165 | * \pre Either the constructor |
ykuroda | 0:13a5d365ba16 | 166 | * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function |
ykuroda | 0:13a5d365ba16 | 167 | * compute(const MatrixType&, const MatrixType& bool) has been called before, and |
ykuroda | 0:13a5d365ba16 | 168 | * \p computeEigenvectors was set to true (the default). |
ykuroda | 0:13a5d365ba16 | 169 | * |
ykuroda | 0:13a5d365ba16 | 170 | * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding |
ykuroda | 0:13a5d365ba16 | 171 | * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The |
ykuroda | 0:13a5d365ba16 | 172 | * eigenvectors are normalized to have (Euclidean) norm equal to one. The |
ykuroda | 0:13a5d365ba16 | 173 | * matrix returned by this function is the matrix \f$ V \f$ in the |
ykuroda | 0:13a5d365ba16 | 174 | * generalized eigendecomposition \f$ A = B V D V^{-1} \f$, if it exists. |
ykuroda | 0:13a5d365ba16 | 175 | * |
ykuroda | 0:13a5d365ba16 | 176 | * \sa eigenvalues() |
ykuroda | 0:13a5d365ba16 | 177 | */ |
ykuroda | 0:13a5d365ba16 | 178 | // EigenvectorsType eigenvectors() const; |
ykuroda | 0:13a5d365ba16 | 179 | |
ykuroda | 0:13a5d365ba16 | 180 | /** \brief Returns an expression of the computed generalized eigenvalues. |
ykuroda | 0:13a5d365ba16 | 181 | * |
ykuroda | 0:13a5d365ba16 | 182 | * \returns An expression of the column vector containing the eigenvalues. |
ykuroda | 0:13a5d365ba16 | 183 | * |
ykuroda | 0:13a5d365ba16 | 184 | * It is a shortcut for \code this->alphas().cwiseQuotient(this->betas()); \endcode |
ykuroda | 0:13a5d365ba16 | 185 | * Not that betas might contain zeros. It is therefore not recommended to use this function, |
ykuroda | 0:13a5d365ba16 | 186 | * but rather directly deal with the alphas and betas vectors. |
ykuroda | 0:13a5d365ba16 | 187 | * |
ykuroda | 0:13a5d365ba16 | 188 | * \pre Either the constructor |
ykuroda | 0:13a5d365ba16 | 189 | * GeneralizedEigenSolver(const MatrixType&,const MatrixType&,bool) or the member function |
ykuroda | 0:13a5d365ba16 | 190 | * compute(const MatrixType&,const MatrixType&,bool) has been called before. |
ykuroda | 0:13a5d365ba16 | 191 | * |
ykuroda | 0:13a5d365ba16 | 192 | * The eigenvalues are repeated according to their algebraic multiplicity, |
ykuroda | 0:13a5d365ba16 | 193 | * so there are as many eigenvalues as rows in the matrix. The eigenvalues |
ykuroda | 0:13a5d365ba16 | 194 | * are not sorted in any particular order. |
ykuroda | 0:13a5d365ba16 | 195 | * |
ykuroda | 0:13a5d365ba16 | 196 | * \sa alphas(), betas(), eigenvectors() |
ykuroda | 0:13a5d365ba16 | 197 | */ |
ykuroda | 0:13a5d365ba16 | 198 | EigenvalueType eigenvalues() const |
ykuroda | 0:13a5d365ba16 | 199 | { |
ykuroda | 0:13a5d365ba16 | 200 | eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); |
ykuroda | 0:13a5d365ba16 | 201 | return EigenvalueType(m_alphas,m_betas); |
ykuroda | 0:13a5d365ba16 | 202 | } |
ykuroda | 0:13a5d365ba16 | 203 | |
ykuroda | 0:13a5d365ba16 | 204 | /** \returns A const reference to the vectors containing the alpha values |
ykuroda | 0:13a5d365ba16 | 205 | * |
ykuroda | 0:13a5d365ba16 | 206 | * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j). |
ykuroda | 0:13a5d365ba16 | 207 | * |
ykuroda | 0:13a5d365ba16 | 208 | * \sa betas(), eigenvalues() */ |
ykuroda | 0:13a5d365ba16 | 209 | ComplexVectorType alphas() const |
ykuroda | 0:13a5d365ba16 | 210 | { |
ykuroda | 0:13a5d365ba16 | 211 | eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); |
ykuroda | 0:13a5d365ba16 | 212 | return m_alphas; |
ykuroda | 0:13a5d365ba16 | 213 | } |
ykuroda | 0:13a5d365ba16 | 214 | |
ykuroda | 0:13a5d365ba16 | 215 | /** \returns A const reference to the vectors containing the beta values |
ykuroda | 0:13a5d365ba16 | 216 | * |
ykuroda | 0:13a5d365ba16 | 217 | * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j). |
ykuroda | 0:13a5d365ba16 | 218 | * |
ykuroda | 0:13a5d365ba16 | 219 | * \sa alphas(), eigenvalues() */ |
ykuroda | 0:13a5d365ba16 | 220 | VectorType betas() const |
ykuroda | 0:13a5d365ba16 | 221 | { |
ykuroda | 0:13a5d365ba16 | 222 | eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); |
ykuroda | 0:13a5d365ba16 | 223 | return m_betas; |
ykuroda | 0:13a5d365ba16 | 224 | } |
ykuroda | 0:13a5d365ba16 | 225 | |
ykuroda | 0:13a5d365ba16 | 226 | /** \brief Computes generalized eigendecomposition of given matrix. |
ykuroda | 0:13a5d365ba16 | 227 | * |
ykuroda | 0:13a5d365ba16 | 228 | * \param[in] A Square matrix whose eigendecomposition is to be computed. |
ykuroda | 0:13a5d365ba16 | 229 | * \param[in] B Square matrix whose eigendecomposition is to be computed. |
ykuroda | 0:13a5d365ba16 | 230 | * \param[in] computeEigenvectors If true, both the eigenvectors and the |
ykuroda | 0:13a5d365ba16 | 231 | * eigenvalues are computed; if false, only the eigenvalues are |
ykuroda | 0:13a5d365ba16 | 232 | * computed. |
ykuroda | 0:13a5d365ba16 | 233 | * \returns Reference to \c *this |
ykuroda | 0:13a5d365ba16 | 234 | * |
ykuroda | 0:13a5d365ba16 | 235 | * This function computes the eigenvalues of the real matrix \p matrix. |
ykuroda | 0:13a5d365ba16 | 236 | * The eigenvalues() function can be used to retrieve them. If |
ykuroda | 0:13a5d365ba16 | 237 | * \p computeEigenvectors is true, then the eigenvectors are also computed |
ykuroda | 0:13a5d365ba16 | 238 | * and can be retrieved by calling eigenvectors(). |
ykuroda | 0:13a5d365ba16 | 239 | * |
ykuroda | 0:13a5d365ba16 | 240 | * The matrix is first reduced to real generalized Schur form using the RealQZ |
ykuroda | 0:13a5d365ba16 | 241 | * class. The generalized Schur decomposition is then used to compute the eigenvalues |
ykuroda | 0:13a5d365ba16 | 242 | * and eigenvectors. |
ykuroda | 0:13a5d365ba16 | 243 | * |
ykuroda | 0:13a5d365ba16 | 244 | * The cost of the computation is dominated by the cost of the |
ykuroda | 0:13a5d365ba16 | 245 | * generalized Schur decomposition. |
ykuroda | 0:13a5d365ba16 | 246 | * |
ykuroda | 0:13a5d365ba16 | 247 | * This method reuses of the allocated data in the GeneralizedEigenSolver object. |
ykuroda | 0:13a5d365ba16 | 248 | */ |
ykuroda | 0:13a5d365ba16 | 249 | GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true); |
ykuroda | 0:13a5d365ba16 | 250 | |
ykuroda | 0:13a5d365ba16 | 251 | ComputationInfo info() const |
ykuroda | 0:13a5d365ba16 | 252 | { |
ykuroda | 0:13a5d365ba16 | 253 | eigen_assert(m_isInitialized && "EigenSolver is not initialized."); |
ykuroda | 0:13a5d365ba16 | 254 | return m_realQZ.info(); |
ykuroda | 0:13a5d365ba16 | 255 | } |
ykuroda | 0:13a5d365ba16 | 256 | |
ykuroda | 0:13a5d365ba16 | 257 | /** Sets the maximal number of iterations allowed. |
ykuroda | 0:13a5d365ba16 | 258 | */ |
ykuroda | 0:13a5d365ba16 | 259 | GeneralizedEigenSolver& setMaxIterations(Index maxIters) |
ykuroda | 0:13a5d365ba16 | 260 | { |
ykuroda | 0:13a5d365ba16 | 261 | m_realQZ.setMaxIterations(maxIters); |
ykuroda | 0:13a5d365ba16 | 262 | return *this; |
ykuroda | 0:13a5d365ba16 | 263 | } |
ykuroda | 0:13a5d365ba16 | 264 | |
ykuroda | 0:13a5d365ba16 | 265 | protected: |
ykuroda | 0:13a5d365ba16 | 266 | |
ykuroda | 0:13a5d365ba16 | 267 | static void check_template_parameters() |
ykuroda | 0:13a5d365ba16 | 268 | { |
ykuroda | 0:13a5d365ba16 | 269 | EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); |
ykuroda | 0:13a5d365ba16 | 270 | EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL); |
ykuroda | 0:13a5d365ba16 | 271 | } |
ykuroda | 0:13a5d365ba16 | 272 | |
ykuroda | 0:13a5d365ba16 | 273 | MatrixType m_eivec; |
ykuroda | 0:13a5d365ba16 | 274 | ComplexVectorType m_alphas; |
ykuroda | 0:13a5d365ba16 | 275 | VectorType m_betas; |
ykuroda | 0:13a5d365ba16 | 276 | bool m_isInitialized; |
ykuroda | 0:13a5d365ba16 | 277 | bool m_eigenvectorsOk; |
ykuroda | 0:13a5d365ba16 | 278 | RealQZ<MatrixType> m_realQZ; |
ykuroda | 0:13a5d365ba16 | 279 | MatrixType m_matS; |
ykuroda | 0:13a5d365ba16 | 280 | |
ykuroda | 0:13a5d365ba16 | 281 | typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; |
ykuroda | 0:13a5d365ba16 | 282 | ColumnVectorType m_tmp; |
ykuroda | 0:13a5d365ba16 | 283 | }; |
ykuroda | 0:13a5d365ba16 | 284 | |
ykuroda | 0:13a5d365ba16 | 285 | //template<typename MatrixType> |
ykuroda | 0:13a5d365ba16 | 286 | //typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const |
ykuroda | 0:13a5d365ba16 | 287 | //{ |
ykuroda | 0:13a5d365ba16 | 288 | // eigen_assert(m_isInitialized && "EigenSolver is not initialized."); |
ykuroda | 0:13a5d365ba16 | 289 | // eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); |
ykuroda | 0:13a5d365ba16 | 290 | // Index n = m_eivec.cols(); |
ykuroda | 0:13a5d365ba16 | 291 | // EigenvectorsType matV(n,n); |
ykuroda | 0:13a5d365ba16 | 292 | // // TODO |
ykuroda | 0:13a5d365ba16 | 293 | // return matV; |
ykuroda | 0:13a5d365ba16 | 294 | //} |
ykuroda | 0:13a5d365ba16 | 295 | |
ykuroda | 0:13a5d365ba16 | 296 | template<typename MatrixType> |
ykuroda | 0:13a5d365ba16 | 297 | GeneralizedEigenSolver<MatrixType>& |
ykuroda | 0:13a5d365ba16 | 298 | GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors) |
ykuroda | 0:13a5d365ba16 | 299 | { |
ykuroda | 0:13a5d365ba16 | 300 | check_template_parameters(); |
ykuroda | 0:13a5d365ba16 | 301 | |
ykuroda | 0:13a5d365ba16 | 302 | using std::sqrt; |
ykuroda | 0:13a5d365ba16 | 303 | using std::abs; |
ykuroda | 0:13a5d365ba16 | 304 | eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows()); |
ykuroda | 0:13a5d365ba16 | 305 | |
ykuroda | 0:13a5d365ba16 | 306 | // Reduce to generalized real Schur form: |
ykuroda | 0:13a5d365ba16 | 307 | // A = Q S Z and B = Q T Z |
ykuroda | 0:13a5d365ba16 | 308 | m_realQZ.compute(A, B, computeEigenvectors); |
ykuroda | 0:13a5d365ba16 | 309 | |
ykuroda | 0:13a5d365ba16 | 310 | if (m_realQZ.info() == Success) |
ykuroda | 0:13a5d365ba16 | 311 | { |
ykuroda | 0:13a5d365ba16 | 312 | m_matS = m_realQZ.matrixS(); |
ykuroda | 0:13a5d365ba16 | 313 | if (computeEigenvectors) |
ykuroda | 0:13a5d365ba16 | 314 | m_eivec = m_realQZ.matrixZ().transpose(); |
ykuroda | 0:13a5d365ba16 | 315 | |
ykuroda | 0:13a5d365ba16 | 316 | // Compute eigenvalues from matS |
ykuroda | 0:13a5d365ba16 | 317 | m_alphas.resize(A.cols()); |
ykuroda | 0:13a5d365ba16 | 318 | m_betas.resize(A.cols()); |
ykuroda | 0:13a5d365ba16 | 319 | Index i = 0; |
ykuroda | 0:13a5d365ba16 | 320 | while (i < A.cols()) |
ykuroda | 0:13a5d365ba16 | 321 | { |
ykuroda | 0:13a5d365ba16 | 322 | if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0)) |
ykuroda | 0:13a5d365ba16 | 323 | { |
ykuroda | 0:13a5d365ba16 | 324 | m_alphas.coeffRef(i) = m_matS.coeff(i, i); |
ykuroda | 0:13a5d365ba16 | 325 | m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i); |
ykuroda | 0:13a5d365ba16 | 326 | ++i; |
ykuroda | 0:13a5d365ba16 | 327 | } |
ykuroda | 0:13a5d365ba16 | 328 | else |
ykuroda | 0:13a5d365ba16 | 329 | { |
ykuroda | 0:13a5d365ba16 | 330 | Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1)); |
ykuroda | 0:13a5d365ba16 | 331 | Scalar z = sqrt(abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1))); |
ykuroda | 0:13a5d365ba16 | 332 | m_alphas.coeffRef(i) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z); |
ykuroda | 0:13a5d365ba16 | 333 | m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z); |
ykuroda | 0:13a5d365ba16 | 334 | |
ykuroda | 0:13a5d365ba16 | 335 | m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i); |
ykuroda | 0:13a5d365ba16 | 336 | m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i); |
ykuroda | 0:13a5d365ba16 | 337 | i += 2; |
ykuroda | 0:13a5d365ba16 | 338 | } |
ykuroda | 0:13a5d365ba16 | 339 | } |
ykuroda | 0:13a5d365ba16 | 340 | } |
ykuroda | 0:13a5d365ba16 | 341 | |
ykuroda | 0:13a5d365ba16 | 342 | m_isInitialized = true; |
ykuroda | 0:13a5d365ba16 | 343 | m_eigenvectorsOk = false;//computeEigenvectors; |
ykuroda | 0:13a5d365ba16 | 344 | |
ykuroda | 0:13a5d365ba16 | 345 | return *this; |
ykuroda | 0:13a5d365ba16 | 346 | } |
ykuroda | 0:13a5d365ba16 | 347 | |
ykuroda | 0:13a5d365ba16 | 348 | } // end namespace Eigen |
ykuroda | 0:13a5d365ba16 | 349 | |
ykuroda | 0:13a5d365ba16 | 350 | #endif // EIGEN_GENERALIZEDEIGENSOLVER_H |