HouseholderSequence.h
Go to the documentation of this file.
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 {
68  ColsAtCompileTime = RowsAtCompileTime,
71  MaxColsAtCompileTime = MaxRowsAtCompileTime,
72  Flags = 0
73  };
74 };
75 
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 {
83 };
84 
85 template<typename VectorsType, typename CoeffsType, int Side>
87 {
90  static 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 {
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 {
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 {
123 
124  public:
125  enum {
130  };
132 
133  typedef HouseholderSequence<
136  VectorsType>::type,
137  typename internal::conditional<NumTraits<Scalar>::IsComplex,
139  CoeffsType>::type,
140  Side
142 
160  HouseholderSequence(const VectorsType& v, const CoeffsType& h)
161  : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()),
162  m_shift(0)
163  {
164  }
165 
168  : m_vectors(other.m_vectors),
169  m_coeffs(other.m_coeffs),
170  m_trans(other.m_trans),
171  m_length(other.m_length),
172  m_shift(other.m_shift)
173  {
174  }
175 
180  Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
181 
186  Index cols() const { return rows(); }
187 
202  const EssentialVectorType essentialVector(Index k) const
203  {
204  eigen_assert(k >= 0 && k < m_length);
206  }
207 
210  {
211  return HouseholderSequence(*this).setTrans(!m_trans);
212  }
213 
216  {
217  return ConjugateReturnType(m_vectors.conjugate(), m_coeffs.conjugate())
218  .setTrans(m_trans)
219  .setLength(m_length)
220  .setShift(m_shift);
221  }
222 
225  {
226  return conjugate().setTrans(!m_trans);
227  }
228 
230  ConjugateReturnType inverse() const { return adjoint(); }
231 
233  template<typename DestType> inline void evalTo(DestType& dst) const
234  {
235  Matrix<Scalar, DestType::RowsAtCompileTime, 1,
236  AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows());
237  evalTo(dst, workspace);
238  }
239 
241  template<typename Dest, typename Workspace>
242  void evalTo(Dest& dst, Workspace& workspace) const
243  {
244  workspace.resize(rows());
245  Index vecs = m_length;
246  if(internal::is_same_dense(dst,m_vectors))
247  {
248  // in-place
249  dst.diagonal().setOnes();
250  dst.template triangularView<StrictlyUpper>().setZero();
251  for(Index k = vecs-1; k >= 0; --k)
252  {
253  Index cornerSize = rows() - k - m_shift;
254  if(m_trans)
255  dst.bottomRightCorner(cornerSize, cornerSize)
256  .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
257  else
258  dst.bottomRightCorner(cornerSize, cornerSize)
259  .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
260 
261  // clear the off diagonal vector
262  dst.col(k).tail(rows()-k-1).setZero();
263  }
264  // clear the remaining columns if needed
265  for(Index k = 0; k<cols()-vecs ; ++k)
266  dst.col(k).tail(rows()-k-1).setZero();
267  }
268  else
269  {
270  dst.setIdentity(rows(), rows());
271  for(Index k = vecs-1; k >= 0; --k)
272  {
273  Index cornerSize = rows() - k - m_shift;
274  if(m_trans)
275  dst.bottomRightCorner(cornerSize, cornerSize)
276  .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
277  else
278  dst.bottomRightCorner(cornerSize, cornerSize)
279  .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
280  }
281  }
282  }
283 
285  template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
286  {
288  applyThisOnTheRight(dst, workspace);
289  }
290 
292  template<typename Dest, typename Workspace>
293  inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const
294  {
295  workspace.resize(dst.rows());
296  for(Index k = 0; k < m_length; ++k)
297  {
298  Index actual_k = m_trans ? m_length-k-1 : k;
299  dst.rightCols(rows()-m_shift-actual_k)
300  .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
301  }
302  }
303 
305  template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
306  {
308  applyThisOnTheLeft(dst, workspace);
309  }
310 
312  template<typename Dest, typename Workspace>
313  inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const
314  {
315  const Index BlockSize = 48;
316  // if the entries are large enough, then apply the reflectors by block
317  if(m_length>=BlockSize && dst.cols()>1)
318  {
319  for(Index i = 0; i < m_length; i+=BlockSize)
320  {
321  Index end = m_trans ? (std::min)(m_length,i+BlockSize) : m_length-i;
322  Index k = m_trans ? i : (std::max)(Index(0),end-BlockSize);
323  Index bs = end-k;
324  Index start = k + m_shift;
325 
326  typedef Block<typename internal::remove_all<VectorsType>::type,Dynamic,Dynamic> SubVectorsType;
327  SubVectorsType sub_vecs1(m_vectors.const_cast_derived(), Side==OnTheRight ? k : start,
328  Side==OnTheRight ? start : k,
329  Side==OnTheRight ? bs : m_vectors.rows()-start,
330  Side==OnTheRight ? m_vectors.cols()-start : bs);
331  typename internal::conditional<Side==OnTheRight, Transpose<SubVectorsType>, SubVectorsType&>::type sub_vecs(sub_vecs1);
332  Block<Dest,Dynamic,Dynamic> sub_dst(dst,dst.rows()-rows()+m_shift+k,0, rows()-m_shift-k,dst.cols());
333  apply_block_householder_on_the_left(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_trans);
334  }
335  }
336  else
337  {
338  workspace.resize(dst.cols());
339  for(Index k = 0; k < m_length; ++k)
340  {
341  Index actual_k = m_trans ? k : m_length-k-1;
342  dst.bottomRows(rows()-m_shift-actual_k)
343  .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
344  }
345  }
346  }
347 
355  template<typename OtherDerived>
357  {
360  applyThisOnTheLeft(res);
361  return res;
362  }
363 
364  template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
365 
376  {
377  m_length = length;
378  return *this;
379  }
380 
393  {
394  m_shift = shift;
395  return *this;
396  }
397 
398  Index length() const { return m_length; }
399  Index shift() const { return m_shift; }
401  /* Necessary for .adjoint() and .conjugate() */
402  template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence;
403 
404  protected:
405 
415  {
416  m_trans = trans;
417  return *this;
418  }
419 
420  bool trans() const { return m_trans; }
422  typename VectorsType::Nested m_vectors;
423  typename CoeffsType::Nested m_coeffs;
424  bool m_trans;
427 };
428 
437 template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
439 {
442  h.applyThisOnTheRight(res);
443  return res;
444 }
445 
450 template<typename VectorsType, typename CoeffsType>
452 {
454 }
455 
462 template<typename VectorsType, typename CoeffsType>
464 {
466 }
467 
468 } // end namespace Eigen
469 
470 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
HouseholderSequence< VectorsType, CoeffsType, OnTheRight > HouseholderSequenceType
SCALAR Scalar
Definition: bench_gemm.cpp:33
#define max(a, b)
Definition: datatypes.h:20
Index length() const
Returns the length of the Householder sequence.
void adjoint(const MatrixType &m)
Definition: adjoint.cpp:67
void apply_block_householder_on_the_left(MatrixType &mat, const VectorsType &vectors, const CoeffsType &hCoeffs, bool forward)
HouseholderSequence< typename internal::conditional< NumTraits< Scalar >::IsComplex, typename internal::remove_all< typename VectorsType::ConjugateReturnType >::type, VectorsType >::type, typename internal::conditional< NumTraits< Scalar >::IsComplex, typename internal::remove_all< typename CoeffsType::ConjugateReturnType >::type, CoeffsType >::type, Side > ConjugateReturnType
bool is_same_dense(const T1 &mat1, const T2 &mat2, typename enable_if< has_direct_access< T1 >::ret &&has_direct_access< T2 >::ret, T1 >::type *=0)
Definition: XprHelper.h:661
Matrix< ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime, 0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime > Type
#define min(a, b)
Definition: datatypes.h:19
HouseholderSequence & setLength(Index length)
Sets the length of the Householder sequence.
ScalarBinaryOpTraits< OtherScalarType, typename MatrixType::Scalar >::ReturnType ResultScalar
internal::hseq_side_dependent_impl< VectorsType, CoeffsType, Side >::EssentialVectorType EssentialVectorType
Index cols() const
Number of columns of transformation viewed as a matrix.
Expression of the transpose of a matrix.
Definition: Transpose.h:52
ArrayXcf v
Definition: Cwise_arg.cpp:1
HouseholderSequence(const HouseholderSequence &other)
Copy constructor.
void applyThisOnTheLeft(Dest &dst) const
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
HouseholderSequence(const VectorsType &v, const CoeffsType &h)
Constructor.
static char trans
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
HouseholderSequence transpose() const
Transpose of the Householder sequence.
Block< const VectorsType, Dynamic, 1 > EssentialVectorType
void applyThisOnTheRight(Dest &dst, Workspace &workspace) const
Sequence of Householder reflections acting on subspaces with decreasing size.
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
static const EssentialVectorType essentialVector(const HouseholderSequenceType &h, Index k)
HouseholderSequence & setTrans(bool trans)
Sets the transpose flag.
EIGEN_DEVICE_FUNC ConjugateReturnType conjugate() const
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.
ConjugateReturnType conjugate() const
Complex conjugate of the Householder sequence.
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
void applyThisOnTheLeft(Dest &dst, Workspace &workspace) const
#define eigen_assert(x)
Definition: Macros.h:579
EIGEN_DEVICE_FUNC NewType cast(const OldType &x)
void evalTo(Dest &dst, Workspace &workspace) const
HouseholderSequence< VectorsType, CoeffsType, OnTheRight > rightHouseholderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Index rows() const
Number of rows of transformation viewed as a matrix.
HouseholderSequence & setShift(Index shift)
Sets the shift of the Householder sequence.
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorUInt128< uint64_t, uint64_t > operator*(const TensorUInt128< HL, LL > &lhs, const TensorUInt128< HR, LR > &rhs)
void applyThisOnTheRight(Dest &dst) const
ConjugateReturnType adjoint() const
Adjoint (conjugate transpose) of the Householder sequence.
void evalTo(DestType &dst) const
const double h
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:103
ConjugateReturnType inverse() const
Inverse of the Householder sequence (equals the adjoint).
const EssentialVectorType essentialVector(Index k) const
Essential part of a Householder vector.
bool trans() const
Returns the transpose flag.
Determines whether the given binary operation of two numeric types is allowed and what the scalar ret...
Definition: XprHelper.h:766
internal::conditional< NumTraits< Scalar >::IsComplex, const CwiseUnaryOp< internal::scalar_conjugate_op< Scalar >, const Derived >, const Derived & >::type ConjugateReturnType
const int Dynamic
Definition: Constants.h:21
Index shift() const
Returns the shift of the Householder sequence.
The matrix class, also used for vectors and row-vectors.
HouseholderSequence< VectorsType, CoeffsType, OnTheLeft > HouseholderSequenceType
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
static const EssentialVectorType essentialVector(const HouseholderSequenceType &h, Index k)
v setZero(3)
internal::traits< HouseholderSequence >::Scalar Scalar


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:42:12