Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
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;
68
71
83 explicit 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 template<typename InputType>
104 explicit RealSchur(const EigenBase<InputType>& matrix, bool computeU = true)
105 : m_matT(matrix.rows(),matrix.cols()),
106 m_matU(matrix.rows(),matrix.cols()),
107 m_workspaceVector(matrix.rows()),
108 m_hess(matrix.rows()),
109 m_isInitialized(false),
110 m_matUisUptodate(false),
111 m_maxIters(-1)
112 {
113 compute(matrix.derived(), computeU);
114 }
115
127 const MatrixType& matrixU() const
128 {
129 eigen_assert(m_isInitialized && "RealSchur is not initialized.");
130 eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the RealSchur decomposition.");
131 return m_matU;
132 }
133
144 const MatrixType& matrixT() const
145 {
146 eigen_assert(m_isInitialized && "RealSchur is not initialized.");
147 return m_matT;
148 }
149
169 template<typename InputType>
170 RealSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true);
171
189 template<typename HessMatrixType, typename OrthMatrixType>
190 RealSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU);
196 {
197 eigen_assert(m_isInitialized && "RealSchur is not initialized.");
198 return m_info;
199 }
200
207 {
208 m_maxIters = maxIters;
209 return *this;
210 }
211
214 {
215 return m_maxIters;
216 }
217
223 static const int m_maxIterationsPerRow = 40;
224
225 private:
226
227 MatrixType m_matT;
228 MatrixType m_matU;
229 ColumnVectorType m_workspaceVector;
231 ComputationInfo m_info;
232 bool m_isInitialized;
233 bool m_matUisUptodate;
234 Index m_maxIters;
235
237
238 Scalar computeNormOfT();
239 Index findSmallSubdiagEntry(Index iu, const Scalar& considerAsZero);
240 void splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift);
241 void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo);
242 void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector);
243 void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace);
244};
245
246
247template<typename MatrixType>
248template<typename InputType>
250{
251 const Scalar considerAsZero = (std::numeric_limits<Scalar>::min)();
252
253 eigen_assert(matrix.cols() == matrix.rows());
254 Index maxIters = m_maxIters;
255 if (maxIters == -1)
256 maxIters = m_maxIterationsPerRow * matrix.rows();
257
258 Scalar scale = matrix.derived().cwiseAbs().maxCoeff();
259 if(scale<considerAsZero)
260 {
261 m_matT.setZero(matrix.rows(),matrix.cols());
262 if(computeU)
263 m_matU.setIdentity(matrix.rows(),matrix.cols());
264 m_info = Success;
265 m_isInitialized = true;
266 m_matUisUptodate = computeU;
267 return *this;
268 }
269
270 // Step 1. Reduce to Hessenberg form
271 m_hess.compute(matrix.derived()/scale);
272
273 // Step 2. Reduce to real Schur form
274 // Note: we copy m_hess.matrixQ() into m_matU here and not in computeFromHessenberg
275 // to be able to pass our working-space buffer for the Householder to Dense evaluation.
276 m_workspaceVector.resize(matrix.cols());
277 if(computeU)
278 m_hess.matrixQ().evalTo(m_matU, m_workspaceVector);
279 computeFromHessenberg(m_hess.matrixH(), m_matU, computeU);
280
281 m_matT *= scale;
282
283 return *this;
284}
285template<typename MatrixType>
286template<typename HessMatrixType, typename OrthMatrixType>
287RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
288{
289 using std::abs;
290
291 m_matT = matrixH;
292 m_workspaceVector.resize(m_matT.cols());
293 if(computeU && !internal::is_same_dense(m_matU,matrixQ))
294 m_matU = matrixQ;
295
296 Index maxIters = m_maxIters;
297 if (maxIters == -1)
298 maxIters = m_maxIterationsPerRow * matrixH.rows();
299 Scalar* workspace = &m_workspaceVector.coeffRef(0);
300
301 // The matrix m_matT is divided in three parts.
302 // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
303 // Rows il,...,iu is the part we are working on (the active window).
304 // Rows iu+1,...,end are already brought in triangular form.
305 Index iu = m_matT.cols() - 1;
306 Index iter = 0; // iteration count for current eigenvalue
307 Index totalIter = 0; // iteration count for whole matrix
308 Scalar exshift(0); // sum of exceptional shifts
309 Scalar norm = computeNormOfT();
310 // sub-diagonal entries smaller than considerAsZero will be treated as zero.
311 // We use eps^2 to enable more precision in small eigenvalues.
312 Scalar considerAsZero = numext::maxi<Scalar>( norm * numext::abs2(NumTraits<Scalar>::epsilon()),
313 (std::numeric_limits<Scalar>::min)() );
314
315 if(norm!=Scalar(0))
316 {
317 while (iu >= 0)
318 {
319 Index il = findSmallSubdiagEntry(iu,considerAsZero);
320
321 // Check for convergence
322 if (il == iu) // One root found
323 {
324 m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
325 if (iu > 0)
326 m_matT.coeffRef(iu, iu-1) = Scalar(0);
327 iu--;
328 iter = 0;
329 }
330 else if (il == iu-1) // Two roots found
331 {
332 splitOffTwoRows(iu, computeU, exshift);
333 iu -= 2;
334 iter = 0;
335 }
336 else // No convergence yet
337 {
338 // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG )
339 Vector3s firstHouseholderVector = Vector3s::Zero(), shiftInfo;
340 computeShift(iu, iter, exshift, shiftInfo);
341 iter = iter + 1;
342 totalIter = totalIter + 1;
343 if (totalIter > maxIters) break;
344 Index im;
345 initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
346 performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
347 }
348 }
349 }
350 if(totalIter <= maxIters)
351 m_info = Success;
352 else
353 m_info = NoConvergence;
354
355 m_isInitialized = true;
356 m_matUisUptodate = computeU;
357 return *this;
358}
359
361template<typename MatrixType>
362inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
363{
364 const Index size = m_matT.cols();
365 // FIXME to be efficient the following would requires a triangular reduxion code
366 // Scalar norm = m_matT.upper().cwiseAbs().sum()
367 // + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
368 Scalar norm(0);
369 for (Index j = 0; j < size; ++j)
370 norm += m_matT.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum();
371 return norm;
372}
373
375template<typename MatrixType>
376inline Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, const Scalar& considerAsZero)
377{
378 using std::abs;
379 Index res = iu;
380 while (res > 0)
381 {
382 Scalar s = abs(m_matT.coeff(res-1,res-1)) + abs(m_matT.coeff(res,res));
383
384 s = numext::maxi<Scalar>(s * NumTraits<Scalar>::epsilon(), considerAsZero);
385
386 if (abs(m_matT.coeff(res,res-1)) <= s)
387 break;
388 res--;
389 }
390 return res;
391}
392
394template<typename MatrixType>
395inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift)
396{
397 using std::sqrt;
398 using std::abs;
399 const Index size = m_matT.cols();
400
401 // The eigenvalues of the 2x2 matrix [a b; c d] are
402 // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc
403 Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
404 Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); // q = tr^2 / 4 - det = discr/4
405 m_matT.coeffRef(iu,iu) += exshift;
406 m_matT.coeffRef(iu-1,iu-1) += exshift;
407
408 if (q >= Scalar(0)) // Two real eigenvalues
409 {
410 Scalar z = sqrt(abs(q));
411 JacobiRotation<Scalar> rot;
412 if (p >= Scalar(0))
413 rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
414 else
415 rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
416
417 m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
418 m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot);
419 m_matT.coeffRef(iu, iu-1) = Scalar(0);
420 if (computeU)
421 m_matU.applyOnTheRight(iu-1, iu, rot);
422 }
423
424 if (iu > 1)
425 m_matT.coeffRef(iu-1, iu-2) = Scalar(0);
426}
427
429template<typename MatrixType>
430inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo)
431{
432 using std::sqrt;
433 using std::abs;
434 shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
435 shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
436 shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
437
438 // Wilkinson's original ad hoc shift
439 if (iter == 10)
440 {
441 exshift += shiftInfo.coeff(0);
442 for (Index i = 0; i <= iu; ++i)
443 m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
444 Scalar s = abs(m_matT.coeff(iu,iu-1)) + abs(m_matT.coeff(iu-1,iu-2));
445 shiftInfo.coeffRef(0) = Scalar(0.75) * s;
446 shiftInfo.coeffRef(1) = Scalar(0.75) * s;
447 shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
448 }
449
450 // MATLAB's new ad hoc shift
451 if (iter == 30)
452 {
453 Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
454 s = s * s + shiftInfo.coeff(2);
455 if (s > Scalar(0))
456 {
457 s = sqrt(s);
458 if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
459 s = -s;
460 s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
461 s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
462 exshift += s;
463 for (Index i = 0; i <= iu; ++i)
464 m_matT.coeffRef(i,i) -= s;
465 shiftInfo.setConstant(Scalar(0.964));
466 }
467 }
468}
469
471template<typename MatrixType>
472inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector)
473{
474 using std::abs;
475 Vector3s& v = firstHouseholderVector; // alias to save typing
476
477 for (im = iu-2; im >= il; --im)
478 {
479 const Scalar Tmm = m_matT.coeff(im,im);
480 const Scalar r = shiftInfo.coeff(0) - Tmm;
481 const Scalar s = shiftInfo.coeff(1) - Tmm;
482 v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
483 v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s;
484 v.coeffRef(2) = m_matT.coeff(im+2,im+1);
485 if (im == il) {
486 break;
487 }
488 const Scalar lhs = m_matT.coeff(im,im-1) * (abs(v.coeff(1)) + abs(v.coeff(2)));
489 const Scalar rhs = v.coeff(0) * (abs(m_matT.coeff(im-1,im-1)) + abs(Tmm) + abs(m_matT.coeff(im+1,im+1)));
490 if (abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
491 break;
492 }
493}
494
496template<typename MatrixType>
497inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace)
498{
499 eigen_assert(im >= il);
500 eigen_assert(im <= iu-2);
501
502 const Index size = m_matT.cols();
503
504 for (Index k = im; k <= iu-2; ++k)
505 {
506 bool firstIteration = (k == im);
507
508 Vector3s v;
509 if (firstIteration)
510 v = firstHouseholderVector;
511 else
512 v = m_matT.template block<3,1>(k,k-1);
513
514 Scalar tau, beta;
515 Matrix<Scalar, 2, 1> ess;
516 v.makeHouseholder(ess, tau, beta);
517
518 if (beta != Scalar(0)) // if v is not zero
519 {
520 if (firstIteration && k > il)
521 m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
522 else if (!firstIteration)
523 m_matT.coeffRef(k,k-1) = beta;
524
525 // These Householder transformations form the O(n^3) part of the algorithm
526 m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
527 m_matT.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
528 if (computeU)
529 m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
530 }
531 }
532
533 Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2);
534 Scalar tau, beta;
535 Matrix<Scalar, 1, 1> ess;
536 v.makeHouseholder(ess, tau, beta);
537
538 if (beta != Scalar(0)) // if v is not zero
539 {
540 m_matT.coeffRef(iu-1, iu-2) = beta;
541 m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
542 m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
543 if (computeU)
544 m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
545 }
546
547 // clean up pollution due to round-off errors
548 for (Index i = im+2; i <= iu; ++i)
549 {
550 m_matT.coeffRef(i,i-2) = Scalar(0);
551 if (i > im+2)
552 m_matT.coeffRef(i,i-3) = Scalar(0);
553 }
554}
555
556} // end namespace Eigen
557
558#endif // EIGEN_REAL_SCHUR_H
Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation.
Definition: HessenbergDecomposition.h:58
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Performs a real Schur decomposition of a square matrix.
Definition: RealSchur.h:55
RealSchur & compute(const EigenBase< InputType > &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealSchur.h:195
const MatrixType & matrixU() const
Returns the orthogonal matrix in the Schur decomposition.
Definition: RealSchur.h:127
static const int m_maxIterationsPerRow
Maximum number of iterations per row.
Definition: RealSchur.h:223
RealSchur(Index size=RowsAtCompileTime==Dynamic ? 1 :RowsAtCompileTime)
Default constructor.
Definition: RealSchur.h:83
Eigen::Index Index
Definition: RealSchur.h:67
const MatrixType & matrixT() const
Returns the quasi-triangular matrix in the Schur decomposition.
Definition: RealSchur.h:144
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: RealSchur.h:213
RealSchur & computeFromHessenberg(const HessMatrixType &matrixH, const OrthMatrixType &matrixQ, bool computeU)
Computes Schur decomposition of a Hessenberg matrix H = Z T Z^T.
RealSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: RealSchur.h:206
RealSchur(const EigenBase< InputType > &matrix, bool computeU=true)
Constructor; computes real Schur decomposition of given matrix.
Definition: RealSchur.h:104
ComputationInfo
Definition: Constants.h:440
@ Success
Definition: Constants.h:442
@ NoConvergence
Definition: Constants.h:446
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const int Dynamic
Definition: Constants.h:22
Definition: EigenBase.h:30
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:63
Derived & derived()
Definition: EigenBase.h:46
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60