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 Eigen::Index Index; ///< \deprecated since Eigen 3.3 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 explicit 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 template<typename InputType> 126 explicit ComplexEigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true) 127 : m_eivec(matrix.rows(),matrix.cols()), 128 m_eivalues(matrix.cols()), 129 m_schur(matrix.rows()), 130 m_isInitialized(false), 131 m_eigenvectorsOk(false), 132 m_matX(matrix.rows(),matrix.cols()) 133 { 134 compute(matrix.derived(), computeEigenvectors); 135 } 136 137 /** \brief Returns the eigenvectors of given matrix. 138 * 139 * \returns A const reference to the matrix whose columns are the eigenvectors. 140 * 141 * \pre Either the constructor 142 * ComplexEigenSolver(const MatrixType& matrix, bool) or the member 143 * function compute(const MatrixType& matrix, bool) has been called before 144 * to compute the eigendecomposition of a matrix, and 145 * \p computeEigenvectors was set to true (the default). 146 * 147 * This function returns a matrix whose columns are the eigenvectors. Column 148 * \f$ k \f$ is an eigenvector corresponding to eigenvalue number \f$ k 149 * \f$ as returned by eigenvalues(). The eigenvectors are normalized to 150 * have (Euclidean) norm equal to one. The matrix returned by this 151 * function is the matrix \f$ V \f$ in the eigendecomposition \f$ A = V D 152 * V^{-1} \f$, if it exists. 153 * 154 * Example: \include ComplexEigenSolver_eigenvectors.cpp 155 * Output: \verbinclude ComplexEigenSolver_eigenvectors.out 156 */ 157 const EigenvectorType& eigenvectors() const 158 { 159 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 160 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 161 return m_eivec; 162 } 163 164 /** \brief Returns the eigenvalues of given matrix. 165 * 166 * \returns A const reference to the column vector containing the eigenvalues. 167 * 168 * \pre Either the constructor 169 * ComplexEigenSolver(const MatrixType& matrix, bool) or the member 170 * function compute(const MatrixType& matrix, bool) has been called before 171 * to compute the eigendecomposition of a matrix. 172 * 173 * This function returns a column vector containing the 174 * eigenvalues. Eigenvalues are repeated according to their 175 * algebraic multiplicity, so there are as many eigenvalues as 176 * rows in the matrix. The eigenvalues are not sorted in any particular 177 * order. 178 * 179 * Example: \include ComplexEigenSolver_eigenvalues.cpp 180 * Output: \verbinclude ComplexEigenSolver_eigenvalues.out 181 */ 182 const EigenvalueType& eigenvalues() const 183 { 184 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 185 return m_eivalues; 186 } 187 188 /** \brief Computes eigendecomposition of given matrix. 189 * 190 * \param[in] matrix Square matrix whose eigendecomposition is to be computed. 191 * \param[in] computeEigenvectors If true, both the eigenvectors and the 192 * eigenvalues are computed; if false, only the eigenvalues are 193 * computed. 194 * \returns Reference to \c *this 195 * 196 * This function computes the eigenvalues of the complex matrix \p matrix. 197 * The eigenvalues() function can be used to retrieve them. If 198 * \p computeEigenvectors is true, then the eigenvectors are also computed 199 * and can be retrieved by calling eigenvectors(). 200 * 201 * The matrix is first reduced to Schur form using the 202 * ComplexSchur class. The Schur decomposition is then used to 203 * compute the eigenvalues and eigenvectors. 204 * 205 * The cost of the computation is dominated by the cost of the 206 * Schur decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$ 207 * is the size of the matrix. 208 * 209 * Example: \include ComplexEigenSolver_compute.cpp 210 * Output: \verbinclude ComplexEigenSolver_compute.out 211 */ 212 template<typename InputType> 213 ComplexEigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true); 214 215 /** \brief Reports whether previous computation was successful. 216 * 217 * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 218 */ 219 ComputationInfo info() const 220 { 221 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 222 return m_schur.info(); 223 } 224 225 /** \brief Sets the maximum number of iterations allowed. */ 226 ComplexEigenSolver& setMaxIterations(Index maxIters) 227 { 228 m_schur.setMaxIterations(maxIters); 229 return *this; 230 } 231 232 /** \brief Returns the maximum number of iterations. */ 233 Index getMaxIterations() 234 { 235 return m_schur.getMaxIterations(); 236 } 237 238 protected: 239 240 static void check_template_parameters() 241 { 242 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 243 } 244 245 EigenvectorType m_eivec; 246 EigenvalueType m_eivalues; 247 ComplexSchur<MatrixType> m_schur; 248 bool m_isInitialized; 249 bool m_eigenvectorsOk; 250 EigenvectorType m_matX; 251 252 private: 253 void doComputeEigenvectors(RealScalar matrixnorm); 254 void sortEigenvalues(bool computeEigenvectors); 255 }; 256 257 258 template<typename MatrixType> 259 template<typename InputType> 260 ComplexEigenSolver<MatrixType>& 261 ComplexEigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors) 262 { 263 check_template_parameters(); 264 265 // this code is inspired from Jampack 266 eigen_assert(matrix.cols() == matrix.rows()); 267 268 // Do a complex Schur decomposition, A = U T U^* 269 // The eigenvalues are on the diagonal of T. 270 m_schur.compute(matrix.derived(), computeEigenvectors); 271 272 if(m_schur.info() == Success) 273 { 274 m_eivalues = m_schur.matrixT().diagonal(); 275 if(computeEigenvectors) 276 doComputeEigenvectors(m_schur.matrixT().norm()); 277 sortEigenvalues(computeEigenvectors); 278 } 279 280 m_isInitialized = true; 281 m_eigenvectorsOk = computeEigenvectors; 282 return *this; 283 } 284 285 286 template<typename MatrixType> 287 void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm) 288 { 289 const Index n = m_eivalues.size(); 290 291 matrixnorm = numext::maxi(matrixnorm,(std::numeric_limits<RealScalar>::min)()); 292 293 // Compute X such that T = X D X^(-1), where D is the diagonal of T. 294 // The matrix X is unit triangular. 295 m_matX = EigenvectorType::Zero(n, n); 296 for(Index k=n-1 ; k>=0 ; k--) 297 { 298 m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0); 299 // Compute X(i,k) using the (i,k) entry of the equation X T = D X 300 for(Index i=k-1 ; i>=0 ; i--) 301 { 302 m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k); 303 if(k-i-1>0) 304 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(); 305 ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k); 306 if(z==ComplexScalar(0)) 307 { 308 // If the i-th and k-th eigenvalue are equal, then z equals 0. 309 // Use a small value instead, to prevent division by zero. 310 numext::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm; 311 } 312 m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z; 313 } 314 } 315 316 // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1) 317 m_eivec.noalias() = m_schur.matrixU() * m_matX; 318 // .. and normalize the eigenvectors 319 for(Index k=0 ; k<n ; k++) 320 { 321 m_eivec.col(k).normalize(); 322 } 323 } 324 325 326 template<typename MatrixType> 327 void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors) 328 { 329 const Index n = m_eivalues.size(); 330 for (Index i=0; i<n; i++) 331 { 332 Index k; 333 m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k); 334 if (k != 0) 335 { 336 k += i; 337 std::swap(m_eivalues[k],m_eivalues[i]); 338 if(computeEigenvectors) 339 m_eivec.col(i).swap(m_eivec.col(k)); 340 } 341 } 342 } 343 344 } // end namespace Eigen 345 346 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H 347
