Eigen  3.2.92
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 
206  RealSchur& setMaxIterations(Index maxIters)
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);
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  eigen_assert(matrix.cols() == matrix.rows());
252  Index maxIters = m_maxIters;
253  if (maxIters == -1)
254  maxIters = m_maxIterationsPerRow * matrix.rows();
255 
256  // Step 1. Reduce to Hessenberg form
257  m_hess.compute(matrix.derived());
258 
259  // Step 2. Reduce to real Schur form
260  computeFromHessenberg(m_hess.matrixH(), m_hess.matrixQ(), computeU);
261 
262  return *this;
263 }
264 template<typename MatrixType>
265 template<typename HessMatrixType, typename OrthMatrixType>
266 RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
267 {
268  m_matT = matrixH;
269  if(computeU)
270  m_matU = matrixQ;
271 
272  Index maxIters = m_maxIters;
273  if (maxIters == -1)
274  maxIters = m_maxIterationsPerRow * matrixH.rows();
275  m_workspaceVector.resize(m_matT.cols());
276  Scalar* workspace = &m_workspaceVector.coeffRef(0);
277 
278  // The matrix m_matT is divided in three parts.
279  // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
280  // Rows il,...,iu is the part we are working on (the active window).
281  // Rows iu+1,...,end are already brought in triangular form.
282  Index iu = m_matT.cols() - 1;
283  Index iter = 0; // iteration count for current eigenvalue
284  Index totalIter = 0; // iteration count for whole matrix
285  Scalar exshift(0); // sum of exceptional shifts
286  Scalar norm = computeNormOfT();
287 
288  if(norm!=0)
289  {
290  while (iu >= 0)
291  {
292  Index il = findSmallSubdiagEntry(iu);
293 
294  // Check for convergence
295  if (il == iu) // One root found
296  {
297  m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
298  if (iu > 0)
299  m_matT.coeffRef(iu, iu-1) = Scalar(0);
300  iu--;
301  iter = 0;
302  }
303  else if (il == iu-1) // Two roots found
304  {
305  splitOffTwoRows(iu, computeU, exshift);
306  iu -= 2;
307  iter = 0;
308  }
309  else // No convergence yet
310  {
311  // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG )
312  Vector3s firstHouseholderVector(0,0,0), shiftInfo;
313  computeShift(iu, iter, exshift, shiftInfo);
314  iter = iter + 1;
315  totalIter = totalIter + 1;
316  if (totalIter > maxIters) break;
317  Index im;
318  initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
319  performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
320  }
321  }
322  }
323  if(totalIter <= maxIters)
324  m_info = Success;
325  else
326  m_info = NoConvergence;
327 
328  m_isInitialized = true;
329  m_matUisUptodate = computeU;
330  return *this;
331 }
332 
334 template<typename MatrixType>
335 inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
336 {
337  const Index size = m_matT.cols();
338  // FIXME to be efficient the following would requires a triangular reduxion code
339  // Scalar norm = m_matT.upper().cwiseAbs().sum()
340  // + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
341  Scalar norm(0);
342  for (Index j = 0; j < size; ++j)
343  norm += m_matT.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum();
344  return norm;
345 }
346 
348 template<typename MatrixType>
349 inline Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu)
350 {
351  using std::abs;
352  Index res = iu;
353  while (res > 0)
354  {
355  Scalar s = abs(m_matT.coeff(res-1,res-1)) + abs(m_matT.coeff(res,res));
356  if (abs(m_matT.coeff(res,res-1)) <= NumTraits<Scalar>::epsilon() * s)
357  break;
358  res--;
359  }
360  return res;
361 }
362 
364 template<typename MatrixType>
365 inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift)
366 {
367  using std::sqrt;
368  using std::abs;
369  const Index size = m_matT.cols();
370 
371  // The eigenvalues of the 2x2 matrix [a b; c d] are
372  // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc
373  Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
374  Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); // q = tr^2 / 4 - det = discr/4
375  m_matT.coeffRef(iu,iu) += exshift;
376  m_matT.coeffRef(iu-1,iu-1) += exshift;
377 
378  if (q >= Scalar(0)) // Two real eigenvalues
379  {
380  Scalar z = sqrt(abs(q));
381  JacobiRotation<Scalar> rot;
382  if (p >= Scalar(0))
383  rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
384  else
385  rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
386 
387  m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
388  m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot);
389  m_matT.coeffRef(iu, iu-1) = Scalar(0);
390  if (computeU)
391  m_matU.applyOnTheRight(iu-1, iu, rot);
392  }
393 
394  if (iu > 1)
395  m_matT.coeffRef(iu-1, iu-2) = Scalar(0);
396 }
397 
399 template<typename MatrixType>
400 inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo)
401 {
402  using std::sqrt;
403  using std::abs;
404  shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
405  shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
406  shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
407 
408  // Wilkinson's original ad hoc shift
409  if (iter == 10)
410  {
411  exshift += shiftInfo.coeff(0);
412  for (Index i = 0; i <= iu; ++i)
413  m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
414  Scalar s = abs(m_matT.coeff(iu,iu-1)) + abs(m_matT.coeff(iu-1,iu-2));
415  shiftInfo.coeffRef(0) = Scalar(0.75) * s;
416  shiftInfo.coeffRef(1) = Scalar(0.75) * s;
417  shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
418  }
419 
420  // MATLAB's new ad hoc shift
421  if (iter == 30)
422  {
423  Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
424  s = s * s + shiftInfo.coeff(2);
425  if (s > Scalar(0))
426  {
427  s = sqrt(s);
428  if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
429  s = -s;
430  s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
431  s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
432  exshift += s;
433  for (Index i = 0; i <= iu; ++i)
434  m_matT.coeffRef(i,i) -= s;
435  shiftInfo.setConstant(Scalar(0.964));
436  }
437  }
438 }
439 
441 template<typename MatrixType>
442 inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector)
443 {
444  using std::abs;
445  Vector3s& v = firstHouseholderVector; // alias to save typing
446 
447  for (im = iu-2; im >= il; --im)
448  {
449  const Scalar Tmm = m_matT.coeff(im,im);
450  const Scalar r = shiftInfo.coeff(0) - Tmm;
451  const Scalar s = shiftInfo.coeff(1) - Tmm;
452  v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
453  v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s;
454  v.coeffRef(2) = m_matT.coeff(im+2,im+1);
455  if (im == il) {
456  break;
457  }
458  const Scalar lhs = m_matT.coeff(im,im-1) * (abs(v.coeff(1)) + abs(v.coeff(2)));
459  const Scalar rhs = v.coeff(0) * (abs(m_matT.coeff(im-1,im-1)) + abs(Tmm) + abs(m_matT.coeff(im+1,im+1)));
460  if (abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
461  break;
462  }
463 }
464 
466 template<typename MatrixType>
467 inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace)
468 {
469  eigen_assert(im >= il);
470  eigen_assert(im <= iu-2);
471 
472  const Index size = m_matT.cols();
473 
474  for (Index k = im; k <= iu-2; ++k)
475  {
476  bool firstIteration = (k == im);
477 
478  Vector3s v;
479  if (firstIteration)
480  v = firstHouseholderVector;
481  else
482  v = m_matT.template block<3,1>(k,k-1);
483 
484  Scalar tau, beta;
485  Matrix<Scalar, 2, 1> ess;
486  v.makeHouseholder(ess, tau, beta);
487 
488  if (beta != Scalar(0)) // if v is not zero
489  {
490  if (firstIteration && k > il)
491  m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
492  else if (!firstIteration)
493  m_matT.coeffRef(k,k-1) = beta;
494 
495  // These Householder transformations form the O(n^3) part of the algorithm
496  m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
497  m_matT.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
498  if (computeU)
499  m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
500  }
501  }
502 
503  Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2);
504  Scalar tau, beta;
505  Matrix<Scalar, 1, 1> ess;
506  v.makeHouseholder(ess, tau, beta);
507 
508  if (beta != Scalar(0)) // if v is not zero
509  {
510  m_matT.coeffRef(iu-1, iu-2) = beta;
511  m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
512  m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
513  if (computeU)
514  m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
515  }
516 
517  // clean up pollution due to round-off errors
518  for (Index i = im+2; i <= iu; ++i)
519  {
520  m_matT.coeffRef(i,i-2) = Scalar(0);
521  if (i > im+2)
522  m_matT.coeffRef(i,i-3) = Scalar(0);
523  }
524 }
525 
526 } // end namespace Eigen
527 
528 #endif // EIGEN_REAL_SCHUR_H
Performs a real Schur decomposition of a square matrix.
Definition: RealSchur.h:54
RealSchur & computeFromHessenberg(const HessMatrixType &matrixH, const OrthMatrixType &matrixQ, bool computeU)
Computes Schur decomposition of a Hessenberg matrix H = Z T Z^T.
const MatrixType & matrixU() const
Returns the orthogonal matrix in the Schur decomposition.
Definition: RealSchur.h:127
Definition: LDLT.h:16
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealSchur.h:195
Derived & derived()
Definition: EigenBase.h:44
Index rows() const
Definition: EigenBase.h:58
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
Definition: EigenBase.h:28
Eigen::Index Index
Definition: RealSchur.h:67
static const int m_maxIterationsPerRow
Maximum number of iterations per row.
Definition: RealSchur.h:223
RealSchur & compute(const EigenBase< InputType > &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
Definition: Constants.h:432
Index cols() const
Definition: EigenBase.h:61
RealSchur(Index size=RowsAtCompileTime==Dynamic?1:RowsAtCompileTime)
Default constructor.
Definition: RealSchur.h:83
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:430
Definition: Constants.h:436