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
14namespace Eigen {
15
56
57namespace internal {
58
59template<typename VectorsType, typename CoeffsType, int Side>
60struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
61{
62 typedef typename VectorsType::Scalar Scalar;
63 typedef typename VectorsType::Index Index;
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
76template<typename VectorsType, typename CoeffsType, int Side>
77struct hseq_side_dependent_impl
78{
79 typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType;
80 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType;
81 typedef typename VectorsType::Index Index;
82 static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
83 {
84 Index start = k+1+h.m_shift;
85 return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
86 }
87};
88
89template<typename VectorsType, typename CoeffsType>
90struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight>
91{
92 typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType;
93 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType;
94 typedef typename VectorsType::Index Index;
95 static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
96 {
97 Index start = k+1+h.m_shift;
98 return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose();
99 }
100};
101
102template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type
103{
104 typedef typename scalar_product_traits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
105 ResultScalar;
106 typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
107 0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
108};
109
110} // end namespace internal
111
112template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence
113 : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
114{
115 enum {
116 RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
117 ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
118 MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
119 MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
120 };
121 typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
122 typedef typename VectorsType::Index Index;
123
124 typedef typename internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType
125 EssentialVectorType;
126
127 public:
128
129 typedef HouseholderSequence<
130 VectorsType,
131 typename internal::conditional<NumTraits<Scalar>::IsComplex,
132 typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
133 CoeffsType>::type,
134 Side
135 > ConjugateReturnType;
136
154 HouseholderSequence(const VectorsType& v, const CoeffsType& h)
155 : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()),
156 m_shift(0)
157 {
158 }
159
161 HouseholderSequence(const HouseholderSequence& other)
162 : m_vectors(other.m_vectors),
163 m_coeffs(other.m_coeffs),
164 m_trans(other.m_trans),
165 m_length(other.m_length),
166 m_shift(other.m_shift)
167 {
168 }
169
174 Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
175
180 Index cols() const { return rows(); }
181
196 const EssentialVectorType essentialVector(Index k) const
197 {
198 eigen_assert(k >= 0 && k < m_length);
199 return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
200 }
201
203 HouseholderSequence transpose() const
204 {
205 return HouseholderSequence(*this).setTrans(!m_trans);
206 }
207
209 ConjugateReturnType conjugate() const
210 {
211 return ConjugateReturnType(m_vectors, m_coeffs.conjugate())
212 .setTrans(m_trans)
213 .setLength(m_length)
214 .setShift(m_shift);
215 }
216
218 ConjugateReturnType adjoint() const
219 {
220 return conjugate().setTrans(!m_trans);
221 }
222
224 ConjugateReturnType inverse() const { return adjoint(); }
225
227 template<typename DestType> inline void evalTo(DestType& dst) const
228 {
229 Matrix<Scalar, DestType::RowsAtCompileTime, 1,
230 AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows());
231 evalTo(dst, workspace);
232 }
233
235 template<typename Dest, typename Workspace>
236 void evalTo(Dest& dst, Workspace& workspace) const
237 {
238 workspace.resize(rows());
239 Index vecs = m_length;
240 if( internal::is_same<typename internal::remove_all<VectorsType>::type,Dest>::value
241 && internal::extract_data(dst) == internal::extract_data(m_vectors))
242 {
243 // in-place
244 dst.diagonal().setOnes();
245 dst.template triangularView<StrictlyUpper>().setZero();
246 for(Index k = vecs-1; k >= 0; --k)
247 {
248 Index cornerSize = rows() - k - m_shift;
249 if(m_trans)
250 dst.bottomRightCorner(cornerSize, cornerSize)
251 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
252 else
253 dst.bottomRightCorner(cornerSize, cornerSize)
254 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
255
256 // clear the off diagonal vector
257 dst.col(k).tail(rows()-k-1).setZero();
258 }
259 // clear the remaining columns if needed
260 for(Index k = 0; k<cols()-vecs ; ++k)
261 dst.col(k).tail(rows()-k-1).setZero();
262 }
263 else
264 {
265 dst.setIdentity(rows(), rows());
266 for(Index k = vecs-1; k >= 0; --k)
267 {
268 Index cornerSize = rows() - k - m_shift;
269 if(m_trans)
270 dst.bottomRightCorner(cornerSize, cornerSize)
271 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
272 else
273 dst.bottomRightCorner(cornerSize, cornerSize)
274 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
275 }
276 }
277 }
278
280 template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
281 {
282 Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
283 applyThisOnTheRight(dst, workspace);
284 }
285
287 template<typename Dest, typename Workspace>
288 inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const
289 {
290 workspace.resize(dst.rows());
291 for(Index k = 0; k < m_length; ++k)
292 {
293 Index actual_k = m_trans ? m_length-k-1 : k;
294 dst.rightCols(rows()-m_shift-actual_k)
295 .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
296 }
297 }
298
300 template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
301 {
302 Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace(dst.cols());
303 applyThisOnTheLeft(dst, workspace);
304 }
305
307 template<typename Dest, typename Workspace>
308 inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const
309 {
310 workspace.resize(dst.cols());
311 for(Index k = 0; k < m_length; ++k)
312 {
313 Index actual_k = m_trans ? k : m_length-k-1;
314 dst.bottomRows(rows()-m_shift-actual_k)
315 .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
316 }
317 }
318
326 template<typename OtherDerived>
327 typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const
328 {
329 typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type
330 res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
331 applyThisOnTheLeft(res);
332 return res;
333 }
334
335 template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
336
346 HouseholderSequence& setLength(Index length)
347 {
348 m_length = length;
349 return *this;
350 }
351
363 HouseholderSequence& setShift(Index shift)
364 {
365 m_shift = shift;
366 return *this;
367 }
368
369 Index length() const { return m_length; }
370 Index shift() const { return m_shift; }
371
372 /* Necessary for .adjoint() and .conjugate() */
373 template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence;
374
375 protected:
376
385 HouseholderSequence& setTrans(bool trans)
386 {
387 m_trans = trans;
388 return *this;
389 }
390
391 bool trans() const { return m_trans; }
392
393 typename VectorsType::Nested m_vectors;
394 typename CoeffsType::Nested m_coeffs;
395 bool m_trans;
396 Index m_length;
397 Index m_shift;
398};
399
408template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
409typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence<VectorsType,CoeffsType,Side>& h)
410{
411 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type
412 res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
413 h.applyThisOnTheRight(res);
414 return res;
415}
416
421template<typename VectorsType, typename CoeffsType>
422HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
423{
424 return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h);
425}
426
433template<typename VectorsType, typename CoeffsType>
434HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h)
435{
436 return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h);
437}
438
439} // end namespace Eigen
440
441#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
Eigen::HouseholderSequence
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition HouseholderSequence.h:114
Eigen::HouseholderSequence::cols
Index cols() const
Number of columns of transformation viewed as a matrix.
Definition HouseholderSequence.h:180
Eigen::HouseholderSequence::essentialVector
const EssentialVectorType essentialVector(Index k) const
Essential part of a Householder vector.
Definition HouseholderSequence.h:196
Eigen::HouseholderSequence::HouseholderSequence
HouseholderSequence(const HouseholderSequence &other)
Copy constructor.
Definition HouseholderSequence.h:161
Eigen::HouseholderSequence::transpose
HouseholderSequence transpose() const
Transpose of the Householder sequence.
Definition HouseholderSequence.h:203
Eigen::HouseholderSequence::setShift
HouseholderSequence & setShift(Index shift)
Sets the shift of the Householder sequence.
Definition HouseholderSequence.h:363
Eigen::HouseholderSequence::HouseholderSequence
HouseholderSequence(const VectorsType &v, const CoeffsType &h)
Constructor.
Definition HouseholderSequence.h:154
Eigen::HouseholderSequence::setLength
HouseholderSequence & setLength(Index length)
Sets the length of the Householder sequence.
Definition HouseholderSequence.h:346
Eigen::ConjugateReturnType::length
Index length() const
Definition HouseholderSequence.h:369
Eigen::HouseholderSequence::inverse
ConjugateReturnType inverse() const
Inverse of the Householder sequence (equals the adjoint).
Definition HouseholderSequence.h:224
Eigen::HouseholderSequence::adjoint
ConjugateReturnType adjoint() const
Adjoint (conjugate transpose) of the Householder sequence.
Definition HouseholderSequence.h:218
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:327
Eigen::HouseholderSequence::conjugate
ConjugateReturnType conjugate() const
Complex conjugate of the Householder sequence.
Definition HouseholderSequence.h:209
Eigen::HouseholderSequence::rows
Index rows() const
Number of rows of transformation viewed as a matrix.
Definition HouseholderSequence.h:174
Eigen::HouseholderSequence::setTrans
HouseholderSequence & setTrans(bool trans)
Sets the transpose flag.
Definition HouseholderSequence.h:385
Eigen::ConjugateReturnType::trans
bool trans() const
Definition HouseholderSequence.h:391
Eigen::ConjugateReturnType::shift
Index shift() const
Definition HouseholderSequence.h:370
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition MatrixBase.h:50
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition Matrix.h:129
Eigen::ColMajor
@ ColMajor
Definition Constants.h:257
Eigen::OnTheLeft
@ OnTheLeft
Definition Constants.h:270
Eigen::OnTheRight
@ OnTheRight
Definition Constants.h:272
Eigen::internal
Definition LDLT.h:18
Eigen::EigenBase
Definition EigenBase.h:27