ComplexSchur.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) 2009 Claire Maurice
5 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
6 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
7 //
8 // This Source Code Form is subject to the terms of the Mozilla
9 // Public License v. 2.0. If a copy of the MPL was not distributed
10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
11 
12 #ifndef EIGEN_COMPLEX_SCHUR_H
13 #define EIGEN_COMPLEX_SCHUR_H
14 
16 
17 namespace Eigen {
18 
19 namespace internal {
20 template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg;
21 }
22 
51 template<typename _MatrixType> class ComplexSchur
52 {
53  public:
54  typedef _MatrixType MatrixType;
55  enum {
56  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
57  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
58  Options = MatrixType::Options,
59  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
60  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
61  };
62 
64  typedef typename MatrixType::Scalar Scalar;
66  typedef typename MatrixType::Index Index;
67 
74  typedef std::complex<RealScalar> ComplexScalar;
75 
82 
94  ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
95  : m_matT(size,size),
96  m_matU(size,size),
97  m_hess(size),
98  m_isInitialized(false),
99  m_matUisUptodate(false),
100  m_maxIters(-1)
101  {}
102 
112  ComplexSchur(const MatrixType& matrix, bool computeU = true)
113  : m_matT(matrix.rows(),matrix.cols()),
114  m_matU(matrix.rows(),matrix.cols()),
115  m_hess(matrix.rows()),
116  m_isInitialized(false),
117  m_matUisUptodate(false),
118  m_maxIters(-1)
119  {
120  compute(matrix, computeU);
121  }
122 
137  const ComplexMatrixType& matrixU() const
138  {
139  eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
140  eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
141  return m_matU;
142  }
143 
161  const ComplexMatrixType& matrixT() const
162  {
163  eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
164  return m_matT;
165  }
166 
189  ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
190 
208  template<typename HessMatrixType, typename OrthMatrixType>
209  ComplexSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU=true);
210 
216  {
217  eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
218  return m_info;
219  }
220 
227  {
228  m_maxIters = maxIters;
229  return *this;
230  }
231 
234  {
235  return m_maxIters;
236  }
237 
243  static const int m_maxIterationsPerRow = 30;
244 
245  protected:
246  ComplexMatrixType m_matT, m_matU;
251  Index m_maxIters;
252 
253  private:
254  bool subdiagonalEntryIsNeglegible(Index i);
255  ComplexScalar computeShift(Index iu, Index iter);
256  void reduceToTriangularForm(bool computeU);
257  friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
258 };
259 
263 template<typename MatrixType>
264 inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
265 {
266  RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1));
267  RealScalar sd = numext::norm1(m_matT.coeff(i+1,i));
269  {
270  m_matT.coeffRef(i+1,i) = ComplexScalar(0);
271  return true;
272  }
273  return false;
274 }
275 
276 
278 template<typename MatrixType>
280 {
281  using std::abs;
282  if (iter == 10 || iter == 20)
283  {
284  // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
285  return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2)));
286  }
287 
288  // compute the shift as one of the eigenvalues of t, the 2x2
289  // diagonal block on the bottom of the active submatrix
290  Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
291  RealScalar normt = t.cwiseAbs().sum();
292  t /= normt; // the normalization by sf is to avoid under/overflow
293 
294  ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
295  ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
296  ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
297  ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
298  ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
299  ComplexScalar eival1 = (trace + disc) / RealScalar(2);
300  ComplexScalar eival2 = (trace - disc) / RealScalar(2);
301 
302  if(numext::norm1(eival1) > numext::norm1(eival2))
303  eival2 = det / eival1;
304  else
305  eival1 = det / eival2;
306 
307  // choose the eigenvalue closest to the bottom entry of the diagonal
308  if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1)))
309  return normt * eival1;
310  else
311  return normt * eival2;
312 }
313 
314 
315 template<typename MatrixType>
316 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
317 {
318  m_matUisUptodate = false;
319  eigen_assert(matrix.cols() == matrix.rows());
320 
321  if(matrix.cols() == 1)
322  {
323  m_matT = matrix.template cast<ComplexScalar>();
324  if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
325  m_info = Success;
326  m_isInitialized = true;
327  m_matUisUptodate = computeU;
328  return *this;
329  }
330 
332  computeFromHessenberg(m_matT, m_matU, computeU);
333  return *this;
334 }
335 
336 template<typename MatrixType>
337 template<typename HessMatrixType, typename OrthMatrixType>
338 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
339 {
340  m_matT = matrixH;
341  if(computeU)
342  m_matU = matrixQ;
343  reduceToTriangularForm(computeU);
344  return *this;
345 }
346 namespace internal {
347 
348 /* Reduce given matrix to Hessenberg form */
349 template<typename MatrixType, bool IsComplex>
350 struct complex_schur_reduce_to_hessenberg
351 {
352  // this is the implementation for the case IsComplex = true
353  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
354  {
355  _this.m_hess.compute(matrix);
356  _this.m_matT = _this.m_hess.matrixH();
357  if(computeU) _this.m_matU = _this.m_hess.matrixQ();
358  }
359 };
360 
361 template<typename MatrixType>
362 struct complex_schur_reduce_to_hessenberg<MatrixType, false>
363 {
364  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
365  {
366  typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
367 
368  // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
369  _this.m_hess.compute(matrix);
370  _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
371  if(computeU)
372  {
373  // This may cause an allocation which seems to be avoidable
374  MatrixType Q = _this.m_hess.matrixQ();
375  _this.m_matU = Q.template cast<ComplexScalar>();
376  }
377  }
378 };
379 
380 } // end namespace internal
381 
382 // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
383 template<typename MatrixType>
385 {
386  Index maxIters = m_maxIters;
387  if (maxIters == -1)
388  maxIters = m_maxIterationsPerRow * m_matT.rows();
389 
390  // The matrix m_matT is divided in three parts.
391  // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
392  // Rows il,...,iu is the part we are working on (the active submatrix).
393  // Rows iu+1,...,end are already brought in triangular form.
394  Index iu = m_matT.cols() - 1;
395  Index il;
396  Index iter = 0; // number of iterations we are working on the (iu,iu) element
397  Index totalIter = 0; // number of iterations for whole matrix
398 
399  while(true)
400  {
401  // find iu, the bottom row of the active submatrix
402  while(iu > 0)
403  {
404  if(!subdiagonalEntryIsNeglegible(iu-1)) break;
405  iter = 0;
406  --iu;
407  }
408 
409  // if iu is zero then we are done; the whole matrix is triangularized
410  if(iu==0) break;
411 
412  // if we spent too many iterations, we give up
413  iter++;
414  totalIter++;
415  if(totalIter > maxIters) break;
416 
417  // find il, the top row of the active submatrix
418  il = iu-1;
419  while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
420  {
421  --il;
422  }
423 
424  /* perform the QR step using Givens rotations. The first rotation
425  creates a bulge; the (il+2,il) element becomes nonzero. This
426  bulge is chased down to the bottom of the active submatrix. */
427 
428  ComplexScalar shift = computeShift(iu, iter);
430  rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
431  m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
432  m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
433  if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
434 
435  for(Index i=il+1 ; i<iu ; i++)
436  {
437  rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
438  m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
439  m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
440  m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
441  if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
442  }
443  }
444 
445  if(totalIter <= maxIters)
446  m_info = Success;
447  else
448  m_info = NoConvergence;
449 
450  m_isInitialized = true;
451  m_matUisUptodate = computeU;
452 }
453 
454 } // end namespace Eigen
455 
456 #endif // EIGEN_COMPLEX_SCHUR_H
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.
NumTraits< Scalar >::Real RealScalar
Definition: ComplexSchur.h:65
#define inline
Definition: Types.h:109
ComplexSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: ComplexSchur.h:226
_MatrixType MatrixType
Definition: ComplexSchur.h:54
void makeGivens(const Scalar &p, const Scalar &q, Scalar *z=0)
Definition: Jacobi.h:148
ComplexSchur & compute(const MatrixType &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
Definition: ComplexSchur.h:316
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
Rotation given by a cosine-sine pair.
HessenbergDecomposition< MatrixType > m_hess
Definition: ComplexSchur.h:247
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
ComplexMatrixType m_matU
Definition: ComplexSchur.h:246
bool isMuchSmallerThan(const Scalar &x, const OtherScalar &y, typename NumTraits< Scalar >::Real precision=NumTraits< Scalar >::dummy_precision())
ComplexSchur & computeFromHessenberg(const HessMatrixType &matrixH, const OrthMatrixType &matrixQ, bool computeU=true)
Compute Schur decomposition from a given Hessenberg matrix.
std::complex< RealScalar > ComplexScalar
Complex scalar type for _MatrixType.
Definition: ComplexSchur.h:74
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > ComplexMatrixType
Type for the matrices in the Schur decomposition.
Definition: ComplexSchur.h:81
MatrixType::Scalar Scalar
Scalar type for matrices of type _MatrixType.
Definition: ComplexSchur.h:64
EIGEN_STRONG_INLINE Index rows() const
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
RealReturnType real() const
static void run(ComplexSchur< MatrixType > &_this, const MatrixType &matrix, bool computeU)
Definition: ComplexSchur.h:364
EIGEN_STRONG_INLINE const Scalar & coeff(Index rowId, Index colId) const
ComputationInfo m_info
Definition: ComplexSchur.h:248
HessenbergDecomposition & compute(const MatrixType &matrix)
Computes Hessenberg decomposition of given matrix.
EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
Provides a generic way to set and pass user-specified options.
Definition: options.hpp:65
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: ComplexSchur.h:233
JacobiRotation adjoint() const
Definition: Jacobi.h:62
ComplexScalar computeShift(Index iu, Index iter)
Definition: ComplexSchur.h:279
const ComplexMatrixType & matrixT() const
Returns the triangular matrix in the Schur decomposition.
Definition: ComplexSchur.h:161
ComplexMatrixType m_matT
Definition: ComplexSchur.h:246
static void run(ComplexSchur< MatrixType > &_this, const MatrixType &matrix, bool computeU)
Definition: ComplexSchur.h:353
ComplexSchur(Index size=RowsAtCompileTime==Dynamic?1:RowsAtCompileTime)
Default constructor.
Definition: ComplexSchur.h:94
ComplexSchur(const MatrixType &matrix, bool computeU=true)
Constructor; computes Schur decomposition of given matrix.
Definition: ComplexSchur.h:112
void reduceToTriangularForm(bool computeU)
Definition: ComplexSchur.h:384
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: ComplexSchur.h:215
MatrixType::Index Index
Definition: ComplexSchur.h:66
EIGEN_STRONG_INLINE Index cols() const
Performs a complex Schur decomposition of a real or complex square matrix.
Definition: ComplexSchur.h:51
#define eigen_assert(x)
const ComplexMatrixType & matrixU() const
Returns the unitary matrix in the Schur decomposition.
Definition: ComplexSchur.h:137


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