Path Tracer
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 namespace Eigen {
15 
57 namespace internal {
58 
59 template<typename VectorsType, typename CoeffsType, int Side>
60 struct 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 
76 struct HouseholderSequenceShape {};
77 
78 template<typename VectorsType, typename CoeffsType, int Side>
79 struct evaluator_traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
80  : public evaluator_traits_base<HouseholderSequence<VectorsType,CoeffsType,Side> >
81 {
82  typedef HouseholderSequenceShape Shape;
83 };
84 
85 template<typename VectorsType, typename CoeffsType, int Side>
86 struct 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 
97 template<typename VectorsType, typename CoeffsType>
98 struct 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 
109 template<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 
119 template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence
120  : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
121 {
122  typedef typename internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType EssentialVectorType;
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
141  > ConjugateReturnType;
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
149  > AdjointReturnType;
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
157  > TransposeReturnType;
158 
159  typedef HouseholderSequence<
160  typename internal::add_const<VectorsType>::type,
161  typename internal::add_const<CoeffsType>::type,
162  Side
163  > ConstHouseholderSequence;
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
191  HouseholderSequence(const HouseholderSequence& other)
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
205  Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
206 
211  EIGEN_DEVICE_FUNC
212  Index cols() const { return rows(); }
213 
228  EIGEN_DEVICE_FUNC
229  const EssentialVectorType essentialVector(Index k) const
230  {
231  eigen_assert(k >= 0 && k < m_length);
232  return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
233  }
234 
236  TransposeReturnType transpose() const
237  {
238  return TransposeReturnType(m_vectors.conjugate(), m_coeffs)
239  .setReverseFlag(!m_reverse)
240  .setLength(m_length)
241  .setShift(m_shift);
242  }
243 
245  ConjugateReturnType conjugate() const
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
259  conjugateIf() const
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 
266  AdjointReturnType adjoint() const
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>
423  typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const
424  {
425  typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type
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
443  HouseholderSequence& setLength(Index length)
444  {
445  m_length = length;
446  return *this;
447  }
448 
460  EIGEN_DEVICE_FUNC
461  HouseholderSequence& setShift(Index shift)
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 
512 template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
513 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence<VectorsType,CoeffsType,Side>& h)
514 {
515  typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type
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 
525 template<typename VectorsType, typename CoeffsType>
526 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
527 {
528  return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h);
529 }
530 
537 template<typename VectorsType, typename CoeffsType>
538 HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h)
539 {
540  return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h);
541 }
542 
543 } // end namespace Eigen
544 
545 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
Eigen::HouseholderSequence::HouseholderSequence
EIGEN_DEVICE_FUNC HouseholderSequence(const VectorsType &v, const CoeffsType &h)
Constructor.
Definition: HouseholderSequence.h:183
Eigen::HouseholderSequence::length
EIGEN_DEVICE_FUNC Index length() const
Returns the length of the Householder sequence.
Definition: HouseholderSequence.h:468
Eigen
Namespace containing all symbols from the Eigen library.
Definition: LDLT.h:16
Eigen::Block
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:105
Eigen::HouseholderSequence::setLength
EIGEN_DEVICE_FUNC HouseholderSequence & setLength(Index length)
Sets the length of the Householder sequence.
Definition: HouseholderSequence.h:443
Eigen::householderSequence
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:526
Eigen::operator*
EIGEN_DEVICE_FUNC const Product< MatrixDerived, PermutationDerived, AliasFreeProduct > operator*(const MatrixBase< MatrixDerived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:515
Eigen::EigenBase
Definition: EigenBase.h:30
Eigen::internal::hseq_side_dependent_impl
Definition: HouseholderSequence.h:87
Eigen::HouseholderSequence::shift
EIGEN_DEVICE_FUNC Index shift() const
Returns the shift of the Householder sequence.
Definition: HouseholderSequence.h:471
Eigen::HouseholderSequence::HouseholderSequence
EIGEN_DEVICE_FUNC HouseholderSequence(const HouseholderSequence &other)
Copy constructor.
Definition: HouseholderSequence.h:191
Eigen::internal::evaluator_traits
Definition: CoreEvaluators.h:80
Eigen::ScalarBinaryOpTraits
Determines whether the given binary operation of two numeric types is allowed and what the scalar ret...
Definition: XprHelper.h:814
Eigen::internal::HouseholderSequenceShape
Definition: HouseholderSequence.h:76
Eigen::HouseholderSequence::cols
EIGEN_DEVICE_FUNC Index cols() const
Number of columns of transformation viewed as a matrix.
Definition: HouseholderSequence.h:212
Eigen::HouseholderSequence::conjugate
ConjugateReturnType conjugate() const
Complex conjugate of the Householder sequence.
Definition: HouseholderSequence.h:245
Eigen::Transpose
Expression of the transpose of a matrix.
Definition: Transpose.h:54
Eigen::HouseholderSequence::conjugateIf
EIGEN_DEVICE_FUNC internal::conditional< Cond, ConjugateReturnType, ConstHouseholderSequence >::type conjugateIf() const
Definition: HouseholderSequence.h:259
Eigen::OnTheLeft
@ OnTheLeft
Definition: Constants.h:331
Eigen::HouseholderSequence::operator*
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::HouseholderSequence::essentialVector
EIGEN_DEVICE_FUNC const EssentialVectorType essentialVector(Index k) const
Essential part of a Householder vector.
Definition: HouseholderSequence.h:229
Eigen::internal::true_type
Definition: Meta.h:63
Eigen::internal::is_identity
Definition: XprHelper.h:687
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:21
Eigen::Diagonal
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:65
Eigen::internal::evaluator_traits_base
Definition: CoreEvaluators.h:71
Eigen::HouseholderSequence::transpose
TransposeReturnType transpose() const
Transpose of the Householder sequence.
Definition: HouseholderSequence.h:236
Eigen::AutoAlign
@ AutoAlign
Definition: Constants.h:322
Eigen::OnTheRight
@ OnTheRight
Definition: Constants.h:333
Eigen::rightHouseholderSequence
HouseholderSequence< VectorsType, CoeffsType, OnTheRight > rightHouseholderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:538
Eigen::internal::traits
Definition: ForwardDeclarations.h:17
Eigen::HouseholderSequence::adjoint
AdjointReturnType adjoint() const
Adjoint (conjugate transpose) of the Householder sequence.
Definition: HouseholderSequence.h:266
Eigen::internal::conditional
Definition: Meta.h:76
Eigen::HouseholderSequence::rows
EIGEN_DEVICE_FUNC Index rows() const
Number of rows of transformation viewed as a matrix.
Definition: HouseholderSequence.h:205
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Eigen::HouseholderSequence::inverse
AdjointReturnType inverse() const
Inverse of the Householder sequence (equals the adjoint).
Definition: HouseholderSequence.h:275
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
Eigen::HouseholderSequence::setShift
EIGEN_DEVICE_FUNC HouseholderSequence & setShift(Index shift)
Sets the shift of the Householder sequence.
Definition: HouseholderSequence.h:461
Eigen::ColMajor
@ ColMajor
Definition: Constants.h:318
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:42
Eigen::HouseholderSequence
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:121
Eigen::internal::matrix_type_times_scalar_type
Definition: HouseholderSequence.h:110
Generated by   doxygen 1.8.19