11 #ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H
12 #define EIGEN_HOUSEHOLDER_SEQUENCE_H
59 template<
typename VectorsType,
typename CoeffsType,
int S
ide>
60 struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
62 typedef typename VectorsType::Scalar Scalar;
63 typedef typename VectorsType::StorageIndex StorageIndex;
64 typedef typename VectorsType::StorageKind StorageKind;
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,
76 struct HouseholderSequenceShape {};
78 template<
typename VectorsType,
typename CoeffsType,
int S
ide>
79 struct evaluator_traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
80 :
public evaluator_traits_base<HouseholderSequence<VectorsType,CoeffsType,Side> >
82 typedef HouseholderSequenceShape Shape;
85 template<
typename VectorsType,
typename CoeffsType,
int S
ide>
86 struct hseq_side_dependent_impl
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)
92 Index start = k+1+h.m_shift;
93 return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
97 template<
typename VectorsType,
typename CoeffsType>
98 struct hseq_side_dependent_impl<VectorsType, CoeffsType,
OnTheRight>
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)
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();
109 template<
typename OtherScalarType,
typename MatrixType>
struct matrix_type_times_scalar_type
111 typedef typename scalar_product_traits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
113 typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
114 0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
119 template<
typename VectorsType,
typename CoeffsType,
int S
ide>
class HouseholderSequence
120 :
public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
122 typedef typename internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType EssentialVectorType;
126 RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
127 ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
128 MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
129 MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
131 typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
134 typename internal::conditional<NumTraits<Scalar>::IsComplex,
135 typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type,
137 typename internal::conditional<NumTraits<Scalar>::IsComplex,
138 typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
141 > ConjugateReturnType;
161 : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()),
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)
204 eigen_assert(k >= 0 && k < m_length);
205 return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*
this, k);
217 return ConjugateReturnType(m_vectors.conjugate(), m_coeffs.conjugate())
233 template<
typename DestType>
inline void evalTo(DestType& dst)
const
235 Matrix<Scalar, DestType::RowsAtCompileTime, 1,
237 evalTo(dst, workspace);
241 template<
typename Dest,
typename Workspace>
242 void evalTo(Dest& dst, Workspace& workspace)
const
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))
250 dst.diagonal().setOnes();
251 dst.template triangularView<StrictlyUpper>().setZero();
252 for(
Index k = vecs-1; k >= 0; --k)
256 dst.bottomRightCorner(cornerSize, cornerSize)
257 .applyHouseholderOnTheRight(
essentialVector(k), m_coeffs.coeff(k), workspace.data());
259 dst.bottomRightCorner(cornerSize, cornerSize)
260 .applyHouseholderOnTheLeft(
essentialVector(k), m_coeffs.coeff(k), workspace.data());
263 dst.col(k).tail(
rows()-k-1).setZero();
267 dst.col(k).tail(
rows()-k-1).setZero();
272 for(
Index k = vecs-1; k >= 0; --k)
276 dst.bottomRightCorner(cornerSize, cornerSize)
277 .applyHouseholderOnTheRight(
essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
279 dst.bottomRightCorner(cornerSize, cornerSize)
280 .applyHouseholderOnTheLeft(
essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
286 template<
typename Dest>
inline void applyThisOnTheRight(Dest& dst)
const
288 Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
289 applyThisOnTheRight(dst, workspace);
293 template<
typename Dest,
typename Workspace>
294 inline void applyThisOnTheRight(Dest& dst, Workspace& workspace)
const
296 workspace.resize(dst.rows());
297 for(
Index k = 0; k < m_length; ++k)
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());
306 template<
typename Dest>
inline void applyThisOnTheLeft(Dest& dst)
const
308 Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace(dst.cols());
309 applyThisOnTheLeft(dst, workspace);
313 template<
typename Dest,
typename Workspace>
314 inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace)
const
316 const Index BlockSize = 48;
318 if(m_length>=BlockSize && dst.cols()>1)
320 for(
Index i = 0; i < m_length; i+=BlockSize)
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);
325 Index start = k + m_shift;
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,
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);
339 workspace.resize(dst.cols());
340 for(
Index k = 0; k < m_length; ++k)
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());
356 template<
typename OtherDerived>
360 res(other.template cast<
typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
361 applyThisOnTheLeft(res);
365 template<
typename _VectorsType,
typename _CoeffsType,
int _S
ide>
friend struct internal::hseq_side_dependent_impl;
403 template <
typename VectorsType2,
typename CoeffsType2,
int S
ide2>
friend class HouseholderSequence;
421 bool trans()
const {
return m_trans; }
423 typename VectorsType::Nested m_vectors;
424 typename CoeffsType::Nested m_coeffs;
438 template<
typename OtherDerived,
typename VectorsType,
typename CoeffsType,
int S
ide>
439 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type operator*(
const MatrixBase<OtherDerived>& other,
const HouseholderSequence<VectorsType,CoeffsType,Side>& h)
441 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type
442 res(other.template cast<
typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
443 h.applyThisOnTheRight(res);
451 template<
typename VectorsType,
typename CoeffsType>
463 template<
typename VectorsType,
typename CoeffsType>
471 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
Index length() const
Returns the length of the Householder sequence.
Definition: HouseholderSequence.h:399
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: Constants.h:327
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:258
HouseholderSequence & setTrans(bool trans)
Sets the transpose flag.
Definition: HouseholderSequence.h:415
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
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:316
bool trans() const
Returns the transpose flag.
Definition: HouseholderSequence.h:421
Definition: Constants.h:312
Index shift() const
Returns the shift of the Householder sequence.
Definition: HouseholderSequence.h:400
Definition: Constants.h:325
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