HessenbergDecomposition.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-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 namespace Eigen {
15 
16 namespace internal {
17 
18 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType;
19 template<typename MatrixType>
21 {
23 };
24 
25 }
26 
57 template<typename _MatrixType> class HessenbergDecomposition
58 {
59  public:
60 
62  typedef _MatrixType MatrixType;
63 
64  enum {
65  Size = MatrixType::RowsAtCompileTime,
66  SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
67  Options = MatrixType::Options,
68  MaxSize = MatrixType::MaxRowsAtCompileTime,
69  MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
70  };
71 
73  typedef typename MatrixType::Scalar Scalar;
74  typedef Eigen::Index Index;
75 
83 
86 
88 
100  explicit HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size)
101  : m_matrix(size,size),
102  m_temp(size),
103  m_isInitialized(false)
104  {
105  if(size>1)
106  m_hCoeffs.resize(size-1);
107  }
108 
118  template<typename InputType>
120  : m_matrix(matrix.derived()),
121  m_temp(matrix.rows()),
122  m_isInitialized(false)
123  {
124  if(matrix.rows()<2)
125  {
126  m_isInitialized = true;
127  return;
128  }
129  m_hCoeffs.resize(matrix.rows()-1,1);
130  _compute(m_matrix, m_hCoeffs, m_temp);
131  m_isInitialized = true;
132  }
133 
151  template<typename InputType>
153  {
154  m_matrix = matrix.derived();
155  if(matrix.rows()<2)
156  {
157  m_isInitialized = true;
158  return *this;
159  }
160  m_hCoeffs.resize(matrix.rows()-1,1);
161  _compute(m_matrix, m_hCoeffs, m_temp);
162  m_isInitialized = true;
163  return *this;
164  }
165 
179  const CoeffVectorType& householderCoefficients() const
180  {
181  eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
182  return m_hCoeffs;
183  }
184 
214  const MatrixType& packedMatrix() const
215  {
216  eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
217  return m_matrix;
218  }
219 
234  HouseholderSequenceType matrixQ() const
235  {
236  eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
237  return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
238  .setLength(m_matrix.rows() - 1)
239  .setShift(1);
240  }
241 
262  MatrixHReturnType matrixH() const
263  {
264  eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
265  return MatrixHReturnType(*this);
266  }
267 
268  private:
269 
272  static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
273 
274  protected:
275  MatrixType m_matrix;
276  CoeffVectorType m_hCoeffs;
277  VectorType m_temp;
279 };
280 
293 template<typename MatrixType>
295 {
296  eigen_assert(matA.rows()==matA.cols());
297  Index n = matA.rows();
298  temp.resize(n);
299  for (Index i = 0; i<n-1; ++i)
300  {
301  // let's consider the vector v = i-th column starting at position i+1
302  Index remainingSize = n-i-1;
303  RealScalar beta;
304  Scalar h;
305  matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
306  matA.col(i).coeffRef(i+1) = beta;
307  hCoeffs.coeffRef(i) = h;
308 
309  // Apply similarity transformation to remaining columns,
310  // i.e., compute A = H A H'
311 
312  // A = H A
313  matA.bottomRightCorner(remainingSize, remainingSize)
314  .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
315 
316  // A = A H'
317  matA.rightCols(remainingSize)
318  .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1), numext::conj(h), &temp.coeffRef(0));
319  }
320 }
321 
322 namespace internal {
323 
339 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
340 : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
341 {
342  public:
348 
354  template <typename ResultType>
355  inline void evalTo(ResultType& result) const
356  {
357  result = m_hess.packedMatrix();
358  Index n = result.rows();
359  if (n>2)
360  result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero();
361  }
362 
363  Index rows() const { return m_hess.packedMatrix().rows(); }
364  Index cols() const { return m_hess.packedMatrix().cols(); }
365 
366  protected:
368 };
369 
370 } // end namespace internal
371 
372 } // end namespace Eigen
373 
374 #endif // EIGEN_HESSENBERGDECOMPOSITION_H
SCALAR Scalar
Definition: bench_gemm.cpp:46
const MatrixType & packedMatrix() const
Returns the internal representation of the decomposition.
static void _compute(MatrixType &matA, CoeffVectorType &hCoeffs, VectorType &temp)
int n
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
MatrixXf MatrixType
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
internal::HessenbergDecompositionMatrixHReturnType< MatrixType > MatrixHReturnType
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
void evalTo(ResultType &result) const
Hessenberg matrix in decomposition.
AnnoyingScalar conj(const AnnoyingScalar &x)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
Sequence of Householder reflections acting on subspaces with decreasing size.
HessenbergDecomposition & compute(const EigenBase< InputType > &matrix)
Computes Hessenberg decomposition of given matrix.
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Values result
NumTraits< Scalar >::Real RealScalar
MatrixHReturnType matrixH() const
Constructs the Hessenberg matrix H in the decomposition.
HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition< MatrixType > &hess)
Constructor.
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
const HessenbergDecomposition< MatrixType > & m_hess
HessenbergDecomposition(const EigenBase< InputType > &matrix)
Constructor; computes Hessenberg decomposition of given matrix.
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
MatrixXf matA(2, 2)
Expression type for return value of HessenbergDecomposition::matrixH()
EIGEN_CONSTEXPR Index size(const T &x)
Definition: Meta.h:479
HessenbergDecomposition(Index size=Size==Dynamic ? 2 :Size)
Default constructor; the decomposition will be computed later.
const double h
A triangularView< Lower >().adjoint().solveInPlace(B)
Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation.
Matrix< Scalar, 1, Size, int(Options)|int(RowMajor), 1, MaxSize > VectorType
HouseholderSequenceType matrixQ() const
Reconstructs the orthogonal matrix Q in the decomposition.
const CoeffVectorType & householderCoefficients() const
Returns the Householder coefficients.
const int Dynamic
Definition: Constants.h:22
HouseholderSequence< MatrixType, typename internal::remove_all< typename CoeffVectorType::ConjugateReturnType >::type > HouseholderSequenceType
Return type of matrixQ()
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Matrix< Scalar, SizeMinusOne, 1, Options &~RowMajor, MaxSizeMinusOne, 1 > CoeffVectorType
Type for vector of Householder coefficients.
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:46
v setZero(3)


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