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