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