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
HouseholderSequence.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
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_HOUSEHOLDER_SEQUENCE_H
12#define EIGEN_HOUSEHOLDER_SEQUENCE_H
13
14namespace Eigen {
15
57namespace internal {
58
59template<typename VectorsType, typename CoeffsType, int Side>
60struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
61{
62 typedef typename VectorsType::Scalar Scalar;
63 typedef typename VectorsType::StorageIndex StorageIndex;
64 typedef typename VectorsType::StorageKind StorageKind;
65 enum {
66 RowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::RowsAtCompileTime
67 : traits<VectorsType>::ColsAtCompileTime,
68 ColsAtCompileTime = RowsAtCompileTime,
69 MaxRowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::MaxRowsAtCompileTime
70 : traits<VectorsType>::MaxColsAtCompileTime,
71 MaxColsAtCompileTime = MaxRowsAtCompileTime,
72 Flags = 0
73 };
74};
75
76struct HouseholderSequenceShape {};
77
78template<typename VectorsType, typename CoeffsType, int Side>
79struct evaluator_traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
80 : public evaluator_traits_base<HouseholderSequence<VectorsType,CoeffsType,Side> >
81{
82 typedef HouseholderSequenceShape Shape;
83};
84
85template<typename VectorsType, typename CoeffsType, int Side>
86struct hseq_side_dependent_impl
87{
88 typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType;
89 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType;
90 static EIGEN_DEVICE_FUNC inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
91 {
92 Index start = k+1+h.m_shift;
93 return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
94 }
95};
96
97template<typename VectorsType, typename CoeffsType>
98struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight>
99{
100 typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType;
101 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType;
102 static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
103 {
104 Index start = k+1+h.m_shift;
105 return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose();
106 }
107};
108
109template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type
110{
111 typedef typename ScalarBinaryOpTraits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
112 ResultScalar;
113 typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
114 0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
115};
116
117} // end namespace internal
118
119template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence
120 : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
121{
123
124 public:
125 enum {
126 RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
127 ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
128 MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
129 MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
130 };
131 typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
132
133 typedef HouseholderSequence<
134 typename internal::conditional<NumTraits<Scalar>::IsComplex,
135 typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type,
136 VectorsType>::type,
137 typename internal::conditional<NumTraits<Scalar>::IsComplex,
138 typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
139 CoeffsType>::type,
140 Side
142
143 typedef HouseholderSequence<
144 VectorsType,
145 typename internal::conditional<NumTraits<Scalar>::IsComplex,
146 typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
147 CoeffsType>::type,
148 Side
150
151 typedef HouseholderSequence<
152 typename internal::conditional<NumTraits<Scalar>::IsComplex,
153 typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type,
154 VectorsType>::type,
155 CoeffsType,
156 Side
158
159 typedef HouseholderSequence<
162 Side
164
182 EIGEN_DEVICE_FUNC
183 HouseholderSequence(const VectorsType& v, const CoeffsType& h)
184 : m_vectors(v), m_coeffs(h), m_reverse(false), m_length(v.diagonalSize()),
185 m_shift(0)
186 {
187 }
188
190 EIGEN_DEVICE_FUNC
192 : m_vectors(other.m_vectors),
193 m_coeffs(other.m_coeffs),
194 m_reverse(other.m_reverse),
195 m_length(other.m_length),
196 m_shift(other.m_shift)
197 {
198 }
199
204 EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
205 Index rows() const EIGEN_NOEXCEPT { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
206
211 EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
212 Index cols() const EIGEN_NOEXCEPT { return rows(); }
213
228 EIGEN_DEVICE_FUNC
230 {
231 eigen_assert(k >= 0 && k < m_length);
232 return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
233 }
234
237 {
238 return TransposeReturnType(m_vectors.conjugate(), m_coeffs)
239 .setReverseFlag(!m_reverse)
240 .setLength(m_length)
241 .setShift(m_shift);
242 }
243
246 {
247 return ConjugateReturnType(m_vectors.conjugate(), m_coeffs.conjugate())
248 .setReverseFlag(m_reverse)
249 .setLength(m_length)
250 .setShift(m_shift);
251 }
252
256 template<bool Cond>
257 EIGEN_DEVICE_FUNC
258 inline typename internal::conditional<Cond,ConjugateReturnType,ConstHouseholderSequence>::type
260 {
261 typedef typename internal::conditional<Cond,ConjugateReturnType,ConstHouseholderSequence>::type ReturnType;
262 return ReturnType(m_vectors.template conjugateIf<Cond>(), m_coeffs.template conjugateIf<Cond>());
263 }
264
267 {
268 return AdjointReturnType(m_vectors, m_coeffs.conjugate())
269 .setReverseFlag(!m_reverse)
270 .setLength(m_length)
271 .setShift(m_shift);
272 }
273
275 AdjointReturnType inverse() const { return adjoint(); }
276
278 template<typename DestType>
279 inline EIGEN_DEVICE_FUNC
280 void evalTo(DestType& dst) const
281 {
282 Matrix<Scalar, DestType::RowsAtCompileTime, 1,
283 AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows());
284 evalTo(dst, workspace);
285 }
286
288 template<typename Dest, typename Workspace>
289 EIGEN_DEVICE_FUNC
290 void evalTo(Dest& dst, Workspace& workspace) const
291 {
292 workspace.resize(rows());
293 Index vecs = m_length;
294 if(internal::is_same_dense(dst,m_vectors))
295 {
296 // in-place
297 dst.diagonal().setOnes();
298 dst.template triangularView<StrictlyUpper>().setZero();
299 for(Index k = vecs-1; k >= 0; --k)
300 {
301 Index cornerSize = rows() - k - m_shift;
302 if(m_reverse)
303 dst.bottomRightCorner(cornerSize, cornerSize)
304 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
305 else
306 dst.bottomRightCorner(cornerSize, cornerSize)
307 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
308
309 // clear the off diagonal vector
310 dst.col(k).tail(rows()-k-1).setZero();
311 }
312 // clear the remaining columns if needed
313 for(Index k = 0; k<cols()-vecs ; ++k)
314 dst.col(k).tail(rows()-k-1).setZero();
315 }
316 else if(m_length>BlockSize)
317 {
318 dst.setIdentity(rows(), rows());
319 if(m_reverse)
320 applyThisOnTheLeft(dst,workspace,true);
321 else
322 applyThisOnTheLeft(dst,workspace,true);
323 }
324 else
325 {
326 dst.setIdentity(rows(), rows());
327 for(Index k = vecs-1; k >= 0; --k)
328 {
329 Index cornerSize = rows() - k - m_shift;
330 if(m_reverse)
331 dst.bottomRightCorner(cornerSize, cornerSize)
332 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
333 else
334 dst.bottomRightCorner(cornerSize, cornerSize)
335 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
336 }
337 }
338 }
339
341 template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
342 {
343 Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
344 applyThisOnTheRight(dst, workspace);
345 }
346
348 template<typename Dest, typename Workspace>
349 inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const
350 {
351 workspace.resize(dst.rows());
352 for(Index k = 0; k < m_length; ++k)
353 {
354 Index actual_k = m_reverse ? m_length-k-1 : k;
355 dst.rightCols(rows()-m_shift-actual_k)
356 .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
357 }
358 }
359
361 template<typename Dest> inline void applyThisOnTheLeft(Dest& dst, bool inputIsIdentity = false) const
362 {
363 Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace;
364 applyThisOnTheLeft(dst, workspace, inputIsIdentity);
365 }
366
368 template<typename Dest, typename Workspace>
369 inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace, bool inputIsIdentity = false) const
370 {
371 if(inputIsIdentity && m_reverse)
372 inputIsIdentity = false;
373 // if the entries are large enough, then apply the reflectors by block
374 if(m_length>=BlockSize && dst.cols()>1)
375 {
376 // Make sure we have at least 2 useful blocks, otherwise it is point-less:
377 Index blockSize = m_length<Index(2*BlockSize) ? (m_length+1)/2 : Index(BlockSize);
378 for(Index i = 0; i < m_length; i+=blockSize)
379 {
380 Index end = m_reverse ? (std::min)(m_length,i+blockSize) : m_length-i;
381 Index k = m_reverse ? i : (std::max)(Index(0),end-blockSize);
382 Index bs = end-k;
383 Index start = k + m_shift;
384
385 typedef Block<typename internal::remove_all<VectorsType>::type,Dynamic,Dynamic> SubVectorsType;
386 SubVectorsType sub_vecs1(m_vectors.const_cast_derived(), Side==OnTheRight ? k : start,
387 Side==OnTheRight ? start : k,
388 Side==OnTheRight ? bs : m_vectors.rows()-start,
389 Side==OnTheRight ? m_vectors.cols()-start : bs);
390 typename internal::conditional<Side==OnTheRight, Transpose<SubVectorsType>, SubVectorsType&>::type sub_vecs(sub_vecs1);
391
392 Index dstStart = dst.rows()-rows()+m_shift+k;
393 Index dstRows = rows()-m_shift-k;
394 Block<Dest,Dynamic,Dynamic> sub_dst(dst,
395 dstStart,
396 inputIsIdentity ? dstStart : 0,
397 dstRows,
398 inputIsIdentity ? dstRows : dst.cols());
399 apply_block_householder_on_the_left(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_reverse);
400 }
401 }
402 else
403 {
404 workspace.resize(dst.cols());
405 for(Index k = 0; k < m_length; ++k)
406 {
407 Index actual_k = m_reverse ? k : m_length-k-1;
408 Index dstStart = rows()-m_shift-actual_k;
409 dst.bottomRightCorner(dstStart, inputIsIdentity ? dstStart : dst.cols())
410 .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
411 }
412 }
413 }
414
422 template<typename OtherDerived>
424 {
426 res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
427 applyThisOnTheLeft(res, internal::is_identity<OtherDerived>::value && res.rows()==res.cols());
428 return res;
429 }
430
431 template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
432
442 EIGEN_DEVICE_FUNC
444 {
445 m_length = length;
446 return *this;
447 }
448
460 EIGEN_DEVICE_FUNC
462 {
463 m_shift = shift;
464 return *this;
465 }
466
467 EIGEN_DEVICE_FUNC
468 Index length() const { return m_length; }
470 EIGEN_DEVICE_FUNC
471 Index shift() const { return m_shift; }
473 /* Necessary for .adjoint() and .conjugate() */
474 template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence;
475
476 protected:
477
488 HouseholderSequence& setReverseFlag(bool reverse)
489 {
490 m_reverse = reverse;
491 return *this;
492 }
493
494 bool reverseFlag() const { return m_reverse; }
496 typename VectorsType::Nested m_vectors;
497 typename CoeffsType::Nested m_coeffs;
498 bool m_reverse;
499 Index m_length;
500 Index m_shift;
501 enum { BlockSize = 48 };
502};
503
512template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
514{
516 res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
517 h.applyThisOnTheRight(res);
518 return res;
519}
520
525template<typename VectorsType, typename CoeffsType>
526HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
527{
529}
530
537template<typename VectorsType, typename CoeffsType>
539{
541}
542
543} // end namespace Eigen
544
545#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:105
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:65
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:121
AdjointReturnType inverse() const
Inverse of the Householder sequence (equals the adjoint).
Definition: HouseholderSequence.h:275
HouseholderSequence & setLength(Index length)
Sets the length of the Householder sequence.
Definition: HouseholderSequence.h:443
Index shift() const
Returns the shift of the Householder sequence.
Definition: HouseholderSequence.h:471
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Number of columns of transformation viewed as a matrix.
Definition: HouseholderSequence.h:212
HouseholderSequence & setShift(Index shift)
Sets the shift of the Householder sequence.
Definition: HouseholderSequence.h:461
internal::conditional< Cond, ConjugateReturnType, ConstHouseholderSequence >::type conjugateIf() const
Definition: HouseholderSequence.h:259
HouseholderSequence(const HouseholderSequence &other)
Copy constructor.
Definition: HouseholderSequence.h:191
internal::matrix_type_times_scalar_type< Scalar, OtherDerived >::Type operator*(const MatrixBase< OtherDerived > &other) const
Computes the product of a Householder sequence with a matrix.
Definition: HouseholderSequence.h:423
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Number of rows of transformation viewed as a matrix.
Definition: HouseholderSequence.h:205
Index length() const
Returns the length of the Householder sequence.
Definition: HouseholderSequence.h:468
ConjugateReturnType conjugate() const
Complex conjugate of the Householder sequence.
Definition: HouseholderSequence.h:245
const EssentialVectorType essentialVector(Index k) const
Essential part of a Householder vector.
Definition: HouseholderSequence.h:229
TransposeReturnType transpose() const
Transpose of the Householder sequence.
Definition: HouseholderSequence.h:236
AdjointReturnType adjoint() const
Adjoint (conjugate transpose) of the Householder sequence.
Definition: HouseholderSequence.h:266
HouseholderSequence(const VectorsType &v, const CoeffsType &h)
Constructor.
Definition: HouseholderSequence.h:183
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:526
HouseholderSequence< VectorsType, CoeffsType, OnTheRight > rightHouseholderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:538
@ ColMajor
Definition: Constants.h:319
@ AutoAlign
Definition: Constants.h:323
@ OnTheLeft
Definition: Constants.h:332
@ OnTheRight
Definition: Constants.h:334
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
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39