![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_REAL_SCHUR_H 00012 #define EIGEN_REAL_SCHUR_H 00013 00014 #include "./HessenbergDecomposition.h" 00015 00016 namespace Eigen { 00017 00054 template<typename _MatrixType> class RealSchur 00055 { 00056 public: 00057 typedef _MatrixType MatrixType; 00058 enum { 00059 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00060 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00061 Options = MatrixType::Options, 00062 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00063 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00064 }; 00065 typedef typename MatrixType::Scalar Scalar; 00066 typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; 00067 typedef Eigen::Index Index; 00068 00069 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType; 00070 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; 00071 00083 explicit RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) 00084 : m_matT(size, size), 00085 m_matU(size, size), 00086 m_workspaceVector(size), 00087 m_hess(size), 00088 m_isInitialized(false), 00089 m_matUisUptodate(false), 00090 m_maxIters(-1) 00091 { } 00092 00103 template<typename InputType> 00104 explicit RealSchur(const EigenBase<InputType>& matrix, bool computeU = true) 00105 : m_matT(matrix.rows(),matrix.cols()), 00106 m_matU(matrix.rows(),matrix.cols()), 00107 m_workspaceVector(matrix.rows()), 00108 m_hess(matrix.rows()), 00109 m_isInitialized(false), 00110 m_matUisUptodate(false), 00111 m_maxIters(-1) 00112 { 00113 compute(matrix.derived(), computeU); 00114 } 00115 00127 const MatrixType& matrixU() const 00128 { 00129 eigen_assert(m_isInitialized && "RealSchur is not initialized."); 00130 eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the RealSchur decomposition."); 00131 return m_matU; 00132 } 00133 00144 const MatrixType& matrixT() const 00145 { 00146 eigen_assert(m_isInitialized && "RealSchur is not initialized."); 00147 return m_matT; 00148 } 00149 00169 template<typename InputType> 00170 RealSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true); 00171 00189 template<typename HessMatrixType, typename OrthMatrixType> 00190 RealSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU); 00195 ComputationInfo info() const 00196 { 00197 eigen_assert(m_isInitialized && "RealSchur is not initialized."); 00198 return m_info; 00199 } 00200 00206 RealSchur& setMaxIterations(Index maxIters) 00207 { 00208 m_maxIters = maxIters; 00209 return *this; 00210 } 00211 00213 Index getMaxIterations() 00214 { 00215 return m_maxIters; 00216 } 00217 00223 static const int m_maxIterationsPerRow = 40; 00224 00225 private: 00226 00227 MatrixType m_matT; 00228 MatrixType m_matU; 00229 ColumnVectorType m_workspaceVector; 00230 HessenbergDecomposition<MatrixType> m_hess; 00231 ComputationInfo m_info; 00232 bool m_isInitialized; 00233 bool m_matUisUptodate; 00234 Index m_maxIters; 00235 00236 typedef Matrix<Scalar,3,1> Vector3s; 00237 00238 Scalar computeNormOfT(); 00239 Index findSmallSubdiagEntry(Index iu); 00240 void splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift); 00241 void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo); 00242 void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector); 00243 void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace); 00244 }; 00245 00246 00247 template<typename MatrixType> 00248 template<typename InputType> 00249 RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU) 00250 { 00251 const Scalar considerAsZero = (std::numeric_limits<Scalar>::min)(); 00252 00253 eigen_assert(matrix.cols() == matrix.rows()); 00254 Index maxIters = m_maxIters; 00255 if (maxIters == -1) 00256 maxIters = m_maxIterationsPerRow * matrix.rows(); 00257 00258 Scalar scale = matrix.derived().cwiseAbs().maxCoeff(); 00259 if(scale<considerAsZero) 00260 { 00261 m_matT.setZero(matrix.rows(),matrix.cols()); 00262 if(computeU) 00263 m_matU.setIdentity(matrix.rows(),matrix.cols()); 00264 m_info = Success; 00265 m_isInitialized = true; 00266 m_matUisUptodate = computeU; 00267 return *this; 00268 } 00269 00270 // Step 1. Reduce to Hessenberg form 00271 m_hess.compute(matrix.derived()/scale); 00272 00273 // Step 2. Reduce to real Schur form 00274 computeFromHessenberg(m_hess.matrixH(), m_hess.matrixQ(), computeU); 00275 00276 m_matT *= scale; 00277 00278 return *this; 00279 } 00280 template<typename MatrixType> 00281 template<typename HessMatrixType, typename OrthMatrixType> 00282 RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU) 00283 { 00284 using std::abs; 00285 00286 m_matT = matrixH; 00287 if(computeU) 00288 m_matU = matrixQ; 00289 00290 Index maxIters = m_maxIters; 00291 if (maxIters == -1) 00292 maxIters = m_maxIterationsPerRow * matrixH.rows(); 00293 m_workspaceVector.resize(m_matT.cols()); 00294 Scalar* workspace = &m_workspaceVector.coeffRef(0); 00295 00296 // The matrix m_matT is divided in three parts. 00297 // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 00298 // Rows il,...,iu is the part we are working on (the active window). 00299 // Rows iu+1,...,end are already brought in triangular form. 00300 Index iu = m_matT.cols() - 1; 00301 Index iter = 0; // iteration count for current eigenvalue 00302 Index totalIter = 0; // iteration count for whole matrix 00303 Scalar exshift(0); // sum of exceptional shifts 00304 Scalar norm = computeNormOfT(); 00305 00306 if(norm!=0) 00307 { 00308 while (iu >= 0) 00309 { 00310 Index il = findSmallSubdiagEntry(iu); 00311 00312 // Check for convergence 00313 if (il == iu) // One root found 00314 { 00315 m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift; 00316 if (iu > 0) 00317 m_matT.coeffRef(iu, iu-1) = Scalar(0); 00318 iu--; 00319 iter = 0; 00320 } 00321 else if (il == iu-1) // Two roots found 00322 { 00323 splitOffTwoRows(iu, computeU, exshift); 00324 iu -= 2; 00325 iter = 0; 00326 } 00327 else // No convergence yet 00328 { 00329 // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG ) 00330 Vector3s firstHouseholderVector(0,0,0), shiftInfo; 00331 computeShift(iu, iter, exshift, shiftInfo); 00332 iter = iter + 1; 00333 totalIter = totalIter + 1; 00334 if (totalIter > maxIters) break; 00335 Index im; 00336 initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector); 00337 performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace); 00338 } 00339 } 00340 } 00341 if(totalIter <= maxIters) 00342 m_info = Success; 00343 else 00344 m_info = NoConvergence; 00345 00346 m_isInitialized = true; 00347 m_matUisUptodate = computeU; 00348 return *this; 00349 } 00350 00352 template<typename MatrixType> 00353 inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT() 00354 { 00355 const Index size = m_matT.cols(); 00356 // FIXME to be efficient the following would requires a triangular reduxion code 00357 // Scalar norm = m_matT.upper().cwiseAbs().sum() 00358 // + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum(); 00359 Scalar norm(0); 00360 for (Index j = 0; j < size; ++j) 00361 norm += m_matT.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum(); 00362 return norm; 00363 } 00364 00366 template<typename MatrixType> 00367 inline Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu) 00368 { 00369 using std::abs; 00370 Index res = iu; 00371 while (res > 0) 00372 { 00373 Scalar s = abs(m_matT.coeff(res-1,res-1)) + abs(m_matT.coeff(res,res)); 00374 if (abs(m_matT.coeff(res,res-1)) <= NumTraits<Scalar>::epsilon() * s) 00375 break; 00376 res--; 00377 } 00378 return res; 00379 } 00380 00382 template<typename MatrixType> 00383 inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift) 00384 { 00385 using std::sqrt; 00386 using std::abs; 00387 const Index size = m_matT.cols(); 00388 00389 // The eigenvalues of the 2x2 matrix [a b; c d] are 00390 // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc 00391 Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu)); 00392 Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); // q = tr^2 / 4 - det = discr/4 00393 m_matT.coeffRef(iu,iu) += exshift; 00394 m_matT.coeffRef(iu-1,iu-1) += exshift; 00395 00396 if (q >= Scalar(0)) // Two real eigenvalues 00397 { 00398 Scalar z = sqrt(abs(q)); 00399 JacobiRotation<Scalar> rot; 00400 if (p >= Scalar(0)) 00401 rot.makeGivens(p + z, m_matT.coeff(iu, iu-1)); 00402 else 00403 rot.makeGivens(p - z, m_matT.coeff(iu, iu-1)); 00404 00405 m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint()); 00406 m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot); 00407 m_matT.coeffRef(iu, iu-1) = Scalar(0); 00408 if (computeU) 00409 m_matU.applyOnTheRight(iu-1, iu, rot); 00410 } 00411 00412 if (iu > 1) 00413 m_matT.coeffRef(iu-1, iu-2) = Scalar(0); 00414 } 00415 00417 template<typename MatrixType> 00418 inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo) 00419 { 00420 using std::sqrt; 00421 using std::abs; 00422 shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu); 00423 shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1); 00424 shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); 00425 00426 // Wilkinson's original ad hoc shift 00427 if (iter == 10) 00428 { 00429 exshift += shiftInfo.coeff(0); 00430 for (Index i = 0; i <= iu; ++i) 00431 m_matT.coeffRef(i,i) -= shiftInfo.coeff(0); 00432 Scalar s = abs(m_matT.coeff(iu,iu-1)) + abs(m_matT.coeff(iu-1,iu-2)); 00433 shiftInfo.coeffRef(0) = Scalar(0.75) * s; 00434 shiftInfo.coeffRef(1) = Scalar(0.75) * s; 00435 shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s; 00436 } 00437 00438 // MATLAB's new ad hoc shift 00439 if (iter == 30) 00440 { 00441 Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0); 00442 s = s * s + shiftInfo.coeff(2); 00443 if (s > Scalar(0)) 00444 { 00445 s = sqrt(s); 00446 if (shiftInfo.coeff(1) < shiftInfo.coeff(0)) 00447 s = -s; 00448 s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0); 00449 s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s; 00450 exshift += s; 00451 for (Index i = 0; i <= iu; ++i) 00452 m_matT.coeffRef(i,i) -= s; 00453 shiftInfo.setConstant(Scalar(0.964)); 00454 } 00455 } 00456 } 00457 00459 template<typename MatrixType> 00460 inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector) 00461 { 00462 using std::abs; 00463 Vector3s& v = firstHouseholderVector; // alias to save typing 00464 00465 for (im = iu-2; im >= il; --im) 00466 { 00467 const Scalar Tmm = m_matT.coeff(im,im); 00468 const Scalar r = shiftInfo.coeff(0) - Tmm; 00469 const Scalar s = shiftInfo.coeff(1) - Tmm; 00470 v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1); 00471 v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s; 00472 v.coeffRef(2) = m_matT.coeff(im+2,im+1); 00473 if (im == il) { 00474 break; 00475 } 00476 const Scalar lhs = m_matT.coeff(im,im-1) * (abs(v.coeff(1)) + abs(v.coeff(2))); 00477 const Scalar rhs = v.coeff(0) * (abs(m_matT.coeff(im-1,im-1)) + abs(Tmm) + abs(m_matT.coeff(im+1,im+1))); 00478 if (abs(lhs) < NumTraits<Scalar>::epsilon() * rhs) 00479 break; 00480 } 00481 } 00482 00484 template<typename MatrixType> 00485 inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace) 00486 { 00487 eigen_assert(im >= il); 00488 eigen_assert(im <= iu-2); 00489 00490 const Index size = m_matT.cols(); 00491 00492 for (Index k = im; k <= iu-2; ++k) 00493 { 00494 bool firstIteration = (k == im); 00495 00496 Vector3s v; 00497 if (firstIteration) 00498 v = firstHouseholderVector; 00499 else 00500 v = m_matT.template block<3,1>(k,k-1); 00501 00502 Scalar tau, beta; 00503 Matrix<Scalar, 2, 1> ess; 00504 v.makeHouseholder(ess, tau, beta); 00505 00506 if (beta != Scalar(0)) // if v is not zero 00507 { 00508 if (firstIteration && k > il) 00509 m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1); 00510 else if (!firstIteration) 00511 m_matT.coeffRef(k,k-1) = beta; 00512 00513 // These Householder transformations form the O(n^3) part of the algorithm 00514 m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace); 00515 m_matT.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace); 00516 if (computeU) 00517 m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace); 00518 } 00519 } 00520 00521 Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2); 00522 Scalar tau, beta; 00523 Matrix<Scalar, 1, 1> ess; 00524 v.makeHouseholder(ess, tau, beta); 00525 00526 if (beta != Scalar(0)) // if v is not zero 00527 { 00528 m_matT.coeffRef(iu-1, iu-2) = beta; 00529 m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace); 00530 m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace); 00531 if (computeU) 00532 m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace); 00533 } 00534 00535 // clean up pollution due to round-off errors 00536 for (Index i = im+2; i <= iu; ++i) 00537 { 00538 m_matT.coeffRef(i,i-2) = Scalar(0); 00539 if (i > im+2) 00540 m_matT.coeffRef(i,i-3) = Scalar(0); 00541 } 00542 } 00543 00544 } // end namespace Eigen 00545 00546 #endif // EIGEN_REAL_SCHUR_H