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