$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 typename MatrixType::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 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 RealSchur(const MatrixType& matrix, bool computeU = true) 00104 : m_matT(matrix.rows(),matrix.cols()), 00105 m_matU(matrix.rows(),matrix.cols()), 00106 m_workspaceVector(matrix.rows()), 00107 m_hess(matrix.rows()), 00108 m_isInitialized(false), 00109 m_matUisUptodate(false), 00110 m_maxIters(-1) 00111 { 00112 compute(matrix, computeU); 00113 } 00114 00126 const MatrixType& matrixU() const 00127 { 00128 eigen_assert(m_isInitialized && "RealSchur is not initialized."); 00129 eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the RealSchur decomposition."); 00130 return m_matU; 00131 } 00132 00143 const MatrixType& matrixT() const 00144 { 00145 eigen_assert(m_isInitialized && "RealSchur is not initialized."); 00146 return m_matT; 00147 } 00148 00168 RealSchur& compute(const MatrixType& matrix, bool computeU = true); 00169 00187 template<typename HessMatrixType, typename OrthMatrixType> 00188 RealSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU); 00193 ComputationInfo info() const 00194 { 00195 eigen_assert(m_isInitialized && "RealSchur is not initialized."); 00196 return m_info; 00197 } 00198 00204 RealSchur& setMaxIterations(Index maxIters) 00205 { 00206 m_maxIters = maxIters; 00207 return *this; 00208 } 00209 00211 Index getMaxIterations() 00212 { 00213 return m_maxIters; 00214 } 00215 00221 static const int m_maxIterationsPerRow = 40; 00222 00223 private: 00224 00225 MatrixType m_matT; 00226 MatrixType m_matU; 00227 ColumnVectorType m_workspaceVector; 00228 HessenbergDecomposition<MatrixType> m_hess; 00229 ComputationInfo m_info; 00230 bool m_isInitialized; 00231 bool m_matUisUptodate; 00232 Index m_maxIters; 00233 00234 typedef Matrix<Scalar,3,1> Vector3s; 00235 00236 Scalar computeNormOfT(); 00237 Index findSmallSubdiagEntry(Index iu); 00238 void splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift); 00239 void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo); 00240 void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector); 00241 void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace); 00242 }; 00243 00244 00245 template<typename MatrixType> 00246 RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU) 00247 { 00248 eigen_assert(matrix.cols() == matrix.rows()); 00249 Index maxIters = m_maxIters; 00250 if (maxIters == -1) 00251 maxIters = m_maxIterationsPerRow * matrix.rows(); 00252 00253 // Step 1. Reduce to Hessenberg form 00254 m_hess.compute(matrix); 00255 00256 // Step 2. Reduce to real Schur form 00257 computeFromHessenberg(m_hess.matrixH(), m_hess.matrixQ(), computeU); 00258 00259 return *this; 00260 } 00261 template<typename MatrixType> 00262 template<typename HessMatrixType, typename OrthMatrixType> 00263 RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU) 00264 { 00265 m_matT = matrixH; 00266 if(computeU) 00267 m_matU = matrixQ; 00268 00269 Index maxIters = m_maxIters; 00270 if (maxIters == -1) 00271 maxIters = m_maxIterationsPerRow * matrixH.rows(); 00272 m_workspaceVector.resize(m_matT.cols()); 00273 Scalar* workspace = &m_workspaceVector.coeffRef(0); 00274 00275 // The matrix m_matT is divided in three parts. 00276 // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 00277 // Rows il,...,iu is the part we are working on (the active window). 00278 // Rows iu+1,...,end are already brought in triangular form. 00279 Index iu = m_matT.cols() - 1; 00280 Index iter = 0; // iteration count for current eigenvalue 00281 Index totalIter = 0; // iteration count for whole matrix 00282 Scalar exshift(0); // sum of exceptional shifts 00283 Scalar norm = computeNormOfT(); 00284 00285 if(norm!=0) 00286 { 00287 while (iu >= 0) 00288 { 00289 Index il = findSmallSubdiagEntry(iu); 00290 00291 // Check for convergence 00292 if (il == iu) // One root found 00293 { 00294 m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift; 00295 if (iu > 0) 00296 m_matT.coeffRef(iu, iu-1) = Scalar(0); 00297 iu--; 00298 iter = 0; 00299 } 00300 else if (il == iu-1) // Two roots found 00301 { 00302 splitOffTwoRows(iu, computeU, exshift); 00303 iu -= 2; 00304 iter = 0; 00305 } 00306 else // No convergence yet 00307 { 00308 // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG ) 00309 Vector3s firstHouseholderVector(0,0,0), shiftInfo; 00310 computeShift(iu, iter, exshift, shiftInfo); 00311 iter = iter + 1; 00312 totalIter = totalIter + 1; 00313 if (totalIter > maxIters) break; 00314 Index im; 00315 initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector); 00316 performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace); 00317 } 00318 } 00319 } 00320 if(totalIter <= maxIters) 00321 m_info = Success; 00322 else 00323 m_info = NoConvergence; 00324 00325 m_isInitialized = true; 00326 m_matUisUptodate = computeU; 00327 return *this; 00328 } 00329 00331 template<typename MatrixType> 00332 inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT() 00333 { 00334 const Index size = m_matT.cols(); 00335 // FIXME to be efficient the following would requires a triangular reduxion code 00336 // Scalar norm = m_matT.upper().cwiseAbs().sum() 00337 // + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum(); 00338 Scalar norm(0); 00339 for (Index j = 0; j < size; ++j) 00340 norm += m_matT.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum(); 00341 return norm; 00342 } 00343 00345 template<typename MatrixType> 00346 inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu) 00347 { 00348 using std::abs; 00349 Index res = iu; 00350 while (res > 0) 00351 { 00352 Scalar s = abs(m_matT.coeff(res-1,res-1)) + abs(m_matT.coeff(res,res)); 00353 if (abs(m_matT.coeff(res,res-1)) <= NumTraits<Scalar>::epsilon() * s) 00354 break; 00355 res--; 00356 } 00357 return res; 00358 } 00359 00361 template<typename MatrixType> 00362 inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift) 00363 { 00364 using std::sqrt; 00365 using std::abs; 00366 const Index size = m_matT.cols(); 00367 00368 // The eigenvalues of the 2x2 matrix [a b; c d] are 00369 // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc 00370 Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu)); 00371 Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); // q = tr^2 / 4 - det = discr/4 00372 m_matT.coeffRef(iu,iu) += exshift; 00373 m_matT.coeffRef(iu-1,iu-1) += exshift; 00374 00375 if (q >= Scalar(0)) // Two real eigenvalues 00376 { 00377 Scalar z = sqrt(abs(q)); 00378 JacobiRotation<Scalar> rot; 00379 if (p >= Scalar(0)) 00380 rot.makeGivens(p + z, m_matT.coeff(iu, iu-1)); 00381 else 00382 rot.makeGivens(p - z, m_matT.coeff(iu, iu-1)); 00383 00384 m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint()); 00385 m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot); 00386 m_matT.coeffRef(iu, iu-1) = Scalar(0); 00387 if (computeU) 00388 m_matU.applyOnTheRight(iu-1, iu, rot); 00389 } 00390 00391 if (iu > 1) 00392 m_matT.coeffRef(iu-1, iu-2) = Scalar(0); 00393 } 00394 00396 template<typename MatrixType> 00397 inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo) 00398 { 00399 using std::sqrt; 00400 using std::abs; 00401 shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu); 00402 shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1); 00403 shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); 00404 00405 // Wilkinson's original ad hoc shift 00406 if (iter == 10) 00407 { 00408 exshift += shiftInfo.coeff(0); 00409 for (Index i = 0; i <= iu; ++i) 00410 m_matT.coeffRef(i,i) -= shiftInfo.coeff(0); 00411 Scalar s = abs(m_matT.coeff(iu,iu-1)) + abs(m_matT.coeff(iu-1,iu-2)); 00412 shiftInfo.coeffRef(0) = Scalar(0.75) * s; 00413 shiftInfo.coeffRef(1) = Scalar(0.75) * s; 00414 shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s; 00415 } 00416 00417 // MATLAB's new ad hoc shift 00418 if (iter == 30) 00419 { 00420 Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0); 00421 s = s * s + shiftInfo.coeff(2); 00422 if (s > Scalar(0)) 00423 { 00424 s = sqrt(s); 00425 if (shiftInfo.coeff(1) < shiftInfo.coeff(0)) 00426 s = -s; 00427 s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0); 00428 s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s; 00429 exshift += s; 00430 for (Index i = 0; i <= iu; ++i) 00431 m_matT.coeffRef(i,i) -= s; 00432 shiftInfo.setConstant(Scalar(0.964)); 00433 } 00434 } 00435 } 00436 00438 template<typename MatrixType> 00439 inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector) 00440 { 00441 using std::abs; 00442 Vector3s& v = firstHouseholderVector; // alias to save typing 00443 00444 for (im = iu-2; im >= il; --im) 00445 { 00446 const Scalar Tmm = m_matT.coeff(im,im); 00447 const Scalar r = shiftInfo.coeff(0) - Tmm; 00448 const Scalar s = shiftInfo.coeff(1) - Tmm; 00449 v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1); 00450 v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s; 00451 v.coeffRef(2) = m_matT.coeff(im+2,im+1); 00452 if (im == il) { 00453 break; 00454 } 00455 const Scalar lhs = m_matT.coeff(im,im-1) * (abs(v.coeff(1)) + abs(v.coeff(2))); 00456 const Scalar rhs = v.coeff(0) * (abs(m_matT.coeff(im-1,im-1)) + abs(Tmm) + abs(m_matT.coeff(im+1,im+1))); 00457 if (abs(lhs) < NumTraits<Scalar>::epsilon() * rhs) 00458 break; 00459 } 00460 } 00461 00463 template<typename MatrixType> 00464 inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace) 00465 { 00466 eigen_assert(im >= il); 00467 eigen_assert(im <= iu-2); 00468 00469 const Index size = m_matT.cols(); 00470 00471 for (Index k = im; k <= iu-2; ++k) 00472 { 00473 bool firstIteration = (k == im); 00474 00475 Vector3s v; 00476 if (firstIteration) 00477 v = firstHouseholderVector; 00478 else 00479 v = m_matT.template block<3,1>(k,k-1); 00480 00481 Scalar tau, beta; 00482 Matrix<Scalar, 2, 1> ess; 00483 v.makeHouseholder(ess, tau, beta); 00484 00485 if (beta != Scalar(0)) // if v is not zero 00486 { 00487 if (firstIteration && k > il) 00488 m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1); 00489 else if (!firstIteration) 00490 m_matT.coeffRef(k,k-1) = beta; 00491 00492 // These Householder transformations form the O(n^3) part of the algorithm 00493 m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace); 00494 m_matT.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace); 00495 if (computeU) 00496 m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace); 00497 } 00498 } 00499 00500 Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2); 00501 Scalar tau, beta; 00502 Matrix<Scalar, 1, 1> ess; 00503 v.makeHouseholder(ess, tau, beta); 00504 00505 if (beta != Scalar(0)) // if v is not zero 00506 { 00507 m_matT.coeffRef(iu-1, iu-2) = beta; 00508 m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace); 00509 m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace); 00510 if (computeU) 00511 m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace); 00512 } 00513 00514 // clean up pollution due to round-off errors 00515 for (Index i = im+2; i <= iu; ++i) 00516 { 00517 m_matT.coeffRef(i,i-2) = Scalar(0); 00518 if (i > im+2) 00519 m_matT.coeffRef(i,i-3) = Scalar(0); 00520 } 00521 } 00522 00523 } // end namespace Eigen 00524 00525 #endif // EIGEN_REAL_SCHUR_H