Eigen  3.3.3
ComplexSchur.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_SCHUR_H
00013 #define EIGEN_COMPLEX_SCHUR_H
00014 
00015 #include "./HessenbergDecomposition.h"
00016 
00017 namespace Eigen { 
00018 
00019 namespace internal {
00020 template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg;
00021 }
00022 
00051 template<typename _MatrixType> class ComplexSchur
00052 {
00053   public:
00054     typedef _MatrixType MatrixType;
00055     enum {
00056       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00057       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00058       Options = MatrixType::Options,
00059       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00060       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00061     };
00062 
00064     typedef typename MatrixType::Scalar Scalar;
00065     typedef typename NumTraits<Scalar>::Real RealScalar;
00066     typedef Eigen::Index Index; 
00067 
00074     typedef std::complex<RealScalar> ComplexScalar;
00075 
00081     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType;
00082 
00094     explicit ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
00095       : m_matT(size,size),
00096         m_matU(size,size),
00097         m_hess(size),
00098         m_isInitialized(false),
00099         m_matUisUptodate(false),
00100         m_maxIters(-1)
00101     {}
00102 
00112     template<typename InputType>
00113     explicit ComplexSchur(const EigenBase<InputType>& matrix, bool computeU = true)
00114       : m_matT(matrix.rows(),matrix.cols()),
00115         m_matU(matrix.rows(),matrix.cols()),
00116         m_hess(matrix.rows()),
00117         m_isInitialized(false),
00118         m_matUisUptodate(false),
00119         m_maxIters(-1)
00120     {
00121       compute(matrix.derived(), computeU);
00122     }
00123 
00138     const ComplexMatrixType& matrixU() const
00139     {
00140       eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
00141       eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
00142       return m_matU;
00143     }
00144 
00162     const ComplexMatrixType& matrixT() const
00163     {
00164       eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
00165       return m_matT;
00166     }
00167 
00190     template<typename InputType>
00191     ComplexSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true);
00192     
00210     template<typename HessMatrixType, typename OrthMatrixType>
00211     ComplexSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ,  bool computeU=true);
00212 
00217     ComputationInfo info() const
00218     {
00219       eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
00220       return m_info;
00221     }
00222 
00228     ComplexSchur& setMaxIterations(Index maxIters)
00229     {
00230       m_maxIters = maxIters;
00231       return *this;
00232     }
00233 
00235     Index getMaxIterations()
00236     {
00237       return m_maxIters;
00238     }
00239 
00245     static const int m_maxIterationsPerRow = 30;
00246 
00247   protected:
00248     ComplexMatrixType m_matT, m_matU;
00249     HessenbergDecomposition<MatrixType> m_hess;
00250     ComputationInfo m_info;
00251     bool m_isInitialized;
00252     bool m_matUisUptodate;
00253     Index m_maxIters;
00254 
00255   private:  
00256     bool subdiagonalEntryIsNeglegible(Index i);
00257     ComplexScalar computeShift(Index iu, Index iter);
00258     void reduceToTriangularForm(bool computeU);
00259     friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
00260 };
00261 
00265 template<typename MatrixType>
00266 inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
00267 {
00268   RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1));
00269   RealScalar sd = numext::norm1(m_matT.coeff(i+1,i));
00270   if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
00271   {
00272     m_matT.coeffRef(i+1,i) = ComplexScalar(0);
00273     return true;
00274   }
00275   return false;
00276 }
00277 
00278 
00280 template<typename MatrixType>
00281 typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter)
00282 {
00283   using std::abs;
00284   if (iter == 10 || iter == 20) 
00285   {
00286     // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
00287     return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2)));
00288   }
00289 
00290   // compute the shift as one of the eigenvalues of t, the 2x2
00291   // diagonal block on the bottom of the active submatrix
00292   Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
00293   RealScalar normt = t.cwiseAbs().sum();
00294   t /= normt;     // the normalization by sf is to avoid under/overflow
00295 
00296   ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
00297   ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
00298   ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
00299   ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
00300   ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
00301   ComplexScalar eival1 = (trace + disc) / RealScalar(2);
00302   ComplexScalar eival2 = (trace - disc) / RealScalar(2);
00303 
00304   if(numext::norm1(eival1) > numext::norm1(eival2))
00305     eival2 = det / eival1;
00306   else
00307     eival1 = det / eival2;
00308 
00309   // choose the eigenvalue closest to the bottom entry of the diagonal
00310   if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1)))
00311     return normt * eival1;
00312   else
00313     return normt * eival2;
00314 }
00315 
00316 
00317 template<typename MatrixType>
00318 template<typename InputType>
00319 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU)
00320 {
00321   m_matUisUptodate = false;
00322   eigen_assert(matrix.cols() == matrix.rows());
00323 
00324   if(matrix.cols() == 1)
00325   {
00326     m_matT = matrix.derived().template cast<ComplexScalar>();
00327     if(computeU)  m_matU = ComplexMatrixType::Identity(1,1);
00328     m_info = Success;
00329     m_isInitialized = true;
00330     m_matUisUptodate = computeU;
00331     return *this;
00332   }
00333 
00334   internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix.derived(), computeU);
00335   computeFromHessenberg(m_matT, m_matU, computeU);
00336   return *this;
00337 }
00338 
00339 template<typename MatrixType>
00340 template<typename HessMatrixType, typename OrthMatrixType>
00341 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
00342 {
00343   m_matT = matrixH;
00344   if(computeU)
00345     m_matU = matrixQ;
00346   reduceToTriangularForm(computeU);
00347   return *this;
00348 }
00349 namespace internal {
00350 
00351 /* Reduce given matrix to Hessenberg form */
00352 template<typename MatrixType, bool IsComplex>
00353 struct complex_schur_reduce_to_hessenberg
00354 {
00355   // this is the implementation for the case IsComplex = true
00356   static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
00357   {
00358     _this.m_hess.compute(matrix);
00359     _this.m_matT = _this.m_hess.matrixH();
00360     if(computeU)  _this.m_matU = _this.m_hess.matrixQ();
00361   }
00362 };
00363 
00364 template<typename MatrixType>
00365 struct complex_schur_reduce_to_hessenberg<MatrixType, false>
00366 {
00367   static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
00368   {
00369     typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
00370 
00371     // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
00372     _this.m_hess.compute(matrix);
00373     _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
00374     if(computeU)  
00375     {
00376       // This may cause an allocation which seems to be avoidable
00377       MatrixType Q = _this.m_hess.matrixQ(); 
00378       _this.m_matU = Q.template cast<ComplexScalar>();
00379     }
00380   }
00381 };
00382 
00383 } // end namespace internal
00384 
00385 // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
00386 template<typename MatrixType>
00387 void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
00388 {  
00389   Index maxIters = m_maxIters;
00390   if (maxIters == -1)
00391     maxIters = m_maxIterationsPerRow * m_matT.rows();
00392 
00393   // The matrix m_matT is divided in three parts. 
00394   // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 
00395   // Rows il,...,iu is the part we are working on (the active submatrix).
00396   // Rows iu+1,...,end are already brought in triangular form.
00397   Index iu = m_matT.cols() - 1;
00398   Index il;
00399   Index iter = 0; // number of iterations we are working on the (iu,iu) element
00400   Index totalIter = 0; // number of iterations for whole matrix
00401 
00402   while(true)
00403   {
00404     // find iu, the bottom row of the active submatrix
00405     while(iu > 0)
00406     {
00407       if(!subdiagonalEntryIsNeglegible(iu-1)) break;
00408       iter = 0;
00409       --iu;
00410     }
00411 
00412     // if iu is zero then we are done; the whole matrix is triangularized
00413     if(iu==0) break;
00414 
00415     // if we spent too many iterations, we give up
00416     iter++;
00417     totalIter++;
00418     if(totalIter > maxIters) break;
00419 
00420     // find il, the top row of the active submatrix
00421     il = iu-1;
00422     while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
00423     {
00424       --il;
00425     }
00426 
00427     /* perform the QR step using Givens rotations. The first rotation
00428        creates a bulge; the (il+2,il) element becomes nonzero. This
00429        bulge is chased down to the bottom of the active submatrix. */
00430 
00431     ComplexScalar shift = computeShift(iu, iter);
00432     JacobiRotation<ComplexScalar> rot;
00433     rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
00434     m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
00435     m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
00436     if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
00437 
00438     for(Index i=il+1 ; i<iu ; i++)
00439     {
00440       rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
00441       m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
00442       m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
00443       m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
00444       if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
00445     }
00446   }
00447 
00448   if(totalIter <= maxIters)
00449     m_info = Success;
00450   else
00451     m_info = NoConvergence;
00452 
00453   m_isInitialized = true;
00454   m_matUisUptodate = computeU;
00455 }
00456 
00457 } // end namespace Eigen
00458 
00459 #endif // EIGEN_COMPLEX_SCHUR_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends