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 Eigen::Index Index;
67 
74  typedef std::complex<RealScalar> ComplexScalar;
75 
82 
94  explicit 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  template<typename InputType>
113  explicit ComplexSchur(const EigenBase<InputType>& matrix, bool computeU = true)
114  : m_matT(matrix.rows(),matrix.cols()),
115  m_matU(matrix.rows(),matrix.cols()),
116  m_hess(matrix.rows()),
117  m_isInitialized(false),
118  m_matUisUptodate(false),
119  m_maxIters(-1)
120  {
121  compute(matrix.derived(), computeU);
122  }
123 
138  const ComplexMatrixType& matrixU() const
139  {
140  eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
141  eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
142  return m_matU;
143  }
144 
162  const ComplexMatrixType& matrixT() const
163  {
164  eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
165  return m_matT;
166  }
167 
190  template<typename InputType>
191  ComplexSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true);
192 
210  template<typename HessMatrixType, typename OrthMatrixType>
211  ComplexSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU=true);
212 
218  {
219  eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
220  return m_info;
221  }
222 
229  {
230  m_maxIters = maxIters;
231  return *this;
232  }
233 
236  {
237  return m_maxIters;
238  }
239 
245  static const int m_maxIterationsPerRow = 30;
246 
247  protected:
248  ComplexMatrixType m_matT, m_matU;
253  Index m_maxIters;
254 
255  private:
256  bool subdiagonalEntryIsNeglegible(Index i);
257  ComplexScalar computeShift(Index iu, Index iter);
258  void reduceToTriangularForm(bool computeU);
260 };
261 
265 template<typename MatrixType>
266 inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
267 {
268  RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1));
269  RealScalar sd = numext::norm1(m_matT.coeff(i+1,i));
271  {
272  m_matT.coeffRef(i+1,i) = ComplexScalar(0);
273  return true;
274  }
275  return false;
276 }
277 
278 
280 template<typename MatrixType>
282 {
283  using std::abs;
284  if (iter == 10 || iter == 20)
285  {
286  // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
287  return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2)));
288  }
289 
290  // compute the shift as one of the eigenvalues of t, the 2x2
291  // diagonal block on the bottom of the active submatrix
292  Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
293  RealScalar normt = t.cwiseAbs().sum();
294  t /= normt; // the normalization by sf is to avoid under/overflow
295 
296  ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
297  ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
298  ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
299  ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
300  ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
301  ComplexScalar eival1 = (trace + disc) / RealScalar(2);
302  ComplexScalar eival2 = (trace - disc) / RealScalar(2);
303 
304  if(numext::norm1(eival1) > numext::norm1(eival2))
305  eival2 = det / eival1;
306  else
307  eival1 = det / eival2;
308 
309  // choose the eigenvalue closest to the bottom entry of the diagonal
310  if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1)))
311  return normt * eival1;
312  else
313  return normt * eival2;
314 }
315 
316 
317 template<typename MatrixType>
318 template<typename InputType>
320 {
321  m_matUisUptodate = false;
322  eigen_assert(matrix.cols() == matrix.rows());
323 
324  if(matrix.cols() == 1)
325  {
326  m_matT = matrix.derived().template cast<ComplexScalar>();
327  if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
328  m_info = Success;
329  m_isInitialized = true;
330  m_matUisUptodate = computeU;
331  return *this;
332  }
333 
335  computeFromHessenberg(m_matT, m_matU, computeU);
336  return *this;
337 }
338 
339 template<typename MatrixType>
340 template<typename HessMatrixType, typename OrthMatrixType>
341 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
342 {
343  m_matT = matrixH;
344  if(computeU)
345  m_matU = matrixQ;
346  reduceToTriangularForm(computeU);
347  return *this;
348 }
349 namespace internal {
350 
351 /* Reduce given matrix to Hessenberg form */
352 template<typename MatrixType, bool IsComplex>
353 struct complex_schur_reduce_to_hessenberg
354 {
355  // this is the implementation for the case IsComplex = true
356  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
357  {
358  _this.m_hess.compute(matrix);
359  _this.m_matT = _this.m_hess.matrixH();
360  if(computeU) _this.m_matU = _this.m_hess.matrixQ();
361  }
362 };
363 
364 template<typename MatrixType>
365 struct complex_schur_reduce_to_hessenberg<MatrixType, false>
366 {
367  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
368  {
369  typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
370 
371  // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
372  _this.m_hess.compute(matrix);
373  _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
374  if(computeU)
375  {
376  // This may cause an allocation which seems to be avoidable
377  MatrixType Q = _this.m_hess.matrixQ();
378  _this.m_matU = Q.template cast<ComplexScalar>();
379  }
380  }
381 };
382 
383 } // end namespace internal
384 
385 // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
386 template<typename MatrixType>
388 {
389  Index maxIters = m_maxIters;
390  if (maxIters == -1)
391  maxIters = m_maxIterationsPerRow * m_matT.rows();
392 
393  // The matrix m_matT is divided in three parts.
394  // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
395  // Rows il,...,iu is the part we are working on (the active submatrix).
396  // Rows iu+1,...,end are already brought in triangular form.
397  Index iu = m_matT.cols() - 1;
398  Index il;
399  Index iter = 0; // number of iterations we are working on the (iu,iu) element
400  Index totalIter = 0; // number of iterations for whole matrix
401 
402  while(true)
403  {
404  // find iu, the bottom row of the active submatrix
405  while(iu > 0)
406  {
407  if(!subdiagonalEntryIsNeglegible(iu-1)) break;
408  iter = 0;
409  --iu;
410  }
411 
412  // if iu is zero then we are done; the whole matrix is triangularized
413  if(iu==0) break;
414 
415  // if we spent too many iterations, we give up
416  iter++;
417  totalIter++;
418  if(totalIter > maxIters) break;
419 
420  // find il, the top row of the active submatrix
421  il = iu-1;
422  while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
423  {
424  --il;
425  }
426 
427  /* perform the QR step using Givens rotations. The first rotation
428  creates a bulge; the (il+2,il) element becomes nonzero. This
429  bulge is chased down to the bottom of the active submatrix. */
430 
431  ComplexScalar shift = computeShift(iu, iter);
433  rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
434  m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
435  m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
436  if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
437 
438  for(Index i=il+1 ; i<iu ; i++)
439  {
440  rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
441  m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
442  m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
443  m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
444  if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
445  }
446  }
447 
448  if(totalIter <= maxIters)
449  m_info = Success;
450  else
451  m_info = NoConvergence;
452 
453  m_isInitialized = true;
454  m_matUisUptodate = computeU;
455 }
456 
457 } // end namespace Eigen
458 
459 #endif // EIGEN_COMPLEX_SCHUR_H
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index cols() const
int EIGEN_BLAS_FUNC() rot(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, RealScalar *ps)
SCALAR Scalar
Definition: bench_gemm.cpp:33
EIGEN_DEVICE_FUNC bool isMuchSmallerThan(const Scalar &x, const OtherScalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
HouseholderSequenceType matrixQ() const
Reconstructs the orthogonal matrix Q in the decomposition.
MatrixHReturnType matrixH() const
Constructs the Hessenberg matrix H in the decomposition.
NumTraits< Scalar >::Real RealScalar
Definition: ComplexSchur.h:65
float real
Definition: datatypes.h:10
ComplexSchur(const EigenBase< InputType > &matrix, bool computeU=true)
Constructor; computes Schur decomposition of given matrix.
Definition: ComplexSchur.h:113
Quaternion Q
Scalar * b
Definition: benchVecAdd.cpp:17
ComplexSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: ComplexSchur.h:228
_MatrixType MatrixType
Definition: ComplexSchur.h:54
#define min(a, b)
Definition: datatypes.h:19
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index rows() const
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
EIGEN_DEVICE_FUNC const SqrtReturnType sqrt() const
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Rotation given by a cosine-sine pair.
HessenbergDecomposition< MatrixType > m_hess
Definition: ComplexSchur.h:249
iterator iter(handle obj)
Definition: pytypes.h:1547
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
ComplexMatrixType m_matU
Definition: ComplexSchur.h:248
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
ComplexSchur & compute(const EigenBase< InputType > &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
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
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
HessenbergDecomposition & compute(const EigenBase< InputType > &matrix)
Computes Hessenberg decomposition of given matrix.
static void run(ComplexSchur< MatrixType > &_this, const MatrixType &matrix, bool computeU)
Definition: ComplexSchur.h:367
ComputationInfo m_info
Definition: ComplexSchur.h:250
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
#define eigen_assert(x)
Definition: Macros.h:579
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: ComplexSchur.h:235
JacobiRotation adjoint() const
Definition: Jacobi.h:62
ComplexScalar computeShift(Index iu, Index iter)
Definition: ComplexSchur.h:281
const ComplexMatrixType & matrixT() const
Returns the triangular matrix in the Schur decomposition.
Definition: ComplexSchur.h:162
ComplexMatrixType m_matT
Definition: ComplexSchur.h:248
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
EIGEN_DEVICE_FUNC Index cols() const
Definition: EigenBase.h:62
static void run(ComplexSchur< MatrixType > &_this, const MatrixType &matrix, bool computeU)
Definition: ComplexSchur.h:356
ComplexSchur(Index size=RowsAtCompileTime==Dynamic?1:RowsAtCompileTime)
Default constructor.
Definition: ComplexSchur.h:94
void reduceToTriangularForm(bool computeU)
Definition: ComplexSchur.h:387
Eigen::Index Index
Definition: ComplexSchur.h:66
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: ComplexSchur.h:217
EIGEN_DEVICE_FUNC Index rows() const
Definition: EigenBase.h:59
const int Dynamic
Definition: Constants.h:21
EIGEN_DONT_INLINE void compute(Solver &solver, const MatrixType &A)
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Performs a complex Schur decomposition of a real or complex square matrix.
Definition: ComplexSchur.h:51
void run(Expr &expr, Dev &dev)
Definition: TensorSyclRun.h:33
#define abs(x)
Definition: datatypes.h:17
ComputationInfo
Definition: Constants.h:430
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:45
const ComplexMatrixType & matrixU() const
Returns the unitary matrix in the Schur decomposition.
Definition: ComplexSchur.h:138
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar & coeff(Index rowId, Index colId) const
Point2 t(10, 10)
void makeGivens(const Scalar &p, const Scalar &q, Scalar *r=0)
Definition: Jacobi.h:148


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:41:49