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 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25 
26 #ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H
27 #define EIGEN_HOUSEHOLDER_SEQUENCE_H
28 
29 namespace Eigen {
30 
72 namespace internal {
73 
74 template<typename VectorsType, typename CoeffsType, int Side>
75 struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
76 {
77  typedef typename VectorsType::Scalar Scalar;
78  typedef typename VectorsType::Index Index;
79  typedef typename VectorsType::StorageKind StorageKind;
80  enum {
81  RowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::RowsAtCompileTime
83  ColsAtCompileTime = RowsAtCompileTime,
84  MaxRowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::MaxRowsAtCompileTime
86  MaxColsAtCompileTime = MaxRowsAtCompileTime,
87  Flags = 0
88  };
89 };
90 
91 template<typename VectorsType, typename CoeffsType, int Side>
92 struct hseq_side_dependent_impl
93 {
94  typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType;
95  typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType;
96  typedef typename VectorsType::Index Index;
97  static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
98  {
99  Index start = k+1+h.m_shift;
100  return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
101  }
102 };
103 
104 template<typename VectorsType, typename CoeffsType>
105 struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight>
106 {
107  typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType;
108  typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType;
109  typedef typename VectorsType::Index Index;
110  static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
111  {
112  Index start = k+1+h.m_shift;
113  return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose();
114  }
115 };
116 
117 template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type
118 {
119  typedef typename scalar_product_traits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
120  ResultScalar;
121  typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
122  0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
123 };
124 
125 } // end namespace internal
126 
127 template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence
128  : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
129 {
130  enum {
131  RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
132  ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
133  MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
134  MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
135  };
136  typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
137  typedef typename VectorsType::Index Index;
138 
141 
142  public:
143 
144  typedef HouseholderSequence<
145  VectorsType,
147  typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
148  CoeffsType>::type,
149  Side
151 
169  HouseholderSequence(const VectorsType& v, const CoeffsType& h)
170  : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()),
171  m_shift(0)
172  {
173  }
174 
177  : m_vectors(other.m_vectors),
178  m_coeffs(other.m_coeffs),
179  m_trans(other.m_trans),
180  m_length(other.m_length),
181  m_shift(other.m_shift)
182  {
183  }
184 
189  Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
190 
195  Index cols() const { return rows(); }
196 
212  {
213  eigen_assert(k >= 0 && k < m_length);
215  }
216 
219  {
220  return HouseholderSequence(*this).setTrans(!m_trans);
221  }
222 
225  {
226  return ConjugateReturnType(m_vectors, m_coeffs.conjugate())
227  .setTrans(m_trans)
229  .setShift(m_shift);
230  }
231 
234  {
235  return conjugate().setTrans(!m_trans);
236  }
237 
239  ConjugateReturnType inverse() const { return adjoint(); }
240 
242  template<typename DestType> inline void evalTo(DestType& dst) const
243  {
244  Matrix<Scalar, DestType::RowsAtCompileTime, 1,
245  AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows());
246  evalTo(dst, workspace);
247  }
248 
250  template<typename Dest, typename Workspace>
251  void evalTo(Dest& dst, Workspace& workspace) const
252  {
253  workspace.resize(rows());
254  Index vecs = m_length;
255  if( internal::is_same<typename internal::remove_all<VectorsType>::type,Dest>::value
257  {
258  // in-place
259  dst.diagonal().setOnes();
260  dst.template triangularView<StrictlyUpper>().setZero();
261  for(Index k = vecs-1; k >= 0; --k)
262  {
263  Index cornerSize = rows() - k - m_shift;
264  if(m_trans)
265  dst.bottomRightCorner(cornerSize, cornerSize)
266  .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
267  else
268  dst.bottomRightCorner(cornerSize, cornerSize)
269  .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
270 
271  // clear the off diagonal vector
272  dst.col(k).tail(rows()-k-1).setZero();
273  }
274  // clear the remaining columns if needed
275  for(Index k = 0; k<cols()-vecs ; ++k)
276  dst.col(k).tail(rows()-k-1).setZero();
277  }
278  else
279  {
280  dst.setIdentity(rows(), rows());
281  for(Index k = vecs-1; k >= 0; --k)
282  {
283  Index cornerSize = rows() - k - m_shift;
284  if(m_trans)
285  dst.bottomRightCorner(cornerSize, cornerSize)
286  .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
287  else
288  dst.bottomRightCorner(cornerSize, cornerSize)
289  .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
290  }
291  }
292  }
293 
295  template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
296  {
298  applyThisOnTheRight(dst, workspace);
299  }
300 
302  template<typename Dest, typename Workspace>
303  inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const
304  {
305  workspace.resize(dst.rows());
306  for(Index k = 0; k < m_length; ++k)
307  {
308  Index actual_k = m_trans ? m_length-k-1 : k;
309  dst.rightCols(rows()-m_shift-actual_k)
310  .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
311  }
312  }
313 
315  template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
316  {
318  applyThisOnTheLeft(dst, workspace);
319  }
320 
322  template<typename Dest, typename Workspace>
323  inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const
324  {
325  workspace.resize(dst.cols());
326  for(Index k = 0; k < m_length; ++k)
327  {
328  Index actual_k = m_trans ? k : m_length-k-1;
329  dst.bottomRows(rows()-m_shift-actual_k)
330  .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
331  }
332  }
333 
341  template<typename OtherDerived>
343  {
345  res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
346  applyThisOnTheLeft(res);
347  return res;
348  }
349 
350  template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
351 
362  {
363  m_length = length;
364  return *this;
365  }
366 
379  {
380  m_shift = shift;
381  return *this;
382  }
383 
384  Index length() const { return m_length; }
385  Index shift() const { return m_shift; }
387  /* Necessary for .adjoint() and .conjugate() */
388  template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence;
389 
390  protected:
391 
401  {
402  m_trans = trans;
403  return *this;
404  }
405 
406  bool trans() const { return m_trans; }
408  typename VectorsType::Nested m_vectors;
409  typename CoeffsType::Nested m_coeffs;
410  bool m_trans;
411  Index m_length;
412  Index m_shift;
413 };
414 
423 template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
425 {
427  res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
428  h.applyThisOnTheRight(res);
429  return res;
430 }
431 
436 template<typename VectorsType, typename CoeffsType>
437 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
438 {
440 }
441 
448 template<typename VectorsType, typename CoeffsType>
450 {
452 }
453 
454 } // end namespace Eigen
455 
456 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H