Eigen  3.2.10
 
Loading...
Searching...
No Matches
RealSchur.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
5// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
6//
7// This Source Code Form is subject to the terms of the Mozilla
8// Public License v. 2.0. If a copy of the MPL was not distributed
9// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11#ifndef EIGEN_REAL_SCHUR_H
12#define EIGEN_REAL_SCHUR_H
13
14#include "./HessenbergDecomposition.h"
15
16namespace Eigen {
17
54template<typename _MatrixType> class RealSchur
55{
56 public:
57 typedef _MatrixType MatrixType;
58 enum {
59 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
60 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
61 Options = MatrixType::Options,
62 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
63 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
64 };
65 typedef typename MatrixType::Scalar Scalar;
66 typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
67 typedef typename MatrixType::Index Index;
68
69 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
70 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
71
83 RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
84 : m_matT(size, size),
85 m_matU(size, size),
86 m_workspaceVector(size),
87 m_hess(size),
88 m_isInitialized(false),
89 m_matUisUptodate(false),
90 m_maxIters(-1)
91 { }
92
103 RealSchur(const MatrixType& matrix, bool computeU = true)
104 : m_matT(matrix.rows(),matrix.cols()),
105 m_matU(matrix.rows(),matrix.cols()),
106 m_workspaceVector(matrix.rows()),
107 m_hess(matrix.rows()),
108 m_isInitialized(false),
109 m_matUisUptodate(false),
110 m_maxIters(-1)
111 {
112 compute(matrix, computeU);
113 }
114
126 const MatrixType& matrixU() const
127 {
128 eigen_assert(m_isInitialized && "RealSchur is not initialized.");
129 eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the RealSchur decomposition.");
130 return m_matU;
131 }
132
143 const MatrixType& matrixT() const
144 {
145 eigen_assert(m_isInitialized && "RealSchur is not initialized.");
146 return m_matT;
147 }
148
168 RealSchur& compute(const MatrixType& matrix, bool computeU = true);
169
187 template<typename HessMatrixType, typename OrthMatrixType>
188 RealSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU);
193 ComputationInfo info() const
194 {
195 eigen_assert(m_isInitialized && "RealSchur is not initialized.");
196 return m_info;
197 }
198
204 RealSchur& setMaxIterations(Index maxIters)
205 {
206 m_maxIters = maxIters;
207 return *this;
208 }
209
211 Index getMaxIterations()
212 {
213 return m_maxIters;
214 }
215
221 static const int m_maxIterationsPerRow = 40;
222
223 private:
224
225 MatrixType m_matT;
226 MatrixType m_matU;
227 ColumnVectorType m_workspaceVector;
228 HessenbergDecomposition<MatrixType> m_hess;
229 ComputationInfo m_info;
230 bool m_isInitialized;
231 bool m_matUisUptodate;
232 Index m_maxIters;
233
234 typedef Matrix<Scalar,3,1> Vector3s;
235
236 Scalar computeNormOfT();
237 Index findSmallSubdiagEntry(Index iu);
238 void splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift);
239 void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo);
240 void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector);
241 void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace);
242};
243
244
245template<typename MatrixType>
246RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
247{
248 eigen_assert(matrix.cols() == matrix.rows());
249 Index maxIters = m_maxIters;
250 if (maxIters == -1)
251 maxIters = m_maxIterationsPerRow * matrix.rows();
252
253 // Step 1. Reduce to Hessenberg form
254 m_hess.compute(matrix);
255
256 // Step 2. Reduce to real Schur form
257 computeFromHessenberg(m_hess.matrixH(), m_hess.matrixQ(), computeU);
258
259 return *this;
260}
261template<typename MatrixType>
262template<typename HessMatrixType, typename OrthMatrixType>
263RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
264{
265 m_matT = matrixH;
266 if(computeU)
267 m_matU = matrixQ;
268
269 Index maxIters = m_maxIters;
270 if (maxIters == -1)
271 maxIters = m_maxIterationsPerRow * matrixH.rows();
272 m_workspaceVector.resize(m_matT.cols());
273 Scalar* workspace = &m_workspaceVector.coeffRef(0);
274
275 // The matrix m_matT is divided in three parts.
276 // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
277 // Rows il,...,iu is the part we are working on (the active window).
278 // Rows iu+1,...,end are already brought in triangular form.
279 Index iu = m_matT.cols() - 1;
280 Index iter = 0; // iteration count for current eigenvalue
281 Index totalIter = 0; // iteration count for whole matrix
282 Scalar exshift(0); // sum of exceptional shifts
283 Scalar norm = computeNormOfT();
284
285 if(norm!=0)
286 {
287 while (iu >= 0)
288 {
289 Index il = findSmallSubdiagEntry(iu);
290
291 // Check for convergence
292 if (il == iu) // One root found
293 {
294 m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
295 if (iu > 0)
296 m_matT.coeffRef(iu, iu-1) = Scalar(0);
297 iu--;
298 iter = 0;
299 }
300 else if (il == iu-1) // Two roots found
301 {
302 splitOffTwoRows(iu, computeU, exshift);
303 iu -= 2;
304 iter = 0;
305 }
306 else // No convergence yet
307 {
308 // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG )
309 Vector3s firstHouseholderVector(0,0,0), shiftInfo;
310 computeShift(iu, iter, exshift, shiftInfo);
311 iter = iter + 1;
312 totalIter = totalIter + 1;
313 if (totalIter > maxIters) break;
314 Index im;
315 initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
316 performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
317 }
318 }
319 }
320 if(totalIter <= maxIters)
321 m_info = Success;
322 else
323 m_info = NoConvergence;
324
325 m_isInitialized = true;
326 m_matUisUptodate = computeU;
327 return *this;
328}
329
331template<typename MatrixType>
332inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
333{
334 const Index size = m_matT.cols();
335 // FIXME to be efficient the following would requires a triangular reduxion code
336 // Scalar norm = m_matT.upper().cwiseAbs().sum()
337 // + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
338 Scalar norm(0);
339 for (Index j = 0; j < size; ++j)
340 norm += m_matT.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum();
341 return norm;
342}
343
345template<typename MatrixType>
346inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu)
347{
348 using std::abs;
349 Index res = iu;
350 while (res > 0)
351 {
352 Scalar s = abs(m_matT.coeff(res-1,res-1)) + abs(m_matT.coeff(res,res));
353 if (abs(m_matT.coeff(res,res-1)) <= NumTraits<Scalar>::epsilon() * s)
354 break;
355 res--;
356 }
357 return res;
358}
359
361template<typename MatrixType>
362inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift)
363{
364 using std::sqrt;
365 using std::abs;
366 const Index size = m_matT.cols();
367
368 // The eigenvalues of the 2x2 matrix [a b; c d] are
369 // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc
370 Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
371 Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); // q = tr^2 / 4 - det = discr/4
372 m_matT.coeffRef(iu,iu) += exshift;
373 m_matT.coeffRef(iu-1,iu-1) += exshift;
374
375 if (q >= Scalar(0)) // Two real eigenvalues
376 {
377 Scalar z = sqrt(abs(q));
378 JacobiRotation<Scalar> rot;
379 if (p >= Scalar(0))
380 rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
381 else
382 rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
383
384 m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
385 m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot);
386 m_matT.coeffRef(iu, iu-1) = Scalar(0);
387 if (computeU)
388 m_matU.applyOnTheRight(iu-1, iu, rot);
389 }
390
391 if (iu > 1)
392 m_matT.coeffRef(iu-1, iu-2) = Scalar(0);
393}
394
396template<typename MatrixType>
397inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo)
398{
399 using std::sqrt;
400 using std::abs;
401 shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
402 shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
403 shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
404
405 // Wilkinson's original ad hoc shift
406 if (iter == 10)
407 {
408 exshift += shiftInfo.coeff(0);
409 for (Index i = 0; i <= iu; ++i)
410 m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
411 Scalar s = abs(m_matT.coeff(iu,iu-1)) + abs(m_matT.coeff(iu-1,iu-2));
412 shiftInfo.coeffRef(0) = Scalar(0.75) * s;
413 shiftInfo.coeffRef(1) = Scalar(0.75) * s;
414 shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
415 }
416
417 // MATLAB's new ad hoc shift
418 if (iter == 30)
419 {
420 Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
421 s = s * s + shiftInfo.coeff(2);
422 if (s > Scalar(0))
423 {
424 s = sqrt(s);
425 if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
426 s = -s;
427 s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
428 s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
429 exshift += s;
430 for (Index i = 0; i <= iu; ++i)
431 m_matT.coeffRef(i,i) -= s;
432 shiftInfo.setConstant(Scalar(0.964));
433 }
434 }
435}
436
438template<typename MatrixType>
439inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector)
440{
441 using std::abs;
442 Vector3s& v = firstHouseholderVector; // alias to save typing
443
444 for (im = iu-2; im >= il; --im)
445 {
446 const Scalar Tmm = m_matT.coeff(im,im);
447 const Scalar r = shiftInfo.coeff(0) - Tmm;
448 const Scalar s = shiftInfo.coeff(1) - Tmm;
449 v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
450 v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s;
451 v.coeffRef(2) = m_matT.coeff(im+2,im+1);
452 if (im == il) {
453 break;
454 }
455 const Scalar lhs = m_matT.coeff(im,im-1) * (abs(v.coeff(1)) + abs(v.coeff(2)));
456 const Scalar rhs = v.coeff(0) * (abs(m_matT.coeff(im-1,im-1)) + abs(Tmm) + abs(m_matT.coeff(im+1,im+1)));
457 if (abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
458 break;
459 }
460}
461
463template<typename MatrixType>
464inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace)
465{
466 eigen_assert(im >= il);
467 eigen_assert(im <= iu-2);
468
469 const Index size = m_matT.cols();
470
471 for (Index k = im; k <= iu-2; ++k)
472 {
473 bool firstIteration = (k == im);
474
475 Vector3s v;
476 if (firstIteration)
477 v = firstHouseholderVector;
478 else
479 v = m_matT.template block<3,1>(k,k-1);
480
481 Scalar tau, beta;
482 Matrix<Scalar, 2, 1> ess;
483 v.makeHouseholder(ess, tau, beta);
484
485 if (beta != Scalar(0)) // if v is not zero
486 {
487 if (firstIteration && k > il)
488 m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
489 else if (!firstIteration)
490 m_matT.coeffRef(k,k-1) = beta;
491
492 // These Householder transformations form the O(n^3) part of the algorithm
493 m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
494 m_matT.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
495 if (computeU)
496 m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
497 }
498 }
499
500 Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2);
501 Scalar tau, beta;
502 Matrix<Scalar, 1, 1> ess;
503 v.makeHouseholder(ess, tau, beta);
504
505 if (beta != Scalar(0)) // if v is not zero
506 {
507 m_matT.coeffRef(iu-1, iu-2) = beta;
508 m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
509 m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
510 if (computeU)
511 m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
512 }
513
514 // clean up pollution due to round-off errors
515 for (Index i = im+2; i <= iu; ++i)
516 {
517 m_matT.coeffRef(i,i-2) = Scalar(0);
518 if (i > im+2)
519 m_matT.coeffRef(i,i-3) = Scalar(0);
520 }
521}
522
523} // end namespace Eigen
524
525#endif // EIGEN_REAL_SCHUR_H
Eigen::HessenbergDecomposition
Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation.
Definition HessenbergDecomposition.h:58
Eigen::JacobiRotation
Rotation given by a cosine-sine pair.
Definition Jacobi.h:35
Eigen::JacobiRotation::makeGivens
void makeGivens(const Scalar &p, const Scalar &q, Scalar *z=0)
Definition Jacobi.h:148
Eigen::MatrixBase::makeHouseholder
void makeHouseholder(EssentialPart &essential, Scalar &tau, RealScalar &beta) const
Definition Householder.h:65
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition Matrix.h:129
Eigen::RealSchur
Performs a real Schur decomposition of a square matrix.
Definition RealSchur.h:55
Eigen::RealSchur::computeFromHessenberg
RealSchur & computeFromHessenberg(const HessMatrixType &matrixH, const OrthMatrixType &matrixQ, bool computeU)
Computes Schur decomposition of a Hessenberg matrix H = Z T Z^T.
Eigen::RealSchur::RealSchur
RealSchur(const MatrixType &matrix, bool computeU=true)
Constructor; computes real Schur decomposition of given matrix.
Definition RealSchur.h:103
Eigen::RealSchur::matrixU
const MatrixType & matrixU() const
Returns the orthogonal matrix in the Schur decomposition.
Definition RealSchur.h:126
Eigen::RealSchur::getMaxIterations
Index getMaxIterations()
Returns the maximum number of iterations.
Definition RealSchur.h:211
Eigen::RealSchur::RealSchur
RealSchur(Index size=RowsAtCompileTime==Dynamic ? 1 :RowsAtCompileTime)
Default constructor.
Definition RealSchur.h:83
Eigen::RealSchur::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition RealSchur.h:193
Eigen::RealSchur::matrixT
const MatrixType & matrixT() const
Returns the quasi-triangular matrix in the Schur decomposition.
Definition RealSchur.h:143
Eigen::RealSchur::compute
RealSchur & compute(const MatrixType &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
Definition RealSchur.h:246
Eigen::RealSchur::setMaxIterations
RealSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition RealSchur.h:204
Eigen::RealSchur::m_maxIterationsPerRow
static const int m_maxIterationsPerRow
Maximum number of iterations per row.
Definition RealSchur.h:221
Eigen::ComputationInfo
ComputationInfo
Definition Constants.h:374
Eigen::NoConvergence
@ NoConvergence
Definition Constants.h:380
Eigen::Success
@ Success
Definition Constants.h:376