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/ComplexEigenSolver.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) 2009 Claire Maurice |
ykuroda | 0:13a5d365ba16 | 5 | // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> |
ykuroda | 0:13a5d365ba16 | 6 | // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> |
ykuroda | 0:13a5d365ba16 | 7 | // |
ykuroda | 0:13a5d365ba16 | 8 | // This Source Code Form is subject to the terms of the Mozilla |
ykuroda | 0:13a5d365ba16 | 9 | // Public License v. 2.0. If a copy of the MPL was not distributed |
ykuroda | 0:13a5d365ba16 | 10 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
ykuroda | 0:13a5d365ba16 | 11 | |
ykuroda | 0:13a5d365ba16 | 12 | #ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H |
ykuroda | 0:13a5d365ba16 | 13 | #define EIGEN_COMPLEX_EIGEN_SOLVER_H |
ykuroda | 0:13a5d365ba16 | 14 | |
ykuroda | 0:13a5d365ba16 | 15 | #include "./ComplexSchur.h" |
ykuroda | 0:13a5d365ba16 | 16 | |
ykuroda | 0:13a5d365ba16 | 17 | namespace Eigen { |
ykuroda | 0:13a5d365ba16 | 18 | |
ykuroda | 0:13a5d365ba16 | 19 | /** \eigenvalues_module \ingroup Eigenvalues_Module |
ykuroda | 0:13a5d365ba16 | 20 | * |
ykuroda | 0:13a5d365ba16 | 21 | * |
ykuroda | 0:13a5d365ba16 | 22 | * \class ComplexEigenSolver |
ykuroda | 0:13a5d365ba16 | 23 | * |
ykuroda | 0:13a5d365ba16 | 24 | * \brief Computes eigenvalues and eigenvectors of general complex matrices |
ykuroda | 0:13a5d365ba16 | 25 | * |
ykuroda | 0:13a5d365ba16 | 26 | * \tparam _MatrixType the type of the matrix of which we are |
ykuroda | 0:13a5d365ba16 | 27 | * computing the eigendecomposition; this is expected to be an |
ykuroda | 0:13a5d365ba16 | 28 | * instantiation of the Matrix class template. |
ykuroda | 0:13a5d365ba16 | 29 | * |
ykuroda | 0:13a5d365ba16 | 30 | * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars |
ykuroda | 0:13a5d365ba16 | 31 | * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v |
ykuroda | 0:13a5d365ba16 | 32 | * \f$. If \f$ D \f$ is a diagonal matrix with the eigenvalues on |
ykuroda | 0:13a5d365ba16 | 33 | * the diagonal, and \f$ V \f$ is a matrix with the eigenvectors as |
ykuroda | 0:13a5d365ba16 | 34 | * its columns, then \f$ A V = V D \f$. The matrix \f$ V \f$ is |
ykuroda | 0:13a5d365ba16 | 35 | * almost always invertible, in which case we have \f$ A = V D V^{-1} |
ykuroda | 0:13a5d365ba16 | 36 | * \f$. This is called the eigendecomposition. |
ykuroda | 0:13a5d365ba16 | 37 | * |
ykuroda | 0:13a5d365ba16 | 38 | * The main function in this class is compute(), which computes the |
ykuroda | 0:13a5d365ba16 | 39 | * eigenvalues and eigenvectors of a given function. The |
ykuroda | 0:13a5d365ba16 | 40 | * documentation for that function contains an example showing the |
ykuroda | 0:13a5d365ba16 | 41 | * main features of the class. |
ykuroda | 0:13a5d365ba16 | 42 | * |
ykuroda | 0:13a5d365ba16 | 43 | * \sa class EigenSolver, class SelfAdjointEigenSolver |
ykuroda | 0:13a5d365ba16 | 44 | */ |
ykuroda | 0:13a5d365ba16 | 45 | template<typename _MatrixType> class ComplexEigenSolver |
ykuroda | 0:13a5d365ba16 | 46 | { |
ykuroda | 0:13a5d365ba16 | 47 | public: |
ykuroda | 0:13a5d365ba16 | 48 | |
ykuroda | 0:13a5d365ba16 | 49 | /** \brief Synonym for the template parameter \p _MatrixType. */ |
ykuroda | 0:13a5d365ba16 | 50 | typedef _MatrixType MatrixType; |
ykuroda | 0:13a5d365ba16 | 51 | |
ykuroda | 0:13a5d365ba16 | 52 | enum { |
ykuroda | 0:13a5d365ba16 | 53 | RowsAtCompileTime = MatrixType::RowsAtCompileTime, |
ykuroda | 0:13a5d365ba16 | 54 | ColsAtCompileTime = MatrixType::ColsAtCompileTime, |
ykuroda | 0:13a5d365ba16 | 55 | Options = MatrixType::Options, |
ykuroda | 0:13a5d365ba16 | 56 | MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, |
ykuroda | 0:13a5d365ba16 | 57 | MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime |
ykuroda | 0:13a5d365ba16 | 58 | }; |
ykuroda | 0:13a5d365ba16 | 59 | |
ykuroda | 0:13a5d365ba16 | 60 | /** \brief Scalar type for matrices of type #MatrixType. */ |
ykuroda | 0:13a5d365ba16 | 61 | typedef typename MatrixType::Scalar Scalar; |
ykuroda | 0:13a5d365ba16 | 62 | typedef typename NumTraits<Scalar>::Real RealScalar; |
ykuroda | 0:13a5d365ba16 | 63 | typedef typename MatrixType::Index Index; |
ykuroda | 0:13a5d365ba16 | 64 | |
ykuroda | 0:13a5d365ba16 | 65 | /** \brief Complex scalar type for #MatrixType. |
ykuroda | 0:13a5d365ba16 | 66 | * |
ykuroda | 0:13a5d365ba16 | 67 | * This is \c std::complex<Scalar> if #Scalar is real (e.g., |
ykuroda | 0:13a5d365ba16 | 68 | * \c float or \c double) and just \c Scalar if #Scalar is |
ykuroda | 0:13a5d365ba16 | 69 | * complex. |
ykuroda | 0:13a5d365ba16 | 70 | */ |
ykuroda | 0:13a5d365ba16 | 71 | typedef std::complex<RealScalar> ComplexScalar; |
ykuroda | 0:13a5d365ba16 | 72 | |
ykuroda | 0:13a5d365ba16 | 73 | /** \brief Type for vector of eigenvalues as returned by eigenvalues(). |
ykuroda | 0:13a5d365ba16 | 74 | * |
ykuroda | 0:13a5d365ba16 | 75 | * This is a column vector with entries of type #ComplexScalar. |
ykuroda | 0:13a5d365ba16 | 76 | * The length of the vector is the size of #MatrixType. |
ykuroda | 0:13a5d365ba16 | 77 | */ |
ykuroda | 0:13a5d365ba16 | 78 | typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1> EigenvalueType; |
ykuroda | 0:13a5d365ba16 | 79 | |
ykuroda | 0:13a5d365ba16 | 80 | /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). |
ykuroda | 0:13a5d365ba16 | 81 | * |
ykuroda | 0:13a5d365ba16 | 82 | * This is a square matrix with entries of type #ComplexScalar. |
ykuroda | 0:13a5d365ba16 | 83 | * The size is the same as the size of #MatrixType. |
ykuroda | 0:13a5d365ba16 | 84 | */ |
ykuroda | 0:13a5d365ba16 | 85 | typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorType; |
ykuroda | 0:13a5d365ba16 | 86 | |
ykuroda | 0:13a5d365ba16 | 87 | /** \brief Default constructor. |
ykuroda | 0:13a5d365ba16 | 88 | * |
ykuroda | 0:13a5d365ba16 | 89 | * The default constructor is useful in cases in which the user intends to |
ykuroda | 0:13a5d365ba16 | 90 | * perform decompositions via compute(). |
ykuroda | 0:13a5d365ba16 | 91 | */ |
ykuroda | 0:13a5d365ba16 | 92 | ComplexEigenSolver() |
ykuroda | 0:13a5d365ba16 | 93 | : m_eivec(), |
ykuroda | 0:13a5d365ba16 | 94 | m_eivalues(), |
ykuroda | 0:13a5d365ba16 | 95 | m_schur(), |
ykuroda | 0:13a5d365ba16 | 96 | m_isInitialized(false), |
ykuroda | 0:13a5d365ba16 | 97 | m_eigenvectorsOk(false), |
ykuroda | 0:13a5d365ba16 | 98 | m_matX() |
ykuroda | 0:13a5d365ba16 | 99 | {} |
ykuroda | 0:13a5d365ba16 | 100 | |
ykuroda | 0:13a5d365ba16 | 101 | /** \brief Default Constructor with memory preallocation |
ykuroda | 0:13a5d365ba16 | 102 | * |
ykuroda | 0:13a5d365ba16 | 103 | * Like the default constructor but with preallocation of the internal data |
ykuroda | 0:13a5d365ba16 | 104 | * according to the specified problem \a size. |
ykuroda | 0:13a5d365ba16 | 105 | * \sa ComplexEigenSolver() |
ykuroda | 0:13a5d365ba16 | 106 | */ |
ykuroda | 0:13a5d365ba16 | 107 | ComplexEigenSolver(Index size) |
ykuroda | 0:13a5d365ba16 | 108 | : m_eivec(size, size), |
ykuroda | 0:13a5d365ba16 | 109 | m_eivalues(size), |
ykuroda | 0:13a5d365ba16 | 110 | m_schur(size), |
ykuroda | 0:13a5d365ba16 | 111 | m_isInitialized(false), |
ykuroda | 0:13a5d365ba16 | 112 | m_eigenvectorsOk(false), |
ykuroda | 0:13a5d365ba16 | 113 | m_matX(size, size) |
ykuroda | 0:13a5d365ba16 | 114 | {} |
ykuroda | 0:13a5d365ba16 | 115 | |
ykuroda | 0:13a5d365ba16 | 116 | /** \brief Constructor; computes eigendecomposition of given matrix. |
ykuroda | 0:13a5d365ba16 | 117 | * |
ykuroda | 0:13a5d365ba16 | 118 | * \param[in] matrix Square matrix whose eigendecomposition is to be computed. |
ykuroda | 0:13a5d365ba16 | 119 | * \param[in] computeEigenvectors If true, both the eigenvectors and the |
ykuroda | 0:13a5d365ba16 | 120 | * eigenvalues are computed; if false, only the eigenvalues are |
ykuroda | 0:13a5d365ba16 | 121 | * computed. |
ykuroda | 0:13a5d365ba16 | 122 | * |
ykuroda | 0:13a5d365ba16 | 123 | * This constructor calls compute() to compute the eigendecomposition. |
ykuroda | 0:13a5d365ba16 | 124 | */ |
ykuroda | 0:13a5d365ba16 | 125 | ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) |
ykuroda | 0:13a5d365ba16 | 126 | : m_eivec(matrix.rows(),matrix.cols()), |
ykuroda | 0:13a5d365ba16 | 127 | m_eivalues(matrix.cols()), |
ykuroda | 0:13a5d365ba16 | 128 | m_schur(matrix.rows()), |
ykuroda | 0:13a5d365ba16 | 129 | m_isInitialized(false), |
ykuroda | 0:13a5d365ba16 | 130 | m_eigenvectorsOk(false), |
ykuroda | 0:13a5d365ba16 | 131 | m_matX(matrix.rows(),matrix.cols()) |
ykuroda | 0:13a5d365ba16 | 132 | { |
ykuroda | 0:13a5d365ba16 | 133 | compute(matrix, computeEigenvectors); |
ykuroda | 0:13a5d365ba16 | 134 | } |
ykuroda | 0:13a5d365ba16 | 135 | |
ykuroda | 0:13a5d365ba16 | 136 | /** \brief Returns the eigenvectors of given matrix. |
ykuroda | 0:13a5d365ba16 | 137 | * |
ykuroda | 0:13a5d365ba16 | 138 | * \returns A const reference to the matrix whose columns are the eigenvectors. |
ykuroda | 0:13a5d365ba16 | 139 | * |
ykuroda | 0:13a5d365ba16 | 140 | * \pre Either the constructor |
ykuroda | 0:13a5d365ba16 | 141 | * ComplexEigenSolver(const MatrixType& matrix, bool) or the member |
ykuroda | 0:13a5d365ba16 | 142 | * function compute(const MatrixType& matrix, bool) has been called before |
ykuroda | 0:13a5d365ba16 | 143 | * to compute the eigendecomposition of a matrix, and |
ykuroda | 0:13a5d365ba16 | 144 | * \p computeEigenvectors was set to true (the default). |
ykuroda | 0:13a5d365ba16 | 145 | * |
ykuroda | 0:13a5d365ba16 | 146 | * This function returns a matrix whose columns are the eigenvectors. Column |
ykuroda | 0:13a5d365ba16 | 147 | * \f$ k \f$ is an eigenvector corresponding to eigenvalue number \f$ k |
ykuroda | 0:13a5d365ba16 | 148 | * \f$ as returned by eigenvalues(). The eigenvectors are normalized to |
ykuroda | 0:13a5d365ba16 | 149 | * have (Euclidean) norm equal to one. The matrix returned by this |
ykuroda | 0:13a5d365ba16 | 150 | * function is the matrix \f$ V \f$ in the eigendecomposition \f$ A = V D |
ykuroda | 0:13a5d365ba16 | 151 | * V^{-1} \f$, if it exists. |
ykuroda | 0:13a5d365ba16 | 152 | * |
ykuroda | 0:13a5d365ba16 | 153 | * Example: \include ComplexEigenSolver_eigenvectors.cpp |
ykuroda | 0:13a5d365ba16 | 154 | * Output: \verbinclude ComplexEigenSolver_eigenvectors.out |
ykuroda | 0:13a5d365ba16 | 155 | */ |
ykuroda | 0:13a5d365ba16 | 156 | const EigenvectorType& eigenvectors() const |
ykuroda | 0:13a5d365ba16 | 157 | { |
ykuroda | 0:13a5d365ba16 | 158 | eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); |
ykuroda | 0:13a5d365ba16 | 159 | eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); |
ykuroda | 0:13a5d365ba16 | 160 | return m_eivec; |
ykuroda | 0:13a5d365ba16 | 161 | } |
ykuroda | 0:13a5d365ba16 | 162 | |
ykuroda | 0:13a5d365ba16 | 163 | /** \brief Returns the eigenvalues of given matrix. |
ykuroda | 0:13a5d365ba16 | 164 | * |
ykuroda | 0:13a5d365ba16 | 165 | * \returns A const reference to the column vector containing the eigenvalues. |
ykuroda | 0:13a5d365ba16 | 166 | * |
ykuroda | 0:13a5d365ba16 | 167 | * \pre Either the constructor |
ykuroda | 0:13a5d365ba16 | 168 | * ComplexEigenSolver(const MatrixType& matrix, bool) or the member |
ykuroda | 0:13a5d365ba16 | 169 | * function compute(const MatrixType& matrix, bool) has been called before |
ykuroda | 0:13a5d365ba16 | 170 | * to compute the eigendecomposition of a matrix. |
ykuroda | 0:13a5d365ba16 | 171 | * |
ykuroda | 0:13a5d365ba16 | 172 | * This function returns a column vector containing the |
ykuroda | 0:13a5d365ba16 | 173 | * eigenvalues. Eigenvalues are repeated according to their |
ykuroda | 0:13a5d365ba16 | 174 | * algebraic multiplicity, so there are as many eigenvalues as |
ykuroda | 0:13a5d365ba16 | 175 | * rows in the matrix. The eigenvalues are not sorted in any particular |
ykuroda | 0:13a5d365ba16 | 176 | * order. |
ykuroda | 0:13a5d365ba16 | 177 | * |
ykuroda | 0:13a5d365ba16 | 178 | * Example: \include ComplexEigenSolver_eigenvalues.cpp |
ykuroda | 0:13a5d365ba16 | 179 | * Output: \verbinclude ComplexEigenSolver_eigenvalues.out |
ykuroda | 0:13a5d365ba16 | 180 | */ |
ykuroda | 0:13a5d365ba16 | 181 | const EigenvalueType& eigenvalues() const |
ykuroda | 0:13a5d365ba16 | 182 | { |
ykuroda | 0:13a5d365ba16 | 183 | eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); |
ykuroda | 0:13a5d365ba16 | 184 | return m_eivalues; |
ykuroda | 0:13a5d365ba16 | 185 | } |
ykuroda | 0:13a5d365ba16 | 186 | |
ykuroda | 0:13a5d365ba16 | 187 | /** \brief Computes eigendecomposition of given matrix. |
ykuroda | 0:13a5d365ba16 | 188 | * |
ykuroda | 0:13a5d365ba16 | 189 | * \param[in] matrix Square matrix whose eigendecomposition is to be computed. |
ykuroda | 0:13a5d365ba16 | 190 | * \param[in] computeEigenvectors If true, both the eigenvectors and the |
ykuroda | 0:13a5d365ba16 | 191 | * eigenvalues are computed; if false, only the eigenvalues are |
ykuroda | 0:13a5d365ba16 | 192 | * computed. |
ykuroda | 0:13a5d365ba16 | 193 | * \returns Reference to \c *this |
ykuroda | 0:13a5d365ba16 | 194 | * |
ykuroda | 0:13a5d365ba16 | 195 | * This function computes the eigenvalues of the complex matrix \p matrix. |
ykuroda | 0:13a5d365ba16 | 196 | * The eigenvalues() function can be used to retrieve them. If |
ykuroda | 0:13a5d365ba16 | 197 | * \p computeEigenvectors is true, then the eigenvectors are also computed |
ykuroda | 0:13a5d365ba16 | 198 | * and can be retrieved by calling eigenvectors(). |
ykuroda | 0:13a5d365ba16 | 199 | * |
ykuroda | 0:13a5d365ba16 | 200 | * The matrix is first reduced to Schur form using the |
ykuroda | 0:13a5d365ba16 | 201 | * ComplexSchur class. The Schur decomposition is then used to |
ykuroda | 0:13a5d365ba16 | 202 | * compute the eigenvalues and eigenvectors. |
ykuroda | 0:13a5d365ba16 | 203 | * |
ykuroda | 0:13a5d365ba16 | 204 | * The cost of the computation is dominated by the cost of the |
ykuroda | 0:13a5d365ba16 | 205 | * Schur decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$ |
ykuroda | 0:13a5d365ba16 | 206 | * is the size of the matrix. |
ykuroda | 0:13a5d365ba16 | 207 | * |
ykuroda | 0:13a5d365ba16 | 208 | * Example: \include ComplexEigenSolver_compute.cpp |
ykuroda | 0:13a5d365ba16 | 209 | * Output: \verbinclude ComplexEigenSolver_compute.out |
ykuroda | 0:13a5d365ba16 | 210 | */ |
ykuroda | 0:13a5d365ba16 | 211 | ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); |
ykuroda | 0:13a5d365ba16 | 212 | |
ykuroda | 0:13a5d365ba16 | 213 | /** \brief Reports whether previous computation was successful. |
ykuroda | 0:13a5d365ba16 | 214 | * |
ykuroda | 0:13a5d365ba16 | 215 | * \returns \c Success if computation was succesful, \c NoConvergence otherwise. |
ykuroda | 0:13a5d365ba16 | 216 | */ |
ykuroda | 0:13a5d365ba16 | 217 | ComputationInfo info() const |
ykuroda | 0:13a5d365ba16 | 218 | { |
ykuroda | 0:13a5d365ba16 | 219 | eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); |
ykuroda | 0:13a5d365ba16 | 220 | return m_schur.info(); |
ykuroda | 0:13a5d365ba16 | 221 | } |
ykuroda | 0:13a5d365ba16 | 222 | |
ykuroda | 0:13a5d365ba16 | 223 | /** \brief Sets the maximum number of iterations allowed. */ |
ykuroda | 0:13a5d365ba16 | 224 | ComplexEigenSolver& setMaxIterations(Index maxIters) |
ykuroda | 0:13a5d365ba16 | 225 | { |
ykuroda | 0:13a5d365ba16 | 226 | m_schur.setMaxIterations(maxIters); |
ykuroda | 0:13a5d365ba16 | 227 | return *this; |
ykuroda | 0:13a5d365ba16 | 228 | } |
ykuroda | 0:13a5d365ba16 | 229 | |
ykuroda | 0:13a5d365ba16 | 230 | /** \brief Returns the maximum number of iterations. */ |
ykuroda | 0:13a5d365ba16 | 231 | Index getMaxIterations() |
ykuroda | 0:13a5d365ba16 | 232 | { |
ykuroda | 0:13a5d365ba16 | 233 | return m_schur.getMaxIterations(); |
ykuroda | 0:13a5d365ba16 | 234 | } |
ykuroda | 0:13a5d365ba16 | 235 | |
ykuroda | 0:13a5d365ba16 | 236 | protected: |
ykuroda | 0:13a5d365ba16 | 237 | |
ykuroda | 0:13a5d365ba16 | 238 | static void check_template_parameters() |
ykuroda | 0:13a5d365ba16 | 239 | { |
ykuroda | 0:13a5d365ba16 | 240 | EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); |
ykuroda | 0:13a5d365ba16 | 241 | } |
ykuroda | 0:13a5d365ba16 | 242 | |
ykuroda | 0:13a5d365ba16 | 243 | EigenvectorType m_eivec; |
ykuroda | 0:13a5d365ba16 | 244 | EigenvalueType m_eivalues; |
ykuroda | 0:13a5d365ba16 | 245 | ComplexSchur<MatrixType> m_schur; |
ykuroda | 0:13a5d365ba16 | 246 | bool m_isInitialized; |
ykuroda | 0:13a5d365ba16 | 247 | bool m_eigenvectorsOk; |
ykuroda | 0:13a5d365ba16 | 248 | EigenvectorType m_matX; |
ykuroda | 0:13a5d365ba16 | 249 | |
ykuroda | 0:13a5d365ba16 | 250 | private: |
ykuroda | 0:13a5d365ba16 | 251 | void doComputeEigenvectors(const RealScalar& matrixnorm); |
ykuroda | 0:13a5d365ba16 | 252 | void sortEigenvalues(bool computeEigenvectors); |
ykuroda | 0:13a5d365ba16 | 253 | }; |
ykuroda | 0:13a5d365ba16 | 254 | |
ykuroda | 0:13a5d365ba16 | 255 | |
ykuroda | 0:13a5d365ba16 | 256 | template<typename MatrixType> |
ykuroda | 0:13a5d365ba16 | 257 | ComplexEigenSolver<MatrixType>& |
ykuroda | 0:13a5d365ba16 | 258 | ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) |
ykuroda | 0:13a5d365ba16 | 259 | { |
ykuroda | 0:13a5d365ba16 | 260 | check_template_parameters(); |
ykuroda | 0:13a5d365ba16 | 261 | |
ykuroda | 0:13a5d365ba16 | 262 | // this code is inspired from Jampack |
ykuroda | 0:13a5d365ba16 | 263 | eigen_assert(matrix.cols() == matrix.rows()); |
ykuroda | 0:13a5d365ba16 | 264 | |
ykuroda | 0:13a5d365ba16 | 265 | // Do a complex Schur decomposition, A = U T U^* |
ykuroda | 0:13a5d365ba16 | 266 | // The eigenvalues are on the diagonal of T. |
ykuroda | 0:13a5d365ba16 | 267 | m_schur.compute(matrix, computeEigenvectors); |
ykuroda | 0:13a5d365ba16 | 268 | |
ykuroda | 0:13a5d365ba16 | 269 | if(m_schur.info() == Success) |
ykuroda | 0:13a5d365ba16 | 270 | { |
ykuroda | 0:13a5d365ba16 | 271 | m_eivalues = m_schur.matrixT().diagonal(); |
ykuroda | 0:13a5d365ba16 | 272 | if(computeEigenvectors) |
ykuroda | 0:13a5d365ba16 | 273 | doComputeEigenvectors(matrix.norm()); |
ykuroda | 0:13a5d365ba16 | 274 | sortEigenvalues(computeEigenvectors); |
ykuroda | 0:13a5d365ba16 | 275 | } |
ykuroda | 0:13a5d365ba16 | 276 | |
ykuroda | 0:13a5d365ba16 | 277 | m_isInitialized = true; |
ykuroda | 0:13a5d365ba16 | 278 | m_eigenvectorsOk = computeEigenvectors; |
ykuroda | 0:13a5d365ba16 | 279 | return *this; |
ykuroda | 0:13a5d365ba16 | 280 | } |
ykuroda | 0:13a5d365ba16 | 281 | |
ykuroda | 0:13a5d365ba16 | 282 | |
ykuroda | 0:13a5d365ba16 | 283 | template<typename MatrixType> |
ykuroda | 0:13a5d365ba16 | 284 | void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(const RealScalar& matrixnorm) |
ykuroda | 0:13a5d365ba16 | 285 | { |
ykuroda | 0:13a5d365ba16 | 286 | const Index n = m_eivalues.size(); |
ykuroda | 0:13a5d365ba16 | 287 | |
ykuroda | 0:13a5d365ba16 | 288 | // Compute X such that T = X D X^(-1), where D is the diagonal of T. |
ykuroda | 0:13a5d365ba16 | 289 | // The matrix X is unit triangular. |
ykuroda | 0:13a5d365ba16 | 290 | m_matX = EigenvectorType::Zero(n, n); |
ykuroda | 0:13a5d365ba16 | 291 | for(Index k=n-1 ; k>=0 ; k--) |
ykuroda | 0:13a5d365ba16 | 292 | { |
ykuroda | 0:13a5d365ba16 | 293 | m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0); |
ykuroda | 0:13a5d365ba16 | 294 | // Compute X(i,k) using the (i,k) entry of the equation X T = D X |
ykuroda | 0:13a5d365ba16 | 295 | for(Index i=k-1 ; i>=0 ; i--) |
ykuroda | 0:13a5d365ba16 | 296 | { |
ykuroda | 0:13a5d365ba16 | 297 | m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k); |
ykuroda | 0:13a5d365ba16 | 298 | if(k-i-1>0) |
ykuroda | 0:13a5d365ba16 | 299 | m_matX.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * m_matX.col(k).segment(i+1,k-i-1)).value(); |
ykuroda | 0:13a5d365ba16 | 300 | ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k); |
ykuroda | 0:13a5d365ba16 | 301 | if(z==ComplexScalar(0)) |
ykuroda | 0:13a5d365ba16 | 302 | { |
ykuroda | 0:13a5d365ba16 | 303 | // If the i-th and k-th eigenvalue are equal, then z equals 0. |
ykuroda | 0:13a5d365ba16 | 304 | // Use a small value instead, to prevent division by zero. |
ykuroda | 0:13a5d365ba16 | 305 | numext::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm; |
ykuroda | 0:13a5d365ba16 | 306 | } |
ykuroda | 0:13a5d365ba16 | 307 | m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z; |
ykuroda | 0:13a5d365ba16 | 308 | } |
ykuroda | 0:13a5d365ba16 | 309 | } |
ykuroda | 0:13a5d365ba16 | 310 | |
ykuroda | 0:13a5d365ba16 | 311 | // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1) |
ykuroda | 0:13a5d365ba16 | 312 | m_eivec.noalias() = m_schur.matrixU() * m_matX; |
ykuroda | 0:13a5d365ba16 | 313 | // .. and normalize the eigenvectors |
ykuroda | 0:13a5d365ba16 | 314 | for(Index k=0 ; k<n ; k++) |
ykuroda | 0:13a5d365ba16 | 315 | { |
ykuroda | 0:13a5d365ba16 | 316 | m_eivec.col(k).normalize(); |
ykuroda | 0:13a5d365ba16 | 317 | } |
ykuroda | 0:13a5d365ba16 | 318 | } |
ykuroda | 0:13a5d365ba16 | 319 | |
ykuroda | 0:13a5d365ba16 | 320 | |
ykuroda | 0:13a5d365ba16 | 321 | template<typename MatrixType> |
ykuroda | 0:13a5d365ba16 | 322 | void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors) |
ykuroda | 0:13a5d365ba16 | 323 | { |
ykuroda | 0:13a5d365ba16 | 324 | const Index n = m_eivalues.size(); |
ykuroda | 0:13a5d365ba16 | 325 | for (Index i=0; i<n; i++) |
ykuroda | 0:13a5d365ba16 | 326 | { |
ykuroda | 0:13a5d365ba16 | 327 | Index k; |
ykuroda | 0:13a5d365ba16 | 328 | m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k); |
ykuroda | 0:13a5d365ba16 | 329 | if (k != 0) |
ykuroda | 0:13a5d365ba16 | 330 | { |
ykuroda | 0:13a5d365ba16 | 331 | k += i; |
ykuroda | 0:13a5d365ba16 | 332 | std::swap(m_eivalues[k],m_eivalues[i]); |
ykuroda | 0:13a5d365ba16 | 333 | if(computeEigenvectors) |
ykuroda | 0:13a5d365ba16 | 334 | m_eivec.col(i).swap(m_eivec.col(k)); |
ykuroda | 0:13a5d365ba16 | 335 | } |
ykuroda | 0:13a5d365ba16 | 336 | } |
ykuroda | 0:13a5d365ba16 | 337 | } |
ykuroda | 0:13a5d365ba16 | 338 | |
ykuroda | 0:13a5d365ba16 | 339 | } // end namespace Eigen |
ykuroda | 0:13a5d365ba16 | 340 | |
ykuroda | 0:13a5d365ba16 | 341 | #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H |