![]() |
Eigen
3.3.3
|
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 Eigen::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 explicit 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 template<typename InputType> 00126 explicit ComplexEigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true) 00127 : m_eivec(matrix.rows(),matrix.cols()), 00128 m_eivalues(matrix.cols()), 00129 m_schur(matrix.rows()), 00130 m_isInitialized(false), 00131 m_eigenvectorsOk(false), 00132 m_matX(matrix.rows(),matrix.cols()) 00133 { 00134 compute(matrix.derived(), computeEigenvectors); 00135 } 00136 00157 const EigenvectorType& eigenvectors() const 00158 { 00159 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 00160 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 00161 return m_eivec; 00162 } 00163 00182 const EigenvalueType& eigenvalues() const 00183 { 00184 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 00185 return m_eivalues; 00186 } 00187 00212 template<typename InputType> 00213 ComplexEigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true); 00214 00219 ComputationInfo info() const 00220 { 00221 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 00222 return m_schur.info(); 00223 } 00224 00226 ComplexEigenSolver& setMaxIterations(Index maxIters) 00227 { 00228 m_schur.setMaxIterations(maxIters); 00229 return *this; 00230 } 00231 00233 Index getMaxIterations() 00234 { 00235 return m_schur.getMaxIterations(); 00236 } 00237 00238 protected: 00239 00240 static void check_template_parameters() 00241 { 00242 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00243 } 00244 00245 EigenvectorType m_eivec; 00246 EigenvalueType m_eivalues; 00247 ComplexSchur<MatrixType> m_schur; 00248 bool m_isInitialized; 00249 bool m_eigenvectorsOk; 00250 EigenvectorType m_matX; 00251 00252 private: 00253 void doComputeEigenvectors(RealScalar matrixnorm); 00254 void sortEigenvalues(bool computeEigenvectors); 00255 }; 00256 00257 00258 template<typename MatrixType> 00259 template<typename InputType> 00260 ComplexEigenSolver<MatrixType>& 00261 ComplexEigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors) 00262 { 00263 check_template_parameters(); 00264 00265 // this code is inspired from Jampack 00266 eigen_assert(matrix.cols() == matrix.rows()); 00267 00268 // Do a complex Schur decomposition, A = U T U^* 00269 // The eigenvalues are on the diagonal of T. 00270 m_schur.compute(matrix.derived(), computeEigenvectors); 00271 00272 if(m_schur.info() == Success) 00273 { 00274 m_eivalues = m_schur.matrixT().diagonal(); 00275 if(computeEigenvectors) 00276 doComputeEigenvectors(m_schur.matrixT().norm()); 00277 sortEigenvalues(computeEigenvectors); 00278 } 00279 00280 m_isInitialized = true; 00281 m_eigenvectorsOk = computeEigenvectors; 00282 return *this; 00283 } 00284 00285 00286 template<typename MatrixType> 00287 void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm) 00288 { 00289 const Index n = m_eivalues.size(); 00290 00291 matrixnorm = numext::maxi(matrixnorm,(std::numeric_limits<RealScalar>::min)()); 00292 00293 // Compute X such that T = X D X^(-1), where D is the diagonal of T. 00294 // The matrix X is unit triangular. 00295 m_matX = EigenvectorType::Zero(n, n); 00296 for(Index k=n-1 ; k>=0 ; k--) 00297 { 00298 m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0); 00299 // Compute X(i,k) using the (i,k) entry of the equation X T = D X 00300 for(Index i=k-1 ; i>=0 ; i--) 00301 { 00302 m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k); 00303 if(k-i-1>0) 00304 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(); 00305 ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k); 00306 if(z==ComplexScalar(0)) 00307 { 00308 // If the i-th and k-th eigenvalue are equal, then z equals 0. 00309 // Use a small value instead, to prevent division by zero. 00310 numext::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm; 00311 } 00312 m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z; 00313 } 00314 } 00315 00316 // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1) 00317 m_eivec.noalias() = m_schur.matrixU() * m_matX; 00318 // .. and normalize the eigenvectors 00319 for(Index k=0 ; k<n ; k++) 00320 { 00321 m_eivec.col(k).normalize(); 00322 } 00323 } 00324 00325 00326 template<typename MatrixType> 00327 void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors) 00328 { 00329 const Index n = m_eivalues.size(); 00330 for (Index i=0; i<n; i++) 00331 { 00332 Index k; 00333 m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k); 00334 if (k != 0) 00335 { 00336 k += i; 00337 std::swap(m_eivalues[k],m_eivalues[i]); 00338 if(computeEigenvectors) 00339 m_eivec.col(i).swap(m_eivec.col(k)); 00340 } 00341 } 00342 } 00343 00344 } // end namespace Eigen 00345 00346 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H