Eigen  3.2.92
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 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 scalar_product_traits<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 
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);
205  return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
206  }
207 
210  {
211  return HouseholderSequence(*this).setTrans(!m_trans);
212  }
213 
215  ConjugateReturnType conjugate() const
216  {
217  return ConjugateReturnType(m_vectors.conjugate(), m_coeffs.conjugate())
218  .setTrans(m_trans)
219  .setLength(m_length)
220  .setShift(m_shift);
221  }
222 
224  ConjugateReturnType adjoint() const
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<typename internal::remove_all<VectorsType>::type,Dest>::value
247  && internal::extract_data(dst) == internal::extract_data(m_vectors))
248  {
249  // in-place
250  dst.diagonal().setOnes();
251  dst.template triangularView<StrictlyUpper>().setZero();
252  for(Index k = vecs-1; k >= 0; --k)
253  {
254  Index cornerSize = rows() - k - m_shift;
255  if(m_trans)
256  dst.bottomRightCorner(cornerSize, cornerSize)
257  .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
258  else
259  dst.bottomRightCorner(cornerSize, cornerSize)
260  .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
261 
262  // clear the off diagonal vector
263  dst.col(k).tail(rows()-k-1).setZero();
264  }
265  // clear the remaining columns if needed
266  for(Index k = 0; k<cols()-vecs ; ++k)
267  dst.col(k).tail(rows()-k-1).setZero();
268  }
269  else
270  {
271  dst.setIdentity(rows(), rows());
272  for(Index k = vecs-1; k >= 0; --k)
273  {
274  Index cornerSize = rows() - k - m_shift;
275  if(m_trans)
276  dst.bottomRightCorner(cornerSize, cornerSize)
277  .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
278  else
279  dst.bottomRightCorner(cornerSize, cornerSize)
280  .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
281  }
282  }
283  }
284 
286  template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
287  {
289  applyThisOnTheRight(dst, workspace);
290  }
291 
293  template<typename Dest, typename Workspace>
294  inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const
295  {
296  workspace.resize(dst.rows());
297  for(Index k = 0; k < m_length; ++k)
298  {
299  Index actual_k = m_trans ? m_length-k-1 : k;
300  dst.rightCols(rows()-m_shift-actual_k)
301  .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
302  }
303  }
304 
306  template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
307  {
309  applyThisOnTheLeft(dst, workspace);
310  }
311 
313  template<typename Dest, typename Workspace>
314  inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const
315  {
316  const Index BlockSize = 48;
317  // if the entries are large enough, then apply the reflectors by block
318  if(m_length>=BlockSize && dst.cols()>1)
319  {
320  for(Index i = 0; i < m_length; i+=BlockSize)
321  {
322  Index end = m_trans ? (std::min)(m_length,i+BlockSize) : m_length-i;
323  Index k = m_trans ? i : (std::max)(Index(0),end-BlockSize);
324  Index bs = end-k;
325  Index start = k + m_shift;
326 
327  typedef Block<typename internal::remove_all<VectorsType>::type,Dynamic,Dynamic> SubVectorsType;
328  SubVectorsType sub_vecs1(m_vectors.const_cast_derived(), Side==OnTheRight ? k : start,
329  Side==OnTheRight ? start : k,
330  Side==OnTheRight ? bs : m_vectors.rows()-start,
331  Side==OnTheRight ? m_vectors.cols()-start : bs);
332  typename internal::conditional<Side==OnTheRight, Transpose<SubVectorsType>, SubVectorsType&>::type sub_vecs(sub_vecs1);
333  Block<Dest,Dynamic,Dynamic> sub_dst(dst,dst.rows()-rows()+m_shift+k,0, rows()-m_shift-k,dst.cols());
334  apply_block_householder_on_the_left(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_trans);
335  }
336  }
337  else
338  {
339  workspace.resize(dst.cols());
340  for(Index k = 0; k < m_length; ++k)
341  {
342  Index actual_k = m_trans ? k : m_length-k-1;
343  dst.bottomRows(rows()-m_shift-actual_k)
344  .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
345  }
346  }
347  }
348 
356  template<typename OtherDerived>
358  {
360  res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
361  applyThisOnTheLeft(res);
362  return res;
363  }
364 
365  template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
366 
377  {
378  m_length = length;
379  return *this;
380  }
381 
394  {
395  m_shift = shift;
396  return *this;
397  }
398 
399  Index length() const { return m_length; }
400  Index shift() const { return m_shift; }
402  /* Necessary for .adjoint() and .conjugate() */
403  template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence;
404 
405  protected:
406 
416  {
417  m_trans = trans;
418  return *this;
419  }
420 
421  bool trans() const { return m_trans; }
423  typename VectorsType::Nested m_vectors;
424  typename CoeffsType::Nested m_coeffs;
425  bool m_trans;
426  Index m_length;
427  Index m_shift;
428 };
429 
438 template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
440 {
442  res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
443  h.applyThisOnTheRight(res);
444  return res;
445 }
446 
451 template<typename VectorsType, typename CoeffsType>
452 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
453 {
455 }
456 
463 template<typename VectorsType, typename CoeffsType>
465 {
467 }
468 
469 } // end namespace Eigen
470 
471 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
Index length() const
Returns the length of the Householder sequence.
Definition: HouseholderSequence.h:399
Definition: Constants.h:333
HouseholderSequence & setLength(Index length)
Sets the length of the Householder sequence.
Definition: HouseholderSequence.h:376
Index cols() const
Number of columns of transformation viewed as a matrix.
Definition: HouseholderSequence.h:186
HouseholderSequence(const HouseholderSequence &other)
Copy constructor.
Definition: HouseholderSequence.h:167
Definition: LDLT.h:16
HouseholderSequence(const VectorsType &v, const CoeffsType &h)
Constructor.
Definition: HouseholderSequence.h:160
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:37
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:452
HouseholderSequence transpose() const
Transpose of the Householder sequence.
Definition: HouseholderSequence.h:209
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: ForwardDeclarations.h:259
Definition: Constants.h:320
HouseholderSequence & setTrans(bool trans)
Sets the transpose flag.
Definition: HouseholderSequence.h:415
Definition: Constants.h:335
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:357
ConjugateReturnType conjugate() const
Complex conjugate of the Householder sequence.
Definition: HouseholderSequence.h:215
HouseholderSequence< VectorsType, CoeffsType, OnTheRight > rightHouseholderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:464
Index rows() const
Number of rows of transformation viewed as a matrix.
Definition: HouseholderSequence.h:180
HouseholderSequence & setShift(Index shift)
Sets the shift of the Householder sequence.
Definition: HouseholderSequence.h:393
ConjugateReturnType adjoint() const
Adjoint (conjugate transpose) of the Householder sequence.
Definition: HouseholderSequence.h:224
Definition: Eigen_Colamd.h:54
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:104
ConjugateReturnType inverse() const
Inverse of the Householder sequence (equals the adjoint).
Definition: HouseholderSequence.h:230
const EssentialVectorType essentialVector(Index k) const
Essential part of a Householder vector.
Definition: HouseholderSequence.h:202
Definition: Constants.h:324
bool trans() const
Returns the transpose flag.
Definition: HouseholderSequence.h:421
Index shift() const
Returns the shift of the Householder sequence.
Definition: HouseholderSequence.h:400
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48