![]() |
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_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