ComplexSchur.h
Go to the documentation of this file.
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_SCHUR_H
00028 #define EIGEN_COMPLEX_SCHUR_H
00029 
00030 #include "./HessenbergDecomposition.h"
00031 
00032 namespace Eigen { 
00033 
00034 namespace internal {
00035 template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg;
00036 }
00037 
00066 template<typename _MatrixType> class ComplexSchur
00067 {
00068   public:
00069     typedef _MatrixType MatrixType;
00070     enum {
00071       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00072       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00073       Options = MatrixType::Options,
00074       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00075       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00076     };
00077 
00079     typedef typename MatrixType::Scalar Scalar;
00080     typedef typename NumTraits<Scalar>::Real RealScalar;
00081     typedef typename MatrixType::Index Index;
00082 
00089     typedef std::complex<RealScalar> ComplexScalar;
00090 
00096     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType;
00097 
00109     ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
00110       : m_matT(size,size),
00111         m_matU(size,size),
00112         m_hess(size),
00113         m_isInitialized(false),
00114         m_matUisUptodate(false)
00115     {}
00116 
00126     ComplexSchur(const MatrixType& matrix, bool computeU = true)
00127             : m_matT(matrix.rows(),matrix.cols()),
00128               m_matU(matrix.rows(),matrix.cols()),
00129               m_hess(matrix.rows()),
00130               m_isInitialized(false),
00131               m_matUisUptodate(false)
00132     {
00133       compute(matrix, computeU);
00134     }
00135 
00150     const ComplexMatrixType& matrixU() const
00151     {
00152       eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
00153       eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
00154       return m_matU;
00155     }
00156 
00174     const ComplexMatrixType& matrixT() const
00175     {
00176       eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
00177       return m_matT;
00178     }
00179 
00199     ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
00200 
00205     ComputationInfo info() const
00206     {
00207       eigen_assert(m_isInitialized && "RealSchur is not initialized.");
00208       return m_info;
00209     }
00210 
00215     static const int m_maxIterations = 30;
00216 
00217   protected:
00218     ComplexMatrixType m_matT, m_matU;
00219     HessenbergDecomposition<MatrixType> m_hess;
00220     ComputationInfo m_info;
00221     bool m_isInitialized;
00222     bool m_matUisUptodate;
00223 
00224   private:  
00225     bool subdiagonalEntryIsNeglegible(Index i);
00226     ComplexScalar computeShift(Index iu, Index iter);
00227     void reduceToTriangularForm(bool computeU);
00228     friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
00229 };
00230 
00234 template<typename MatrixType>
00235 inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
00236 {
00237   RealScalar d = internal::norm1(m_matT.coeff(i,i)) + internal::norm1(m_matT.coeff(i+1,i+1));
00238   RealScalar sd = internal::norm1(m_matT.coeff(i+1,i));
00239   if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
00240   {
00241     m_matT.coeffRef(i+1,i) = ComplexScalar(0);
00242     return true;
00243   }
00244   return false;
00245 }
00246 
00247 
00249 template<typename MatrixType>
00250 typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter)
00251 {
00252   if (iter == 10 || iter == 20) 
00253   {
00254     // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
00255     return internal::abs(internal::real(m_matT.coeff(iu,iu-1))) + internal::abs(internal::real(m_matT.coeff(iu-1,iu-2)));
00256   }
00257 
00258   // compute the shift as one of the eigenvalues of t, the 2x2
00259   // diagonal block on the bottom of the active submatrix
00260   Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
00261   RealScalar normt = t.cwiseAbs().sum();
00262   t /= normt;     // the normalization by sf is to avoid under/overflow
00263 
00264   ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
00265   ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
00266   ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
00267   ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
00268   ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
00269   ComplexScalar eival1 = (trace + disc) / RealScalar(2);
00270   ComplexScalar eival2 = (trace - disc) / RealScalar(2);
00271 
00272   if(internal::norm1(eival1) > internal::norm1(eival2))
00273     eival2 = det / eival1;
00274   else
00275     eival1 = det / eival2;
00276 
00277   // choose the eigenvalue closest to the bottom entry of the diagonal
00278   if(internal::norm1(eival1-t.coeff(1,1)) < internal::norm1(eival2-t.coeff(1,1)))
00279     return normt * eival1;
00280   else
00281     return normt * eival2;
00282 }
00283 
00284 
00285 template<typename MatrixType>
00286 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
00287 {
00288   m_matUisUptodate = false;
00289   eigen_assert(matrix.cols() == matrix.rows());
00290 
00291   if(matrix.cols() == 1)
00292   {
00293     m_matT = matrix.template cast<ComplexScalar>();
00294     if(computeU)  m_matU = ComplexMatrixType::Identity(1,1);
00295     m_info = Success;
00296     m_isInitialized = true;
00297     m_matUisUptodate = computeU;
00298     return *this;
00299   }
00300 
00301   internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU);
00302   reduceToTriangularForm(computeU);
00303   return *this;
00304 }
00305 
00306 namespace internal {
00307 
00308 /* Reduce given matrix to Hessenberg form */
00309 template<typename MatrixType, bool IsComplex>
00310 struct complex_schur_reduce_to_hessenberg
00311 {
00312   // this is the implementation for the case IsComplex = true
00313   static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
00314   {
00315     _this.m_hess.compute(matrix);
00316     _this.m_matT = _this.m_hess.matrixH();
00317     if(computeU)  _this.m_matU = _this.m_hess.matrixQ();
00318   }
00319 };
00320 
00321 template<typename MatrixType>
00322 struct complex_schur_reduce_to_hessenberg<MatrixType, false>
00323 {
00324   static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
00325   {
00326     typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
00327     typedef typename ComplexSchur<MatrixType>::ComplexMatrixType ComplexMatrixType;
00328 
00329     // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
00330     _this.m_hess.compute(matrix);
00331     _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
00332     if(computeU)  
00333     {
00334       // This may cause an allocation which seems to be avoidable
00335       MatrixType Q = _this.m_hess.matrixQ(); 
00336       _this.m_matU = Q.template cast<ComplexScalar>();
00337     }
00338   }
00339 };
00340 
00341 } // end namespace internal
00342 
00343 // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
00344 template<typename MatrixType>
00345 void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
00346 {  
00347   // The matrix m_matT is divided in three parts. 
00348   // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 
00349   // Rows il,...,iu is the part we are working on (the active submatrix).
00350   // Rows iu+1,...,end are already brought in triangular form.
00351   Index iu = m_matT.cols() - 1;
00352   Index il;
00353   Index iter = 0; // number of iterations we are working on the (iu,iu) element
00354 
00355   while(true)
00356   {
00357     // find iu, the bottom row of the active submatrix
00358     while(iu > 0)
00359     {
00360       if(!subdiagonalEntryIsNeglegible(iu-1)) break;
00361       iter = 0;
00362       --iu;
00363     }
00364 
00365     // if iu is zero then we are done; the whole matrix is triangularized
00366     if(iu==0) break;
00367 
00368     // if we spent too many iterations on the current element, we give up
00369     iter++;
00370     if(iter > m_maxIterations) break;
00371 
00372     // find il, the top row of the active submatrix
00373     il = iu-1;
00374     while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
00375     {
00376       --il;
00377     }
00378 
00379     /* perform the QR step using Givens rotations. The first rotation
00380        creates a bulge; the (il+2,il) element becomes nonzero. This
00381        bulge is chased down to the bottom of the active submatrix. */
00382 
00383     ComplexScalar shift = computeShift(iu, iter);
00384     JacobiRotation<ComplexScalar> rot;
00385     rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
00386     m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
00387     m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
00388     if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
00389 
00390     for(Index i=il+1 ; i<iu ; i++)
00391     {
00392       rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
00393       m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
00394       m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
00395       m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
00396       if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
00397     }
00398   }
00399 
00400   if(iter <= m_maxIterations) 
00401     m_info = Success;
00402   else
00403     m_info = NoConvergence;
00404 
00405   m_isInitialized = true;
00406   m_matUisUptodate = computeU;
00407 }
00408 
00409 } // end namespace Eigen
00410 
00411 #endif // EIGEN_COMPLEX_SCHUR_H