Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.3.9
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 
16 namespace Eigen {
17 
54 template<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 Eigen::Index Index;
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 
247 template<typename MatrixType>
248 template<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  computeFromHessenberg(m_hess.matrixH(), m_hess.matrixQ(), computeU);
275 
276  m_matT *= scale;
277 
278  return *this;
279 }
280 template<typename MatrixType>
281 template<typename HessMatrixType, typename OrthMatrixType>
282 RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
283 {
284  using std::abs;
285 
286  m_matT = matrixH;
287  if(computeU)
288  m_matU = matrixQ;
289 
290  Index maxIters = m_maxIters;
291  if (maxIters == -1)
292  maxIters = m_maxIterationsPerRow * matrixH.rows();
293  m_workspaceVector.resize(m_matT.cols());
294  Scalar* workspace = &m_workspaceVector.coeffRef(0);
295 
296  // The matrix m_matT is divided in three parts.
297  // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
298  // Rows il,...,iu is the part we are working on (the active window).
299  // Rows iu+1,...,end are already brought in triangular form.
300  Index iu = m_matT.cols() - 1;
301  Index iter = 0; // iteration count for current eigenvalue
302  Index totalIter = 0; // iteration count for whole matrix
303  Scalar exshift(0); // sum of exceptional shifts
304  Scalar norm = computeNormOfT();
305  // sub-diagonal entries smaller than considerAsZero will be treated as zero.
306  // We use eps^2 to enable more precision in small eigenvalues.
307  Scalar considerAsZero = numext::maxi<Scalar>( norm * numext::abs2(NumTraits<Scalar>::epsilon()),
308  (std::numeric_limits<Scalar>::min)() );
309 
310  if(norm!=Scalar(0))
311  {
312  while (iu >= 0)
313  {
314  Index il = findSmallSubdiagEntry(iu,considerAsZero);
315 
316  // Check for convergence
317  if (il == iu) // One root found
318  {
319  m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
320  if (iu > 0)
321  m_matT.coeffRef(iu, iu-1) = Scalar(0);
322  iu--;
323  iter = 0;
324  }
325  else if (il == iu-1) // Two roots found
326  {
327  splitOffTwoRows(iu, computeU, exshift);
328  iu -= 2;
329  iter = 0;
330  }
331  else // No convergence yet
332  {
333  // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG )
334  Vector3s firstHouseholderVector = Vector3s::Zero(), shiftInfo;
335  computeShift(iu, iter, exshift, shiftInfo);
336  iter = iter + 1;
337  totalIter = totalIter + 1;
338  if (totalIter > maxIters) break;
339  Index im;
340  initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
341  performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
342  }
343  }
344  }
345  if(totalIter <= maxIters)
346  m_info = Success;
347  else
348  m_info = NoConvergence;
349 
350  m_isInitialized = true;
351  m_matUisUptodate = computeU;
352  return *this;
353 }
354 
356 template<typename MatrixType>
357 inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
358 {
359  const Index size = m_matT.cols();
360  // FIXME to be efficient the following would requires a triangular reduxion code
361  // Scalar norm = m_matT.upper().cwiseAbs().sum()
362  // + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
363  Scalar norm(0);
364  for (Index j = 0; j < size; ++j)
365  norm += m_matT.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum();
366  return norm;
367 }
368 
370 template<typename MatrixType>
371 inline Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, const Scalar& considerAsZero)
372 {
373  using std::abs;
374  Index res = iu;
375  while (res > 0)
376  {
377  Scalar s = abs(m_matT.coeff(res-1,res-1)) + abs(m_matT.coeff(res,res));
378 
379  s = numext::maxi<Scalar>(s * NumTraits<Scalar>::epsilon(), considerAsZero);
380 
381  if (abs(m_matT.coeff(res,res-1)) <= s)
382  break;
383  res--;
384  }
385  return res;
386 }
387 
389 template<typename MatrixType>
390 inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift)
391 {
392  using std::sqrt;
393  using std::abs;
394  const Index size = m_matT.cols();
395 
396  // The eigenvalues of the 2x2 matrix [a b; c d] are
397  // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc
398  Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
399  Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); // q = tr^2 / 4 - det = discr/4
400  m_matT.coeffRef(iu,iu) += exshift;
401  m_matT.coeffRef(iu-1,iu-1) += exshift;
402 
403  if (q >= Scalar(0)) // Two real eigenvalues
404  {
405  Scalar z = sqrt(abs(q));
406  JacobiRotation<Scalar> rot;
407  if (p >= Scalar(0))
408  rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
409  else
410  rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
411 
412  m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
413  m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot);
414  m_matT.coeffRef(iu, iu-1) = Scalar(0);
415  if (computeU)
416  m_matU.applyOnTheRight(iu-1, iu, rot);
417  }
418 
419  if (iu > 1)
420  m_matT.coeffRef(iu-1, iu-2) = Scalar(0);
421 }
422 
424 template<typename MatrixType>
425 inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo)
426 {
427  using std::sqrt;
428  using std::abs;
429  shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
430  shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
431  shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
432 
433  // Wilkinson's original ad hoc shift
434  if (iter == 10)
435  {
436  exshift += shiftInfo.coeff(0);
437  for (Index i = 0; i <= iu; ++i)
438  m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
439  Scalar s = abs(m_matT.coeff(iu,iu-1)) + abs(m_matT.coeff(iu-1,iu-2));
440  shiftInfo.coeffRef(0) = Scalar(0.75) * s;
441  shiftInfo.coeffRef(1) = Scalar(0.75) * s;
442  shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
443  }
444 
445  // MATLAB's new ad hoc shift
446  if (iter == 30)
447  {
448  Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
449  s = s * s + shiftInfo.coeff(2);
450  if (s > Scalar(0))
451  {
452  s = sqrt(s);
453  if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
454  s = -s;
455  s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
456  s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
457  exshift += s;
458  for (Index i = 0; i <= iu; ++i)
459  m_matT.coeffRef(i,i) -= s;
460  shiftInfo.setConstant(Scalar(0.964));
461  }
462  }
463 }
464 
466 template<typename MatrixType>
467 inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector)
468 {
469  using std::abs;
470  Vector3s& v = firstHouseholderVector; // alias to save typing
471 
472  for (im = iu-2; im >= il; --im)
473  {
474  const Scalar Tmm = m_matT.coeff(im,im);
475  const Scalar r = shiftInfo.coeff(0) - Tmm;
476  const Scalar s = shiftInfo.coeff(1) - Tmm;
477  v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
478  v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s;
479  v.coeffRef(2) = m_matT.coeff(im+2,im+1);
480  if (im == il) {
481  break;
482  }
483  const Scalar lhs = m_matT.coeff(im,im-1) * (abs(v.coeff(1)) + abs(v.coeff(2)));
484  const Scalar rhs = v.coeff(0) * (abs(m_matT.coeff(im-1,im-1)) + abs(Tmm) + abs(m_matT.coeff(im+1,im+1)));
485  if (abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
486  break;
487  }
488 }
489 
491 template<typename MatrixType>
492 inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace)
493 {
494  eigen_assert(im >= il);
495  eigen_assert(im <= iu-2);
496 
497  const Index size = m_matT.cols();
498 
499  for (Index k = im; k <= iu-2; ++k)
500  {
501  bool firstIteration = (k == im);
502 
503  Vector3s v;
504  if (firstIteration)
505  v = firstHouseholderVector;
506  else
507  v = m_matT.template block<3,1>(k,k-1);
508 
509  Scalar tau, beta;
510  Matrix<Scalar, 2, 1> ess;
511  v.makeHouseholder(ess, tau, beta);
512 
513  if (beta != Scalar(0)) // if v is not zero
514  {
515  if (firstIteration && k > il)
516  m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
517  else if (!firstIteration)
518  m_matT.coeffRef(k,k-1) = beta;
519 
520  // These Householder transformations form the O(n^3) part of the algorithm
521  m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
522  m_matT.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
523  if (computeU)
524  m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
525  }
526  }
527 
528  Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2);
529  Scalar tau, beta;
530  Matrix<Scalar, 1, 1> ess;
531  v.makeHouseholder(ess, tau, beta);
532 
533  if (beta != Scalar(0)) // if v is not zero
534  {
535  m_matT.coeffRef(iu-1, iu-2) = beta;
536  m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
537  m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
538  if (computeU)
539  m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
540  }
541 
542  // clean up pollution due to round-off errors
543  for (Index i = im+2; i <= iu; ++i)
544  {
545  m_matT.coeffRef(i,i-2) = Scalar(0);
546  if (i > im+2)
547  m_matT.coeffRef(i,i-3) = Scalar(0);
548  }
549 }
550 
551 } // end namespace Eigen
552 
553 #endif // EIGEN_REAL_SCHUR_H
Eigen::RealSchur::Index
Eigen::Index Index
Definition: RealSchur.h:67
Eigen::EigenBase::cols
Index cols() const
Definition: EigenBase.h:62
Eigen
Namespace containing all symbols from the Eigen library.
Definition: Core:309
Eigen::EigenBase
Definition: EigenBase.h:30
Eigen::RealSchur::matrixT
const MatrixType & matrixT() const
Returns the quasi-triangular matrix in the Schur decomposition.
Definition: RealSchur.h:144
Eigen::RealSchur::RealSchur
RealSchur(const EigenBase< InputType > &matrix, bool computeU=true)
Constructor; computes real Schur decomposition of given matrix.
Definition: RealSchur.h:104
Eigen::Success
@ Success
Definition: Constants.h:432
Eigen::RealSchur::compute
RealSchur & compute(const EigenBase< InputType > &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
Eigen::RealSchur::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealSchur.h:195
Eigen::RealSchur
Performs a real Schur decomposition of a square matrix.
Definition: RealSchur.h:55
Eigen::NoConvergence
@ NoConvergence
Definition: Constants.h:436
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::m_maxIterationsPerRow
static const int m_maxIterationsPerRow
Maximum number of iterations per row.
Definition: RealSchur.h:223
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:21
Eigen::RealSchur::getMaxIterations
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: RealSchur.h:213
Eigen::EigenBase::derived
Derived & derived()
Definition: EigenBase.h:45
Eigen::RealSchur::RealSchur
RealSchur(Index size=RowsAtCompileTime==Dynamic ? 1 :RowsAtCompileTime)
Default constructor.
Definition: RealSchur.h:83
Eigen::Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 >
Eigen::EigenBase::rows
Index rows() const
Definition: EigenBase.h:59
Eigen::RealSchur::matrixU
const MatrixType & matrixU() const
Returns the orthogonal matrix in the Schur decomposition.
Definition: RealSchur.h:127
Eigen::RealSchur::setMaxIterations
RealSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: RealSchur.h:206
Eigen::ComputationInfo
ComputationInfo
Definition: Constants.h:430
Eigen::HessenbergDecomposition< MatrixType >
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33