00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
00028 #define EIGEN_COMPLEX_EIGEN_SOLVER_H
00029
00030 #include "./ComplexSchur.h"
00031
00032 namespace Eigen {
00033
00060 template<typename _MatrixType> class ComplexEigenSolver
00061 {
00062 public:
00063
00065 typedef _MatrixType MatrixType;
00066
00067 enum {
00068 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00069 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00070 Options = MatrixType::Options,
00071 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00072 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00073 };
00074
00076 typedef typename MatrixType::Scalar Scalar;
00077 typedef typename NumTraits<Scalar>::Real RealScalar;
00078 typedef typename MatrixType::Index Index;
00079
00086 typedef std::complex<RealScalar> ComplexScalar;
00087
00093 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1> EigenvalueType;
00094
00100 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorType;
00101
00107 ComplexEigenSolver()
00108 : m_eivec(),
00109 m_eivalues(),
00110 m_schur(),
00111 m_isInitialized(false),
00112 m_eigenvectorsOk(false),
00113 m_matX()
00114 {}
00115
00122 ComplexEigenSolver(Index size)
00123 : m_eivec(size, size),
00124 m_eivalues(size),
00125 m_schur(size),
00126 m_isInitialized(false),
00127 m_eigenvectorsOk(false),
00128 m_matX(size, size)
00129 {}
00130
00140 ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
00141 : m_eivec(matrix.rows(),matrix.cols()),
00142 m_eivalues(matrix.cols()),
00143 m_schur(matrix.rows()),
00144 m_isInitialized(false),
00145 m_eigenvectorsOk(false),
00146 m_matX(matrix.rows(),matrix.cols())
00147 {
00148 compute(matrix, computeEigenvectors);
00149 }
00150
00171 const EigenvectorType& eigenvectors() const
00172 {
00173 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00174 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00175 return m_eivec;
00176 }
00177
00196 const EigenvalueType& eigenvalues() const
00197 {
00198 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00199 return m_eivalues;
00200 }
00201
00226 ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
00227
00232 ComputationInfo info() const
00233 {
00234 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00235 return m_schur.info();
00236 }
00237
00238 protected:
00239 EigenvectorType m_eivec;
00240 EigenvalueType m_eivalues;
00241 ComplexSchur<MatrixType> m_schur;
00242 bool m_isInitialized;
00243 bool m_eigenvectorsOk;
00244 EigenvectorType m_matX;
00245
00246 private:
00247 void doComputeEigenvectors(RealScalar matrixnorm);
00248 void sortEigenvalues(bool computeEigenvectors);
00249 };
00250
00251
00252 template<typename MatrixType>
00253 ComplexEigenSolver<MatrixType>& ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
00254 {
00255
00256 assert(matrix.cols() == matrix.rows());
00257
00258
00259
00260 m_schur.compute(matrix, computeEigenvectors);
00261
00262 if(m_schur.info() == Success)
00263 {
00264 m_eivalues = m_schur.matrixT().diagonal();
00265 if(computeEigenvectors)
00266 doComputeEigenvectors(matrix.norm());
00267 sortEigenvalues(computeEigenvectors);
00268 }
00269
00270 m_isInitialized = true;
00271 m_eigenvectorsOk = computeEigenvectors;
00272 return *this;
00273 }
00274
00275
00276 template<typename MatrixType>
00277 void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm)
00278 {
00279 const Index n = m_eivalues.size();
00280
00281
00282
00283 m_matX = EigenvectorType::Zero(n, n);
00284 for(Index k=n-1 ; k>=0 ; k--)
00285 {
00286 m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
00287
00288 for(Index i=k-1 ; i>=0 ; i--)
00289 {
00290 m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
00291 if(k-i-1>0)
00292 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();
00293 ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
00294 if(z==ComplexScalar(0))
00295 {
00296
00297
00298 internal::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
00299 }
00300 m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
00301 }
00302 }
00303
00304
00305 m_eivec.noalias() = m_schur.matrixU() * m_matX;
00306
00307 for(Index k=0 ; k<n ; k++)
00308 {
00309 m_eivec.col(k).normalize();
00310 }
00311 }
00312
00313
00314 template<typename MatrixType>
00315 void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors)
00316 {
00317 const Index n = m_eivalues.size();
00318 for (Index i=0; i<n; i++)
00319 {
00320 Index k;
00321 m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
00322 if (k != 0)
00323 {
00324 k += i;
00325 std::swap(m_eivalues[k],m_eivalues[i]);
00326 if(computeEigenvectors)
00327 m_eivec.col(i).swap(m_eivec.col(k));
00328 }
00329 }
00330 }
00331
00332 }
00333
00334 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H