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


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:36