Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
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 
16 namespace Eigen {
17 
59 namespace internal {
60 
61 template<typename VectorsType, typename CoeffsType, int Side>
62 struct 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 
78 struct HouseholderSequenceShape {};
79 
80 template<typename VectorsType, typename CoeffsType, int Side>
81 struct evaluator_traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
82  : public evaluator_traits_base<HouseholderSequence<VectorsType,CoeffsType,Side> >
83 {
84  typedef HouseholderSequenceShape Shape;
85 };
86 
87 template<typename VectorsType, typename CoeffsType, int Side>
88 struct 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 
99 template<typename VectorsType, typename CoeffsType>
100 struct 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 
111 template<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 
121 template<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  std::conditional_t<NumTraits<Scalar>::IsComplex,
137  internal::remove_all_t<typename VectorsType::ConjugateReturnType>,
138  VectorsType>,
139  std::conditional_t<NumTraits<Scalar>::IsComplex,
140  internal::remove_all_t<typename CoeffsType::ConjugateReturnType>,
141  CoeffsType>,
142  Side
144 
145  typedef HouseholderSequence<
146  VectorsType,
147  std::conditional_t<NumTraits<Scalar>::IsComplex,
148  internal::remove_all_t<typename CoeffsType::ConjugateReturnType>,
149  CoeffsType>,
150  Side
152 
153  typedef HouseholderSequence<
154  std::conditional_t<NumTraits<Scalar>::IsComplex,
155  internal::remove_all_t<typename VectorsType::ConjugateReturnType>,
156  VectorsType>,
157  CoeffsType,
158  Side
160 
161  typedef HouseholderSequence<
162  std::add_const_t<VectorsType>,
163  std::add_const_t<CoeffsType>,
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 std::conditional_t<Cond,ConjugateReturnType,ConstHouseholderSequence>
261  conjugateIf() const
262  {
263  typedef std::conditional_t<Cond,ConjugateReturnType,ConstHouseholderSequence> 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<internal::remove_all_t<VectorsType>,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  std::conditional_t<Side==OnTheRight, Transpose<SubVectorsType>, SubVectorsType&> sub_vecs(sub_vecs1);
393 
394  Index dstRows = rows()-m_shift-k;
395 
396  if (inputIsIdentity)
397  {
398  Block<Dest, Dynamic, Dynamic> sub_dst = dst.bottomRightCorner(dstRows, dstRows);
399  apply_block_householder_on_the_left(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_reverse);
400  }
401  else
402  {
403  auto sub_dst = dst.bottomRows(dstRows);
404  apply_block_householder_on_the_left(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_reverse);
405  }
406  }
407  }
408  else
409  {
410  workspace.resize(dst.cols());
411  for(Index k = 0; k < m_length; ++k)
412  {
413  Index actual_k = m_reverse ? k : m_length-k-1;
414  Index dstRows = rows()-m_shift-actual_k;
415 
416  if (inputIsIdentity)
417  {
418  Block<Dest, Dynamic, Dynamic> sub_dst = dst.bottomRightCorner(dstRows, dstRows);
419  sub_dst.applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
420  }
421  else
422  {
423  auto sub_dst = dst.bottomRows(dstRows);
424  sub_dst.applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
425  }
426  }
427  }
428  }
429 
437  template<typename OtherDerived>
439  {
441  res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
442  applyThisOnTheLeft(res, internal::is_identity<OtherDerived>::value && res.rows()==res.cols());
443  return res;
444  }
445 
446  template<typename VectorsType_, typename CoeffsType_, int Side_> friend struct internal::hseq_side_dependent_impl;
447 
457  EIGEN_DEVICE_FUNC
459  {
460  m_length = length;
461  return *this;
462  }
463 
475  EIGEN_DEVICE_FUNC
477  {
478  m_shift = shift;
479  return *this;
480  }
481 
482  EIGEN_DEVICE_FUNC
483  Index length() const { return m_length; }
485  EIGEN_DEVICE_FUNC
486  Index shift() const { return m_shift; }
488  /* Necessary for .adjoint() and .conjugate() */
489  template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence;
490 
491  protected:
492 
503  HouseholderSequence& setReverseFlag(bool reverse)
504  {
505  m_reverse = reverse;
506  return *this;
507  }
508 
509  bool reverseFlag() const { return m_reverse; }
511  typename VectorsType::Nested m_vectors;
512  typename CoeffsType::Nested m_coeffs;
513  bool m_reverse;
514  Index m_length;
515  Index m_shift;
516  enum { BlockSize = 48 };
517 };
518 
527 template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
529 {
531  res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
532  h.applyThisOnTheRight(res);
533  return res;
534 }
535 
540 template<typename VectorsType, typename CoeffsType>
541 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
542 {
544 }
545 
552 template<typename VectorsType, typename CoeffsType>
554 {
556 }
557 
558 } // end namespace Eigen
559 
560 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:107
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 & setShift(Index shift)
Sets the shift of the Householder sequence.
Definition: HouseholderSequence.h:476
HouseholderSequence & setLength(Index length)
Sets the length of the Householder sequence.
Definition: HouseholderSequence.h:458
Index shift() const
Returns the shift of the Householder sequence.
Definition: HouseholderSequence.h:486
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Number of columns of transformation viewed as a matrix.
Definition: HouseholderSequence.h:214
std::conditional_t< Cond, ConjugateReturnType, ConstHouseholderSequence > conjugateIf() const
Definition: HouseholderSequence.h:261
HouseholderSequence(const HouseholderSequence &other)
Copy constructor.
Definition: HouseholderSequence.h:193
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Number of rows of transformation viewed as a matrix.
Definition: HouseholderSequence.h:207
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:438
Index length() const
Returns the length of the Householder sequence.
Definition: HouseholderSequence.h:483
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:541
HouseholderSequence< VectorsType, CoeffsType, OnTheRight > rightHouseholderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:553
@ 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: Core:139
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