Michael Ernst Peter / Eigen

src/Eigenvalues/ComplexEigenSolver.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) 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