Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
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
14namespace Eigen {
15
16namespace internal {
17
18template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType;
19template<typename MatrixType>
20struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> >
21{
22 typedef MatrixType ReturnType;
23};
24
25}
26
57template<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;
75
83
86
87 typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType;
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
180 {
181 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
182 return m_hCoeffs;
183 }
184
215 {
216 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
217 return m_matrix;
218 }
219
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
270 typedef Matrix<Scalar, 1, Size, int(Options) | int(RowMajor), 1, MaxSize> VectorType;
271 typedef typename NumTraits<Scalar>::Real RealScalar;
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;
278 bool m_isInitialized;
279};
280
293template<typename MatrixType>
294void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp)
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
322namespace internal {
323
339template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
340: public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
341{
342 public:
347 HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { }
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:
367 const HessenbergDecomposition<MatrixType>& m_hess;
368};
369
370} // end namespace internal
371
372} // end namespace Eigen
373
374#endif // EIGEN_HESSENBERGDECOMPOSITION_H
Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation.
Definition: HessenbergDecomposition.h:58
HessenbergDecomposition & compute(const EigenBase< InputType > &matrix)
Computes Hessenberg decomposition of given matrix.
Definition: HessenbergDecomposition.h:152
HouseholderSequenceType matrixQ() const
Reconstructs the orthogonal matrix Q in the decomposition.
Definition: HessenbergDecomposition.h:234
const MatrixType & packedMatrix() const
Returns the internal representation of the decomposition.
Definition: HessenbergDecomposition.h:214
Matrix< Scalar, SizeMinusOne, 1, Options &~RowMajor, MaxSizeMinusOne, 1 > CoeffVectorType
Type for vector of Householder coefficients.
Definition: HessenbergDecomposition.h:82
HouseholderSequence< MatrixType, typename internal::remove_all< typename CoeffVectorType::ConjugateReturnType >::type > HouseholderSequenceType
Return type of matrixQ()
Definition: HessenbergDecomposition.h:85
Eigen::Index Index
Definition: HessenbergDecomposition.h:74
MatrixHReturnType matrixH() const
Constructs the Hessenberg matrix H in the decomposition.
Definition: HessenbergDecomposition.h:262
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
Definition: HessenbergDecomposition.h:62
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: HessenbergDecomposition.h:73
HessenbergDecomposition(const EigenBase< InputType > &matrix)
Constructor; computes Hessenberg decomposition of given matrix.
Definition: HessenbergDecomposition.h:119
HessenbergDecomposition(Index size=Size==Dynamic ? 2 :Size)
Default constructor; the decomposition will be computed later.
Definition: HessenbergDecomposition.h:100
const CoeffVectorType & householderCoefficients() const
Returns the Householder coefficients.
Definition: HessenbergDecomposition.h:179
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:121
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
@ RowMajor
Definition: Constants.h:321
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const int Dynamic
Definition: Constants.h:22
Definition: EigenBase.h:30
Derived & derived()
Definition: EigenBase.h:46
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60