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  RealScalar eival1_norm = numext::norm1(eival1);
304  RealScalar eival2_norm = numext::norm1(eival2);
305  // A division by zero can only occur if eival1==eival2==0.
306  // In this case, det==0, and all we have to do is checking that eival2_norm!=0
307  if(eival1_norm > eival2_norm)
308  eival2 = det / eival1;
309  else if(eival2_norm!=RealScalar(0))
310  eival1 = det / eival2;
311 
312  // choose the eigenvalue closest to the bottom entry of the diagonal
313  if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1)))
314  return normt * eival1;
315  else
316  return normt * eival2;
317 }
318 
319 
320 template<typename MatrixType>
321 template<typename InputType>
323 {
324  m_matUisUptodate = false;
325  eigen_assert(matrix.cols() == matrix.rows());
326 
327  if(matrix.cols() == 1)
328  {
329  m_matT = matrix.derived().template cast<ComplexScalar>();
330  if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
331  m_info = Success;
332  m_isInitialized = true;
333  m_matUisUptodate = computeU;
334  return *this;
335  }
336 
338  computeFromHessenberg(m_matT, m_matU, computeU);
339  return *this;
340 }
341 
342 template<typename MatrixType>
343 template<typename HessMatrixType, typename OrthMatrixType>
344 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
345 {
346  m_matT = matrixH;
347  if(computeU)
348  m_matU = matrixQ;
349  reduceToTriangularForm(computeU);
350  return *this;
351 }
352 namespace internal {
353 
354 /* Reduce given matrix to Hessenberg form */
355 template<typename MatrixType, bool IsComplex>
356 struct complex_schur_reduce_to_hessenberg
357 {
358  // this is the implementation for the case IsComplex = true
359  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
360  {
361  _this.m_hess.compute(matrix);
362  _this.m_matT = _this.m_hess.matrixH();
363  if(computeU) _this.m_matU = _this.m_hess.matrixQ();
364  }
365 };
366 
367 template<typename MatrixType>
368 struct complex_schur_reduce_to_hessenberg<MatrixType, false>
369 {
370  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
371  {
372  typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
373 
374  // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
375  _this.m_hess.compute(matrix);
376  _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
377  if(computeU)
378  {
379  // This may cause an allocation which seems to be avoidable
380  MatrixType Q = _this.m_hess.matrixQ();
381  _this.m_matU = Q.template cast<ComplexScalar>();
382  }
383  }
384 };
385 
386 } // end namespace internal
387 
388 // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
389 template<typename MatrixType>
391 {
392  Index maxIters = m_maxIters;
393  if (maxIters == -1)
394  maxIters = m_maxIterationsPerRow * m_matT.rows();
395 
396  // The matrix m_matT is divided in three parts.
397  // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
398  // Rows il,...,iu is the part we are working on (the active submatrix).
399  // Rows iu+1,...,end are already brought in triangular form.
400  Index iu = m_matT.cols() - 1;
401  Index il;
402  Index iter = 0; // number of iterations we are working on the (iu,iu) element
403  Index totalIter = 0; // number of iterations for whole matrix
404 
405  while(true)
406  {
407  // find iu, the bottom row of the active submatrix
408  while(iu > 0)
409  {
410  if(!subdiagonalEntryIsNeglegible(iu-1)) break;
411  iter = 0;
412  --iu;
413  }
414 
415  // if iu is zero then we are done; the whole matrix is triangularized
416  if(iu==0) break;
417 
418  // if we spent too many iterations, we give up
419  iter++;
420  totalIter++;
421  if(totalIter > maxIters) break;
422 
423  // find il, the top row of the active submatrix
424  il = iu-1;
425  while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
426  {
427  --il;
428  }
429 
430  /* perform the QR step using Givens rotations. The first rotation
431  creates a bulge; the (il+2,il) element becomes nonzero. This
432  bulge is chased down to the bottom of the active submatrix. */
433 
434  ComplexScalar shift = computeShift(iu, iter);
436  rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
437  m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
438  m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
439  if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
440 
441  for(Index i=il+1 ; i<iu ; i++)
442  {
443  rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
444  m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
445  m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
446  m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
447  if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
448  }
449  }
450 
451  if(totalIter <= maxIters)
452  m_info = Success;
453  else
454  m_info = NoConvergence;
455 
456  m_isInitialized = true;
457  m_matUisUptodate = computeU;
458 }
459 
460 } // end namespace Eigen
461 
462 #endif // EIGEN_COMPLEX_SCHUR_H
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:46
EIGEN_DEVICE_FUNC bool isMuchSmallerThan(const Scalar &x, const OtherScalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
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
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
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:2273
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
ComplexMatrixType m_matU
Definition: ComplexSchur.h:248
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:63
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
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.
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: ComplexSchur.h:217
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
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:370
ComputationInfo m_info
Definition: ComplexSchur.h:250
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
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar & coeff(Index rowId, Index colId) const
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
#define eigen_assert(x)
Definition: Macros.h:1037
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: ComplexSchur.h:235
ComplexSchur(Index size=RowsAtCompileTime==Dynamic ? 1 :RowsAtCompileTime)
Default constructor.
Definition: ComplexSchur.h:94
ComplexScalar computeShift(Index iu, Index iter)
Definition: ComplexSchur.h:281
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
ComplexMatrixType m_matT
Definition: ComplexSchur.h:248
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
EIGEN_DEVICE_FUNC JacobiRotation adjoint() const
Definition: Jacobi.h:67
static void run(ComplexSchur< MatrixType > &_this, const MatrixType &matrix, bool computeU)
Definition: ComplexSchur.h:359
EIGEN_CONSTEXPR Index size(const T &x)
Definition: Meta.h:479
void reduceToTriangularForm(bool computeU)
Definition: ComplexSchur.h:390
const ComplexMatrixType & matrixU() const
Returns the unitary matrix in the Schur decomposition.
Definition: ComplexSchur.h:138
Eigen::Index Index
Definition: ComplexSchur.h:66
HouseholderSequenceType matrixQ() const
Reconstructs the orthogonal matrix Q in the decomposition.
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
const int Dynamic
Definition: Constants.h:22
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
#define abs(x)
Definition: datatypes.h:17
ComputationInfo
Definition: Constants.h:440
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:46
const ComplexMatrixType & matrixT() const
Returns the triangular matrix in the Schur decomposition.
Definition: ComplexSchur.h:162
Point2 t(10, 10)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:03