Eigne Matrix Class Library
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
ComplexEigenSolver.h
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2009 Claire Maurice 00005 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 00006 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> 00007 // 00008 // This Source Code Form is subject to the terms of the Mozilla 00009 // Public License v. 2.0. If a copy of the MPL was not distributed 00010 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00011 00012 #ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H 00013 #define EIGEN_COMPLEX_EIGEN_SOLVER_H 00014 00015 #include "./ComplexSchur.h" 00016 00017 namespace Eigen { 00018 00019 /** \eigenvalues_module \ingroup Eigenvalues_Module 00020 * 00021 * 00022 * \class ComplexEigenSolver 00023 * 00024 * \brief Computes eigenvalues and eigenvectors of general complex matrices 00025 * 00026 * \tparam _MatrixType the type of the matrix of which we are 00027 * computing the eigendecomposition; this is expected to be an 00028 * instantiation of the Matrix class template. 00029 * 00030 * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars 00031 * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v 00032 * \f$. If \f$ D \f$ is a diagonal matrix with the eigenvalues on 00033 * the diagonal, and \f$ V \f$ is a matrix with the eigenvectors as 00034 * its columns, then \f$ A V = V D \f$. The matrix \f$ V \f$ is 00035 * almost always invertible, in which case we have \f$ A = V D V^{-1} 00036 * \f$. This is called the eigendecomposition. 00037 * 00038 * The main function in this class is compute(), which computes the 00039 * eigenvalues and eigenvectors of a given function. The 00040 * documentation for that function contains an example showing the 00041 * main features of the class. 00042 * 00043 * \sa class EigenSolver, class SelfAdjointEigenSolver 00044 */ 00045 template<typename _MatrixType> class ComplexEigenSolver 00046 { 00047 public: 00048 00049 /** \brief Synonym for the template parameter \p _MatrixType. */ 00050 typedef _MatrixType MatrixType; 00051 00052 enum { 00053 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00054 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00055 Options = MatrixType::Options, 00056 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00057 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00058 }; 00059 00060 /** \brief Scalar type for matrices of type #MatrixType. */ 00061 typedef typename MatrixType::Scalar Scalar; 00062 typedef typename NumTraits<Scalar>::Real RealScalar; 00063 typedef typename MatrixType::Index Index; 00064 00065 /** \brief Complex scalar type for #MatrixType. 00066 * 00067 * This is \c std::complex<Scalar> if #Scalar is real (e.g., 00068 * \c float or \c double) and just \c Scalar if #Scalar is 00069 * complex. 00070 */ 00071 typedef std::complex<RealScalar> ComplexScalar; 00072 00073 /** \brief Type for vector of eigenvalues as returned by eigenvalues(). 00074 * 00075 * This is a column vector with entries of type #ComplexScalar. 00076 * The length of the vector is the size of #MatrixType. 00077 */ 00078 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1> EigenvalueType; 00079 00080 /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). 00081 * 00082 * This is a square matrix with entries of type #ComplexScalar. 00083 * The size is the same as the size of #MatrixType. 00084 */ 00085 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorType; 00086 00087 /** \brief Default constructor. 00088 * 00089 * The default constructor is useful in cases in which the user intends to 00090 * perform decompositions via compute(). 00091 */ 00092 ComplexEigenSolver() 00093 : m_eivec(), 00094 m_eivalues(), 00095 m_schur(), 00096 m_isInitialized(false), 00097 m_eigenvectorsOk(false), 00098 m_matX() 00099 {} 00100 00101 /** \brief Default Constructor with memory preallocation 00102 * 00103 * Like the default constructor but with preallocation of the internal data 00104 * according to the specified problem \a size. 00105 * \sa ComplexEigenSolver() 00106 */ 00107 ComplexEigenSolver(Index size) 00108 : m_eivec(size, size), 00109 m_eivalues(size), 00110 m_schur(size), 00111 m_isInitialized(false), 00112 m_eigenvectorsOk(false), 00113 m_matX(size, size) 00114 {} 00115 00116 /** \brief Constructor; computes eigendecomposition of given matrix. 00117 * 00118 * \param[in] matrix Square matrix whose eigendecomposition is to be computed. 00119 * \param[in] computeEigenvectors If true, both the eigenvectors and the 00120 * eigenvalues are computed; if false, only the eigenvalues are 00121 * computed. 00122 * 00123 * This constructor calls compute() to compute the eigendecomposition. 00124 */ 00125 ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) 00126 : m_eivec(matrix.rows(),matrix.cols()), 00127 m_eivalues(matrix.cols()), 00128 m_schur(matrix.rows()), 00129 m_isInitialized(false), 00130 m_eigenvectorsOk(false), 00131 m_matX(matrix.rows(),matrix.cols()) 00132 { 00133 compute(matrix, computeEigenvectors); 00134 } 00135 00136 /** \brief Returns the eigenvectors of given matrix. 00137 * 00138 * \returns A const reference to the matrix whose columns are the eigenvectors. 00139 * 00140 * \pre Either the constructor 00141 * ComplexEigenSolver(const MatrixType& matrix, bool) or the member 00142 * function compute(const MatrixType& matrix, bool) has been called before 00143 * to compute the eigendecomposition of a matrix, and 00144 * \p computeEigenvectors was set to true (the default). 00145 * 00146 * This function returns a matrix whose columns are the eigenvectors. Column 00147 * \f$ k \f$ is an eigenvector corresponding to eigenvalue number \f$ k 00148 * \f$ as returned by eigenvalues(). The eigenvectors are normalized to 00149 * have (Euclidean) norm equal to one. The matrix returned by this 00150 * function is the matrix \f$ V \f$ in the eigendecomposition \f$ A = V D 00151 * V^{-1} \f$, if it exists. 00152 * 00153 * Example: \include ComplexEigenSolver_eigenvectors.cpp 00154 * Output: \verbinclude ComplexEigenSolver_eigenvectors.out 00155 */ 00156 const EigenvectorType & eigenvectors() const 00157 { 00158 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 00159 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 00160 return m_eivec; 00161 } 00162 00163 /** \brief Returns the eigenvalues of given matrix. 00164 * 00165 * \returns A const reference to the column vector containing the eigenvalues. 00166 * 00167 * \pre Either the constructor 00168 * ComplexEigenSolver(const MatrixType& matrix, bool) or the member 00169 * function compute(const MatrixType& matrix, bool) has been called before 00170 * to compute the eigendecomposition of a matrix. 00171 * 00172 * This function returns a column vector containing the 00173 * eigenvalues. Eigenvalues are repeated according to their 00174 * algebraic multiplicity, so there are as many eigenvalues as 00175 * rows in the matrix. The eigenvalues are not sorted in any particular 00176 * order. 00177 * 00178 * Example: \include ComplexEigenSolver_eigenvalues.cpp 00179 * Output: \verbinclude ComplexEigenSolver_eigenvalues.out 00180 */ 00181 const EigenvalueType & eigenvalues() const 00182 { 00183 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 00184 return m_eivalues; 00185 } 00186 00187 /** \brief Computes eigendecomposition of given matrix. 00188 * 00189 * \param[in] matrix Square matrix whose eigendecomposition is to be computed. 00190 * \param[in] computeEigenvectors If true, both the eigenvectors and the 00191 * eigenvalues are computed; if false, only the eigenvalues are 00192 * computed. 00193 * \returns Reference to \c *this 00194 * 00195 * This function computes the eigenvalues of the complex matrix \p matrix. 00196 * The eigenvalues() function can be used to retrieve them. If 00197 * \p computeEigenvectors is true, then the eigenvectors are also computed 00198 * and can be retrieved by calling eigenvectors(). 00199 * 00200 * The matrix is first reduced to Schur form using the 00201 * ComplexSchur class. The Schur decomposition is then used to 00202 * compute the eigenvalues and eigenvectors. 00203 * 00204 * The cost of the computation is dominated by the cost of the 00205 * Schur decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$ 00206 * is the size of the matrix. 00207 * 00208 * Example: \include ComplexEigenSolver_compute.cpp 00209 * Output: \verbinclude ComplexEigenSolver_compute.out 00210 */ 00211 ComplexEigenSolver & compute(const MatrixType& matrix, bool computeEigenvectors = true); 00212 00213 /** \brief Reports whether previous computation was successful. 00214 * 00215 * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 00216 */ 00217 ComputationInfo info() const 00218 { 00219 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 00220 return m_schur.info(); 00221 } 00222 00223 /** \brief Sets the maximum number of iterations allowed. */ 00224 ComplexEigenSolver & setMaxIterations(Index maxIters) 00225 { 00226 m_schur.setMaxIterations(maxIters); 00227 return *this; 00228 } 00229 00230 /** \brief Returns the maximum number of iterations. */ 00231 Index getMaxIterations() 00232 { 00233 return m_schur.getMaxIterations(); 00234 } 00235 00236 protected: 00237 00238 static void check_template_parameters() 00239 { 00240 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00241 } 00242 00243 EigenvectorType m_eivec; 00244 EigenvalueType m_eivalues; 00245 ComplexSchur<MatrixType> m_schur; 00246 bool m_isInitialized; 00247 bool m_eigenvectorsOk; 00248 EigenvectorType m_matX; 00249 00250 private: 00251 void doComputeEigenvectors(const RealScalar& matrixnorm); 00252 void sortEigenvalues(bool computeEigenvectors); 00253 }; 00254 00255 00256 template<typename MatrixType> 00257 ComplexEigenSolver<MatrixType>& 00258 ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) 00259 { 00260 check_template_parameters(); 00261 00262 // this code is inspired from Jampack 00263 eigen_assert(matrix.cols() == matrix.rows()); 00264 00265 // Do a complex Schur decomposition, A = U T U^* 00266 // The eigenvalues are on the diagonal of T. 00267 m_schur.compute(matrix, computeEigenvectors); 00268 00269 if(m_schur.info() == Success) 00270 { 00271 m_eivalues = m_schur.matrixT().diagonal(); 00272 if(computeEigenvectors) 00273 doComputeEigenvectors(matrix.norm()); 00274 sortEigenvalues(computeEigenvectors); 00275 } 00276 00277 m_isInitialized = true; 00278 m_eigenvectorsOk = computeEigenvectors; 00279 return *this; 00280 } 00281 00282 00283 template<typename MatrixType> 00284 void ComplexEigenSolver<MatrixType>::doComputeEigenvectors (const RealScalar& matrixnorm) 00285 { 00286 const Index n = m_eivalues.size(); 00287 00288 // Compute X such that T = X D X^(-1), where D is the diagonal of T. 00289 // The matrix X is unit triangular. 00290 m_matX = EigenvectorType::Zero(n, n); 00291 for(Index k=n-1 ; k>=0 ; k--) 00292 { 00293 m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0); 00294 // Compute X(i,k) using the (i,k) entry of the equation X T = D X 00295 for(Index i=k-1 ; i>=0 ; i--) 00296 { 00297 m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k); 00298 if(k-i-1>0) 00299 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(); 00300 ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k); 00301 if(z==ComplexScalar(0)) 00302 { 00303 // If the i-th and k-th eigenvalue are equal, then z equals 0. 00304 // Use a small value instead, to prevent division by zero. 00305 numext::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm; 00306 } 00307 m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z; 00308 } 00309 } 00310 00311 // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1) 00312 m_eivec.noalias() = m_schur.matrixU() * m_matX; 00313 // .. and normalize the eigenvectors 00314 for(Index k=0 ; k<n ; k++) 00315 { 00316 m_eivec.col(k).normalize(); 00317 } 00318 } 00319 00320 00321 template<typename MatrixType> 00322 void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors) 00323 { 00324 const Index n = m_eivalues.size(); 00325 for (Index i=0; i<n; i++) 00326 { 00327 Index k; 00328 m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k); 00329 if (k != 0) 00330 { 00331 k += i; 00332 std::swap(m_eivalues[k],m_eivalues[i]); 00333 if(computeEigenvectors) 00334 m_eivec.col(i).swap(m_eivec.col(k)); 00335 } 00336 } 00337 } 00338 00339 } // end namespace Eigen 00340 00341 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H
Generated on Tue Jul 12 2022 17:46:51 by 1.7.2