RealSchur.h
Go to the documentation of this file.
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 
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 typename MatrixType::Index Index;
68 
71 
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);
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 
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;
232  Index m_maxIters;
233 
235 
236  Scalar computeNormOfT();
237  Index findSmallSubdiagEntry(Index iu, const Scalar& norm);
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 
245 template<typename MatrixType>
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
258 
259  return *this;
260 }
261 template<typename MatrixType>
262 template<typename HessMatrixType, typename OrthMatrixType>
263 RealSchur<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();
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, norm);
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
324 
325  m_isInitialized = true;
326  m_matUisUptodate = computeU;
327  return *this;
328 }
329 
331 template<typename MatrixType>
332 inline 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 
345 template<typename MatrixType>
346 inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, const Scalar& norm)
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 (s == 0.0)
354  s = norm;
355  if (abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
356  break;
357  res--;
358  }
359  return res;
360 }
361 
363 template<typename MatrixType>
364 inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift)
365 {
366  using std::sqrt;
367  using std::abs;
368  const Index size = m_matT.cols();
369 
370  // The eigenvalues of the 2x2 matrix [a b; c d] are
371  // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc
372  Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
373  Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); // q = tr^2 / 4 - det = discr/4
374  m_matT.coeffRef(iu,iu) += exshift;
375  m_matT.coeffRef(iu-1,iu-1) += exshift;
376 
377  if (q >= Scalar(0)) // Two real eigenvalues
378  {
379  Scalar z = sqrt(abs(q));
381  if (p >= Scalar(0))
382  rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
383  else
384  rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
385 
386  m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
387  m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot);
388  m_matT.coeffRef(iu, iu-1) = Scalar(0);
389  if (computeU)
390  m_matU.applyOnTheRight(iu-1, iu, rot);
391  }
392 
393  if (iu > 1)
394  m_matT.coeffRef(iu-1, iu-2) = Scalar(0);
395 }
396 
398 template<typename MatrixType>
399 inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo)
400 {
401  using std::sqrt;
402  using std::abs;
403  shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
404  shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
405  shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
406 
407  // Wilkinson's original ad hoc shift
408  if (iter == 10)
409  {
410  exshift += shiftInfo.coeff(0);
411  for (Index i = 0; i <= iu; ++i)
412  m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
413  Scalar s = abs(m_matT.coeff(iu,iu-1)) + abs(m_matT.coeff(iu-1,iu-2));
414  shiftInfo.coeffRef(0) = Scalar(0.75) * s;
415  shiftInfo.coeffRef(1) = Scalar(0.75) * s;
416  shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
417  }
418 
419  // MATLAB's new ad hoc shift
420  if (iter == 30)
421  {
422  Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
423  s = s * s + shiftInfo.coeff(2);
424  if (s > Scalar(0))
425  {
426  s = sqrt(s);
427  if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
428  s = -s;
429  s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
430  s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
431  exshift += s;
432  for (Index i = 0; i <= iu; ++i)
433  m_matT.coeffRef(i,i) -= s;
434  shiftInfo.setConstant(Scalar(0.964));
435  }
436  }
437 }
438 
440 template<typename MatrixType>
441 inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector)
442 {
443  using std::abs;
444  Vector3s& v = firstHouseholderVector; // alias to save typing
445 
446  for (im = iu-2; im >= il; --im)
447  {
448  const Scalar Tmm = m_matT.coeff(im,im);
449  const Scalar r = shiftInfo.coeff(0) - Tmm;
450  const Scalar s = shiftInfo.coeff(1) - Tmm;
451  v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
452  v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s;
453  v.coeffRef(2) = m_matT.coeff(im+2,im+1);
454  if (im == il) {
455  break;
456  }
457  const Scalar lhs = m_matT.coeff(im,im-1) * (abs(v.coeff(1)) + abs(v.coeff(2)));
458  const Scalar rhs = v.coeff(0) * (abs(m_matT.coeff(im-1,im-1)) + abs(Tmm) + abs(m_matT.coeff(im+1,im+1)));
459  if (abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
460  {
461  break;
462  }
463  }
464 }
465 
467 template<typename MatrixType>
468 inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace)
469 {
470  eigen_assert(im >= il);
471  eigen_assert(im <= iu-2);
472 
473  const Index size = m_matT.cols();
474 
475  for (Index k = im; k <= iu-2; ++k)
476  {
477  bool firstIteration = (k == im);
478 
479  Vector3s v;
480  if (firstIteration)
481  v = firstHouseholderVector;
482  else
483  v = m_matT.template block<3,1>(k,k-1);
484 
485  Scalar tau, beta;
487  v.makeHouseholder(ess, tau, beta);
488 
489  if (beta != Scalar(0)) // if v is not zero
490  {
491  if (firstIteration && k > il)
492  m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
493  else if (!firstIteration)
494  m_matT.coeffRef(k,k-1) = beta;
495 
496  // These Householder transformations form the O(n^3) part of the algorithm
497  m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
498  m_matT.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
499  if (computeU)
500  m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
501  }
502  }
503 
504  Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2);
505  Scalar tau, beta;
507  v.makeHouseholder(ess, tau, beta);
508 
509  if (beta != Scalar(0)) // if v is not zero
510  {
511  m_matT.coeffRef(iu-1, iu-2) = beta;
512  m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
513  m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
514  if (computeU)
515  m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
516  }
517 
518  // clean up pollution due to round-off errors
519  for (Index i = im+2; i <= iu; ++i)
520  {
521  m_matT.coeffRef(i,i-2) = Scalar(0);
522  if (i > im+2)
523  m_matT.coeffRef(i,i-3) = Scalar(0);
524  }
525 }
526 
527 } // end namespace Eigen
528 
529 #endif // EIGEN_REAL_SCHUR_H
MatrixType::Index Index
Definition: RealSchur.h:67
bool m_matUisUptodate
Definition: RealSchur.h:231
HouseholderSequenceType matrixQ() const
Reconstructs the orthogonal matrix Q in the decomposition.
IntermediateState sqrt(const Expression &arg)
MatrixHReturnType matrixH() const
Constructs the Hessenberg matrix H in the decomposition.
ComputationInfo m_info
Definition: RealSchur.h:229
Performs a real Schur decomposition of a square matrix.
Definition: RealSchur.h:54
RealSchur(const MatrixType &matrix, bool computeU=true)
Constructor; computes real Schur decomposition of given matrix.
Definition: RealSchur.h:103
void makeGivens(const Scalar &p, const Scalar &q, Scalar *z=0)
Definition: Jacobi.h:148
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:126
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
Rotation given by a cosine-sine pair.
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealSchur.h:193
const MatrixType & matrixT() const
Returns the quasi-triangular matrix in the Schur decomposition.
Definition: RealSchur.h:143
MatrixType m_matT
Definition: RealSchur.h:225
Index findSmallSubdiagEntry(Index iu, const Scalar &norm)
Definition: RealSchur.h:346
Scalar computeNormOfT()
Definition: RealSchur.h:332
void computeShift(Index iu, Index iter, Scalar &exshift, Vector3s &shiftInfo)
Definition: RealSchur.h:399
MatrixType m_matU
Definition: RealSchur.h:226
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: RealSchur.h:211
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
HessenbergDecomposition< MatrixType > m_hess
Definition: RealSchur.h:228
bool m_isInitialized
Definition: RealSchur.h:230
EIGEN_STRONG_INLINE const Scalar & coeff(Index rowId, Index colId) const
_MatrixType MatrixType
Definition: RealSchur.h:57
HessenbergDecomposition & compute(const MatrixType &matrix)
Computes Hessenberg decomposition of given matrix.
MatrixType::Scalar Scalar
Definition: RealSchur.h:65
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > EigenvalueType
Definition: RealSchur.h:69
EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
EIGEN_STRONG_INLINE void resize(Index nbRows, Index nbCols)
static const int m_maxIterationsPerRow
Maximum number of iterations per row.
Definition: RealSchur.h:221
Provides a generic way to set and pass user-specified options.
Definition: options.hpp:65
JacobiRotation adjoint() const
Definition: Jacobi.h:62
void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s &firstHouseholderVector, Scalar *workspace)
Definition: RealSchur.h:468
Derived & setConstant(Index size, const Scalar &value)
RealSchur & compute(const MatrixType &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
Definition: RealSchur.h:246
#define v
void rhs(const real_t *x, real_t *f)
Matrix< Scalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > ColumnVectorType
Definition: RealSchur.h:70
Matrix< Scalar, 3, 1 > Vector3s
Definition: RealSchur.h:234
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:204
void initFrancisQRStep(Index il, Index iu, const Vector3s &shiftInfo, Index &im, Vector3s &firstHouseholderVector)
Definition: RealSchur.h:441
void splitOffTwoRows(Index iu, bool computeU, const Scalar &exshift)
Definition: RealSchur.h:364
std::complex< typename NumTraits< Scalar >::Real > ComplexScalar
Definition: RealSchur.h:66
ColumnVectorType m_workspaceVector
Definition: RealSchur.h:227
#define eigen_assert(x)


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:35:03