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;
67  typedef Eigen::Index Index;
68 
71 
84  : m_matT(size, size),
85  m_matU(size, 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;
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  // 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)
280 
281  m_matT *= scale;
282 
283  return *this;
284 }
285 template<typename MatrixType>
286 template<typename HessMatrixType, typename OrthMatrixType>
287 RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
288 {
289  using std::abs;
290 
291  m_matT = matrixH;
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()),
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
354 
355  m_isInitialized = true;
356  m_matUisUptodate = computeU;
357  return *this;
358 }
359 
361 template<typename MatrixType>
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 
375 template<typename MatrixType>
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 
394 template<typename MatrixType>
395 inline 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));
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 
429 template<typename MatrixType>
430 inline 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 
471 template<typename MatrixType>
472 inline 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 
496 template<typename MatrixType>
497 inline 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;
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;
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
bool m_matUisUptodate
Definition: RealSchur.h:233
int EIGEN_BLAS_FUNC() rot(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, RealScalar *ps)
internal::traits< Derived >::Scalar Scalar
ComputationInfo m_info
Definition: RealSchur.h:231
Performs a real Schur decomposition of a square matrix.
Definition: RealSchur.h:54
#define min(a, b)
Definition: datatypes.h:19
RealSchur & computeFromHessenberg(const HessMatrixType &matrixH, const OrthMatrixType &matrixQ, bool computeU)
Computes Schur decomposition of a Hessenberg matrix H = Z T Z^T.
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Rotation given by a cosine-sine pair.
iterator iter(handle obj)
Definition: pytypes.h:2273
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:63
RealSchur(Index size=RowsAtCompileTime==Dynamic ? 1 :RowsAtCompileTime)
Default constructor.
Definition: RealSchur.h:83
MatrixType m_matT
Definition: RealSchur.h:227
Scalar computeNormOfT()
Definition: RealSchur.h:362
void computeShift(Index iu, Index iter, Scalar &exshift, Vector3s &shiftInfo)
Definition: RealSchur.h:430
static double epsilon
Definition: testRot3.cpp:37
MatrixType m_matU
Definition: RealSchur.h:228
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: RealSchur.h:213
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
HessenbergDecomposition< MatrixType > m_hess
Definition: RealSchur.h:230
Eigen::Index Index
Definition: RealSchur.h:67
bool m_isInitialized
Definition: RealSchur.h:232
HessenbergDecomposition & compute(const EigenBase< InputType > &matrix)
Computes Hessenberg decomposition of given matrix.
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
MatrixHReturnType matrixH() const
Constructs the Hessenberg matrix H in the decomposition.
EIGEN_DEVICE_FUNC void makeGivens(const Scalar &p, const Scalar &q, Scalar *r=0)
Definition: Jacobi.h:162
_MatrixType MatrixType
Definition: RealSchur.h:57
MatrixType::Scalar Scalar
Definition: RealSchur.h:65
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar & coeff(Index rowId, Index colId) const
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > EigenvalueType
Definition: RealSchur.h:69
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
static const int m_maxIterationsPerRow
Maximum number of iterations per row.
Definition: RealSchur.h:223
#define eigen_assert(x)
Definition: Macros.h:1037
Array< int, Dynamic, 1 > v
RealScalar s
EIGEN_DEVICE_FUNC const Scalar & q
EIGEN_DEVICE_FUNC bool is_same_dense(const T1 &mat1, const T2 &mat2, typename enable_if< possibly_same_dense< T1, T2 >::value >::type *=0)
Definition: XprHelper.h:695
EIGEN_DEVICE_FUNC Derived & setConstant(Index size, const Scalar &val)
void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s &firstHouseholderVector, Scalar *workspace)
Definition: RealSchur.h:497
RealSchur & compute(const EigenBase< InputType > &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
EIGEN_DEVICE_FUNC JacobiRotation adjoint() const
Definition: Jacobi.h:67
const MatrixType & matrixT() const
Returns the quasi-triangular matrix in the Schur decomposition.
Definition: RealSchur.h:144
Matrix< Scalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > ColumnVectorType
Definition: RealSchur.h:70
Matrix< Scalar, 3, 1 > Vector3s
Definition: RealSchur.h:236
mp::number< mp::cpp_dec_float< 100 >, mp::et_on > Real
const MatrixType & matrixU() const
Returns the orthogonal matrix in the Schur decomposition.
Definition: RealSchur.h:127
RealSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: RealSchur.h:206
EIGEN_DEVICE_FUNC void evalTo(DestType &dst) const
Index findSmallSubdiagEntry(Index iu, const Scalar &considerAsZero)
Definition: RealSchur.h:376
float * p
HouseholderSequenceType matrixQ() const
Reconstructs the orthogonal matrix Q in the decomposition.
void initFrancisQRStep(Index il, Index iu, const Vector3s &shiftInfo, Index &im, Vector3s &firstHouseholderVector)
Definition: RealSchur.h:472
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics scale
void splitOffTwoRows(Index iu, bool computeU, const Scalar &exshift)
Definition: RealSchur.h:395
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
std::complex< typename NumTraits< Scalar >::Real > ComplexScalar
Definition: RealSchur.h:66
const int Dynamic
Definition: Constants.h:22
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealSchur.h:195
ColumnVectorType m_workspaceVector
Definition: RealSchur.h:229
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
RealSchur(const EigenBase< InputType > &matrix, bool computeU=true)
Constructor; computes real Schur decomposition of given matrix.
Definition: RealSchur.h:104
#define abs(x)
Definition: datatypes.h:17
ComputationInfo
Definition: Constants.h:440
EIGEN_DEVICE_FUNC bool abs2(bool x)
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:46
std::ptrdiff_t j


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:30