$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 00045 template<typename _MatrixType> class ComplexEigenSolver 00046 { 00047 public: 00048 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 00061 typedef typename MatrixType::Scalar Scalar; 00062 typedef typename NumTraits<Scalar>::Real RealScalar; 00063 typedef typename MatrixType::Index Index; 00064 00071 typedef std::complex<RealScalar> ComplexScalar; 00072 00078 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1> EigenvalueType; 00079 00085 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorType; 00086 00092 ComplexEigenSolver() 00093 : m_eivec(), 00094 m_eivalues(), 00095 m_schur(), 00096 m_isInitialized(false), 00097 m_eigenvectorsOk(false), 00098 m_matX() 00099 {} 00100 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 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 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 00181 const EigenvalueType& eigenvalues() const 00182 { 00183 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 00184 return m_eivalues; 00185 } 00186 00211 ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); 00212 00217 ComputationInfo info() const 00218 { 00219 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 00220 return m_schur.info(); 00221 } 00222 00224 ComplexEigenSolver& setMaxIterations(Index maxIters) 00225 { 00226 m_schur.setMaxIterations(maxIters); 00227 return *this; 00228 } 00229 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