Michael Ernst Peter / Eigen

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?

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) 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