Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
HessenbergDecomposition.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2010 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_HESSENBERGDECOMPOSITION_H
12 #define EIGEN_HESSENBERGDECOMPOSITION_H
13 
14 #include "./InternalHeaderCheck.h"
15 
16 namespace Eigen {
17 
18 namespace internal {
19 
20 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType;
21 template<typename MatrixType>
22 struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> >
23 {
24  typedef MatrixType ReturnType;
25 };
26 
27 }
28 
59 template<typename MatrixType_> class HessenbergDecomposition
60 {
61  public:
62 
64  typedef MatrixType_ MatrixType;
65 
66  enum {
67  Size = MatrixType::RowsAtCompileTime,
68  SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
69  Options = MatrixType::Options,
70  MaxSize = MatrixType::MaxRowsAtCompileTime,
71  MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
72  };
73 
75  typedef typename MatrixType::Scalar Scalar;
76  typedef Eigen::Index Index;
77 
85 
88 
89  typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType;
90 
102  explicit HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size)
103  : m_matrix(size,size),
104  m_temp(size),
105  m_isInitialized(false)
106  {
107  if(size>1)
108  m_hCoeffs.resize(size-1);
109  }
110 
120  template<typename InputType>
122  : m_matrix(matrix.derived()),
123  m_temp(matrix.rows()),
124  m_isInitialized(false)
125  {
126  if(matrix.rows()<2)
127  {
128  m_isInitialized = true;
129  return;
130  }
131  m_hCoeffs.resize(matrix.rows()-1,1);
132  _compute(m_matrix, m_hCoeffs, m_temp);
133  m_isInitialized = true;
134  }
135 
153  template<typename InputType>
155  {
156  m_matrix = matrix.derived();
157  if(matrix.rows()<2)
158  {
159  m_isInitialized = true;
160  return *this;
161  }
162  m_hCoeffs.resize(matrix.rows()-1,1);
163  _compute(m_matrix, m_hCoeffs, m_temp);
164  m_isInitialized = true;
165  return *this;
166  }
167 
182  {
183  eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
184  return m_hCoeffs;
185  }
186 
216  const MatrixType& packedMatrix() const
217  {
218  eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
219  return m_matrix;
220  }
221 
237  {
238  eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
239  return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
240  .setLength(m_matrix.rows() - 1)
241  .setShift(1);
242  }
243 
264  MatrixHReturnType matrixH() const
265  {
266  eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
267  return MatrixHReturnType(*this);
268  }
269 
270  private:
271 
272  typedef Matrix<Scalar, 1, Size, int(Options) | int(RowMajor), 1, MaxSize> VectorType;
273  typedef typename NumTraits<Scalar>::Real RealScalar;
274  static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
275 
276  protected:
277  MatrixType m_matrix;
278  CoeffVectorType m_hCoeffs;
279  VectorType m_temp;
280  bool m_isInitialized;
281 };
282 
295 template<typename MatrixType>
296 void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp)
297 {
298  eigen_assert(matA.rows()==matA.cols());
299  Index n = matA.rows();
300  temp.resize(n);
301  for (Index i = 0; i<n-1; ++i)
302  {
303  // let's consider the vector v = i-th column starting at position i+1
304  Index remainingSize = n-i-1;
305  RealScalar beta;
306  Scalar h;
307  matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
308  matA.col(i).coeffRef(i+1) = beta;
309  hCoeffs.coeffRef(i) = h;
310 
311  // Apply similarity transformation to remaining columns,
312  // i.e., compute A = H A H'
313 
314  // A = H A
315  matA.bottomRightCorner(remainingSize, remainingSize)
316  .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
317 
318  // A = A H'
319  matA.rightCols(remainingSize)
320  .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1), numext::conj(h), &temp.coeffRef(0));
321  }
322 }
323 
324 namespace internal {
325 
341 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
342 : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
343 {
344  public:
349  HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { }
350 
356  template <typename ResultType>
357  inline void evalTo(ResultType& result) const
358  {
359  result = m_hess.packedMatrix();
360  Index n = result.rows();
361  if (n>2)
362  result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero();
363  }
364 
365  Index rows() const { return m_hess.packedMatrix().rows(); }
366  Index cols() const { return m_hess.packedMatrix().cols(); }
367 
368  protected:
369  const HessenbergDecomposition<MatrixType>& m_hess;
370 };
371 
372 } // end namespace internal
373 
374 } // end namespace Eigen
375 
376 #endif // EIGEN_HESSENBERGDECOMPOSITION_H
Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation.
Definition: HessenbergDecomposition.h:60
Matrix< Scalar, SizeMinusOne, 1, Options &~RowMajor, MaxSizeMinusOne, 1 > CoeffVectorType
Type for vector of Householder coefficients.
Definition: HessenbergDecomposition.h:84
const CoeffVectorType & householderCoefficients() const
Returns the Householder coefficients.
Definition: HessenbergDecomposition.h:181
HessenbergDecomposition(Index size=Size==Dynamic ? 2 :Size)
Default constructor; the decomposition will be computed later.
Definition: HessenbergDecomposition.h:102
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: HessenbergDecomposition.h:75
HessenbergDecomposition & compute(const EigenBase< InputType > &matrix)
Computes Hessenberg decomposition of given matrix.
Definition: HessenbergDecomposition.h:154
HessenbergDecomposition(const EigenBase< InputType > &matrix)
Constructor; computes Hessenberg decomposition of given matrix.
Definition: HessenbergDecomposition.h:121
MatrixHReturnType matrixH() const
Constructs the Hessenberg matrix H in the decomposition.
Definition: HessenbergDecomposition.h:264
HouseholderSequenceType matrixQ() const
Reconstructs the orthogonal matrix Q in the decomposition.
Definition: HessenbergDecomposition.h:236
const MatrixType & packedMatrix() const
Returns the internal representation of the decomposition.
Definition: HessenbergDecomposition.h:216
HouseholderSequence< MatrixType, internal::remove_all_t< typename CoeffVectorType::ConjugateReturnType > > HouseholderSequenceType
Return type of matrixQ()
Definition: HessenbergDecomposition.h:87
Eigen::Index Index
Definition: HessenbergDecomposition.h:76
MatrixType_ MatrixType
Synonym for the template parameter MatrixType_.
Definition: HessenbergDecomposition.h:64
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:123
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:283
@ RowMajor
Definition: Constants.h:323
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const int Dynamic
Definition: Constants.h:24
Definition: EigenBase.h:32
Derived & derived()
Definition: EigenBase.h:48
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:62
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231