Eigen  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
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
14#include "./InternalHeaderCheck.h"
15
16namespace Eigen {
17
59namespace internal {
60
61template<typename VectorsType, typename CoeffsType, int Side>
62struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
63{
64 typedef typename VectorsType::Scalar Scalar;
65 typedef typename VectorsType::StorageIndex StorageIndex;
66 typedef typename VectorsType::StorageKind StorageKind;
67 enum {
68 RowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::RowsAtCompileTime
69 : traits<VectorsType>::ColsAtCompileTime,
70 ColsAtCompileTime = RowsAtCompileTime,
71 MaxRowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::MaxRowsAtCompileTime
72 : traits<VectorsType>::MaxColsAtCompileTime,
73 MaxColsAtCompileTime = MaxRowsAtCompileTime,
74 Flags = 0
75 };
76};
77
78struct HouseholderSequenceShape {};
79
80template<typename VectorsType, typename CoeffsType, int Side>
81struct evaluator_traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
82 : public evaluator_traits_base<HouseholderSequence<VectorsType,CoeffsType,Side> >
83{
84 typedef HouseholderSequenceShape Shape;
85};
86
87template<typename VectorsType, typename CoeffsType, int Side>
88struct hseq_side_dependent_impl
89{
90 typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType;
91 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType;
92 static EIGEN_DEVICE_FUNC inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
93 {
94 Index start = k+1+h.m_shift;
95 return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
96 }
97};
98
99template<typename VectorsType, typename CoeffsType>
100struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight>
101{
102 typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType;
103 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType;
104 static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
105 {
106 Index start = k+1+h.m_shift;
107 return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose();
108 }
109};
110
111template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type
112{
113 typedef typename ScalarBinaryOpTraits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
114 ResultScalar;
115 typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
116 0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
117};
118
119} // end namespace internal
120
121template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence
122 : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
123{
125
126 public:
127 enum {
128 RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
129 ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
130 MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
131 MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
132 };
133 typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
134
135 typedef HouseholderSequence<
136 typename internal::conditional<NumTraits<Scalar>::IsComplex,
137 typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type,
138 VectorsType>::type,
139 typename internal::conditional<NumTraits<Scalar>::IsComplex,
140 typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
141 CoeffsType>::type,
142 Side
144
145 typedef HouseholderSequence<
146 VectorsType,
147 typename internal::conditional<NumTraits<Scalar>::IsComplex,
148 typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
149 CoeffsType>::type,
150 Side
152
153 typedef HouseholderSequence<
154 typename internal::conditional<NumTraits<Scalar>::IsComplex,
155 typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type,
156 VectorsType>::type,
157 CoeffsType,
158 Side
160
161 typedef HouseholderSequence<
164 Side
166
184 EIGEN_DEVICE_FUNC
185 HouseholderSequence(const VectorsType& v, const CoeffsType& h)
186 : m_vectors(v), m_coeffs(h), m_reverse(false), m_length(v.diagonalSize()),
187 m_shift(0)
188 {
189 }
190
192 EIGEN_DEVICE_FUNC
194 : m_vectors(other.m_vectors),
195 m_coeffs(other.m_coeffs),
196 m_reverse(other.m_reverse),
197 m_length(other.m_length),
198 m_shift(other.m_shift)
199 {
200 }
201
206 EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
207 Index rows() const EIGEN_NOEXCEPT { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
208
213 EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
214 Index cols() const EIGEN_NOEXCEPT { return rows(); }
215
230 EIGEN_DEVICE_FUNC
232 {
233 eigen_assert(k >= 0 && k < m_length);
234 return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
235 }
236
239 {
240 return TransposeReturnType(m_vectors.conjugate(), m_coeffs)
241 .setReverseFlag(!m_reverse)
242 .setLength(m_length)
243 .setShift(m_shift);
244 }
245
248 {
249 return ConjugateReturnType(m_vectors.conjugate(), m_coeffs.conjugate())
250 .setReverseFlag(m_reverse)
251 .setLength(m_length)
252 .setShift(m_shift);
253 }
254
258 template<bool Cond>
259 EIGEN_DEVICE_FUNC
260 inline typename internal::conditional<Cond,ConjugateReturnType,ConstHouseholderSequence>::type
262 {
263 typedef typename internal::conditional<Cond,ConjugateReturnType,ConstHouseholderSequence>::type ReturnType;
264 return ReturnType(m_vectors.template conjugateIf<Cond>(), m_coeffs.template conjugateIf<Cond>());
265 }
266
269 {
270 return AdjointReturnType(m_vectors, m_coeffs.conjugate())
271 .setReverseFlag(!m_reverse)
272 .setLength(m_length)
273 .setShift(m_shift);
274 }
275
277 AdjointReturnType inverse() const { return adjoint(); }
278
280 template<typename DestType>
281 inline EIGEN_DEVICE_FUNC
282 void evalTo(DestType& dst) const
283 {
284 Matrix<Scalar, DestType::RowsAtCompileTime, 1,
285 AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows());
286 evalTo(dst, workspace);
287 }
288
290 template<typename Dest, typename Workspace>
291 EIGEN_DEVICE_FUNC
292 void evalTo(Dest& dst, Workspace& workspace) const
293 {
294 workspace.resize(rows());
295 Index vecs = m_length;
296 if(internal::is_same_dense(dst,m_vectors))
297 {
298 // in-place
299 dst.diagonal().setOnes();
300 dst.template triangularView<StrictlyUpper>().setZero();
301 for(Index k = vecs-1; k >= 0; --k)
302 {
303 Index cornerSize = rows() - k - m_shift;
304 if(m_reverse)
305 dst.bottomRightCorner(cornerSize, cornerSize)
306 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
307 else
308 dst.bottomRightCorner(cornerSize, cornerSize)
309 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
310
311 // clear the off diagonal vector
312 dst.col(k).tail(rows()-k-1).setZero();
313 }
314 // clear the remaining columns if needed
315 for(Index k = 0; k<cols()-vecs ; ++k)
316 dst.col(k).tail(rows()-k-1).setZero();
317 }
318 else if(m_length>BlockSize)
319 {
320 dst.setIdentity(rows(), rows());
321 if(m_reverse)
322 applyThisOnTheLeft(dst,workspace,true);
323 else
324 applyThisOnTheLeft(dst,workspace,true);
325 }
326 else
327 {
328 dst.setIdentity(rows(), rows());
329 for(Index k = vecs-1; k >= 0; --k)
330 {
331 Index cornerSize = rows() - k - m_shift;
332 if(m_reverse)
333 dst.bottomRightCorner(cornerSize, cornerSize)
334 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
335 else
336 dst.bottomRightCorner(cornerSize, cornerSize)
337 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
338 }
339 }
340 }
341
343 template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
344 {
345 Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
346 applyThisOnTheRight(dst, workspace);
347 }
348
350 template<typename Dest, typename Workspace>
351 inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const
352 {
353 workspace.resize(dst.rows());
354 for(Index k = 0; k < m_length; ++k)
355 {
356 Index actual_k = m_reverse ? m_length-k-1 : k;
357 dst.rightCols(rows()-m_shift-actual_k)
358 .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
359 }
360 }
361
363 template<typename Dest> inline void applyThisOnTheLeft(Dest& dst, bool inputIsIdentity = false) const
364 {
365 Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace;
366 applyThisOnTheLeft(dst, workspace, inputIsIdentity);
367 }
368
370 template<typename Dest, typename Workspace>
371 inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace, bool inputIsIdentity = false) const
372 {
373 if(inputIsIdentity && m_reverse)
374 inputIsIdentity = false;
375 // if the entries are large enough, then apply the reflectors by block
376 if(m_length>=BlockSize && dst.cols()>1)
377 {
378 // Make sure we have at least 2 useful blocks, otherwise it is point-less:
379 Index blockSize = m_length<Index(2*BlockSize) ? (m_length+1)/2 : Index(BlockSize);
380 for(Index i = 0; i < m_length; i+=blockSize)
381 {
382 Index end = m_reverse ? (std::min)(m_length,i+blockSize) : m_length-i;
383 Index k = m_reverse ? i : (std::max)(Index(0),end-blockSize);
384 Index bs = end-k;
385 Index start = k + m_shift;
386
387 typedef Block<typename internal::remove_all<VectorsType>::type,Dynamic,Dynamic> SubVectorsType;
388 SubVectorsType sub_vecs1(m_vectors.const_cast_derived(), Side==OnTheRight ? k : start,
389 Side==OnTheRight ? start : k,
390 Side==OnTheRight ? bs : m_vectors.rows()-start,
391 Side==OnTheRight ? m_vectors.cols()-start : bs);
392 typename internal::conditional<Side==OnTheRight, Transpose<SubVectorsType>, SubVectorsType&>::type sub_vecs(sub_vecs1);
393
394 Index dstStart = dst.rows()-rows()+m_shift+k;
395 Index dstRows = rows()-m_shift-k;
396 Block<Dest,Dynamic,Dynamic> sub_dst(dst,
397 dstStart,
398 inputIsIdentity ? dstStart : 0,
399 dstRows,
400 inputIsIdentity ? dstRows : dst.cols());
401 apply_block_householder_on_the_left(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_reverse);
402 }
403 }
404 else
405 {
406 workspace.resize(dst.cols());
407 for(Index k = 0; k < m_length; ++k)
408 {
409 Index actual_k = m_reverse ? k : m_length-k-1;
410 Index dstStart = rows()-m_shift-actual_k;
411 dst.bottomRightCorner(dstStart, inputIsIdentity ? dstStart : dst.cols())
412 .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
413 }
414 }
415 }
416
424 template<typename OtherDerived>
426 {
428 res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
429 applyThisOnTheLeft(res, internal::is_identity<OtherDerived>::value && res.rows()==res.cols());
430 return res;
431 }
432
433 template<typename VectorsType_, typename CoeffsType_, int Side_> friend struct internal::hseq_side_dependent_impl;
434
444 EIGEN_DEVICE_FUNC
446 {
447 m_length = length;
448 return *this;
449 }
450
462 EIGEN_DEVICE_FUNC
464 {
465 m_shift = shift;
466 return *this;
467 }
468
469 EIGEN_DEVICE_FUNC
470 Index length() const { return m_length; }
472 EIGEN_DEVICE_FUNC
473 Index shift() const { return m_shift; }
475 /* Necessary for .adjoint() and .conjugate() */
476 template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence;
477
478 protected:
479
490 HouseholderSequence& setReverseFlag(bool reverse)
491 {
492 m_reverse = reverse;
493 return *this;
494 }
495
496 bool reverseFlag() const { return m_reverse; }
498 typename VectorsType::Nested m_vectors;
499 typename CoeffsType::Nested m_coeffs;
500 bool m_reverse;
501 Index m_length;
502 Index m_shift;
503 enum { BlockSize = 48 };
504};
505
514template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
516{
518 res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
519 h.applyThisOnTheRight(res);
520 return res;
521}
522
527template<typename VectorsType, typename CoeffsType>
528HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
529{
531}
532
539template<typename VectorsType, typename CoeffsType>
541{
543}
544
545} // end namespace Eigen
546
547#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:107
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:67
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:123
AdjointReturnType inverse() const
Inverse of the Householder sequence (equals the adjoint).
Definition: HouseholderSequence.h:277
HouseholderSequence & setLength(Index length)
Sets the length of the Householder sequence.
Definition: HouseholderSequence.h:445
Index shift() const
Returns the shift of the Householder sequence.
Definition: HouseholderSequence.h:473
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Number of columns of transformation viewed as a matrix.
Definition: HouseholderSequence.h:214
HouseholderSequence & setShift(Index shift)
Sets the shift of the Householder sequence.
Definition: HouseholderSequence.h:463
internal::conditional< Cond, ConjugateReturnType, ConstHouseholderSequence >::type conjugateIf() const
Definition: HouseholderSequence.h:261
HouseholderSequence(const HouseholderSequence &other)
Copy constructor.
Definition: HouseholderSequence.h:193
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:425
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Number of rows of transformation viewed as a matrix.
Definition: HouseholderSequence.h:207
Index length() const
Returns the length of the Householder sequence.
Definition: HouseholderSequence.h:470
ConjugateReturnType conjugate() const
Complex conjugate of the Householder sequence.
Definition: HouseholderSequence.h:247
const EssentialVectorType essentialVector(Index k) const
Essential part of a Householder vector.
Definition: HouseholderSequence.h:231
TransposeReturnType transpose() const
Transpose of the Householder sequence.
Definition: HouseholderSequence.h:238
AdjointReturnType adjoint() const
Adjoint (conjugate transpose) of the Householder sequence.
Definition: HouseholderSequence.h:268
HouseholderSequence(const VectorsType &v, const CoeffsType &h)
Constructor.
Definition: HouseholderSequence.h:185
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:182
static const lastp1_t end
Definition: IndexedViewHelper.h:183
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:528
HouseholderSequence< VectorsType, CoeffsType, OnTheRight > rightHouseholderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:540
@ ColMajor
Definition: Constants.h:321
@ AutoAlign
Definition: Constants.h:325
@ OnTheLeft
Definition: Constants.h:334
@ OnTheRight
Definition: Constants.h:336
Namespace containing all symbols from the Eigen library.
Definition: B01_Experimental.dox:1
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Product< MatrixDerived, PermutationDerived, AliasFreeProduct > operator*(const MatrixBase< MatrixDerived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:517
const int Dynamic
Definition: Constants.h:24
Definition: EigenBase.h:32
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41