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 EIGEN_DEVICE_FUNC
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 ScalarBinaryOpTraits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
113 typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
114 0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
120 :
public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
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,
145 typename internal::conditional<NumTraits<Scalar>::IsComplex,
146 typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
152 typename internal::conditional<NumTraits<Scalar>::IsComplex,
153 typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type,
160 typename internal::add_const<VectorsType>::type,
161 typename internal::add_const<CoeffsType>::type,
184 : m_vectors(v), m_coeffs(h), m_reverse(false), m_length(v.diagonalSize()),
192 : m_vectors(other.m_vectors),
193 m_coeffs(other.m_coeffs),
194 m_reverse(other.m_reverse),
195 m_length(other.m_length),
196 m_shift(other.m_shift)
204 EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
211 EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
231 eigen_assert(k >= 0 && k < m_length);
232 return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*
this, k);
239 .setReverseFlag(!m_reverse)
248 .setReverseFlag(m_reverse)
258 inline typename internal::conditional<Cond,ConjugateReturnType,ConstHouseholderSequence>::type
261 typedef typename internal::conditional<Cond,ConjugateReturnType,ConstHouseholderSequence>::type ReturnType;
262 return ReturnType(m_vectors.template conjugateIf<Cond>(), m_coeffs.template conjugateIf<Cond>());
269 .setReverseFlag(!m_reverse)
278 template<
typename DestType>
279 inline EIGEN_DEVICE_FUNC
280 void evalTo(DestType& dst)
const
282 Matrix<Scalar, DestType::RowsAtCompileTime, 1,
284 evalTo(dst, workspace);
288 template<
typename Dest,
typename Workspace>
290 void evalTo(Dest& dst, Workspace& workspace)
const
292 workspace.resize(
rows());
293 Index vecs = m_length;
294 if(internal::is_same_dense(dst,m_vectors))
297 dst.diagonal().setOnes();
298 dst.template triangularView<StrictlyUpper>().setZero();
299 for(
Index k = vecs-1; k >= 0; --k)
303 dst.bottomRightCorner(cornerSize, cornerSize)
304 .applyHouseholderOnTheRight(
essentialVector(k), m_coeffs.coeff(k), workspace.data());
306 dst.bottomRightCorner(cornerSize, cornerSize)
307 .applyHouseholderOnTheLeft(
essentialVector(k), m_coeffs.coeff(k), workspace.data());
310 dst.col(k).tail(
rows()-k-1).setZero();
314 dst.col(k).tail(
rows()-k-1).setZero();
316 else if(m_length>BlockSize)
320 applyThisOnTheLeft(dst,workspace,
true);
322 applyThisOnTheLeft(dst,workspace,
true);
327 for(
Index k = vecs-1; k >= 0; --k)
331 dst.bottomRightCorner(cornerSize, cornerSize)
332 .applyHouseholderOnTheRight(
essentialVector(k), m_coeffs.coeff(k), workspace.data());
334 dst.bottomRightCorner(cornerSize, cornerSize)
335 .applyHouseholderOnTheLeft(
essentialVector(k), m_coeffs.coeff(k), workspace.data());
341 template<
typename Dest>
inline void applyThisOnTheRight(Dest& dst)
const
343 Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
344 applyThisOnTheRight(dst, workspace);
348 template<
typename Dest,
typename Workspace>
349 inline void applyThisOnTheRight(Dest& dst, Workspace& workspace)
const
351 workspace.resize(dst.rows());
352 for(
Index k = 0; k < m_length; ++k)
354 Index actual_k = m_reverse ? m_length-k-1 : k;
355 dst.rightCols(
rows()-m_shift-actual_k)
356 .applyHouseholderOnTheRight(
essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
361 template<
typename Dest>
inline void applyThisOnTheLeft(Dest& dst,
bool inputIsIdentity =
false)
const
363 Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace;
364 applyThisOnTheLeft(dst, workspace, inputIsIdentity);
368 template<
typename Dest,
typename Workspace>
369 inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace,
bool inputIsIdentity =
false)
const
371 if(inputIsIdentity && m_reverse)
372 inputIsIdentity =
false;
374 if(m_length>=BlockSize && dst.cols()>1)
377 Index blockSize = m_length<
Index(2*BlockSize) ? (m_length+1)/2 :
Index(BlockSize);
378 for(
Index i = 0; i < m_length; i+=blockSize)
380 Index end = m_reverse ? (std::min)(m_length,i+blockSize) : m_length-i;
381 Index k = m_reverse ? i : (std::max)(
Index(0),end-blockSize);
383 Index start = k + m_shift;
385 typedef Block<typename internal::remove_all<VectorsType>::type,
Dynamic,
Dynamic> SubVectorsType;
386 SubVectorsType sub_vecs1(m_vectors.const_cast_derived(), Side==
OnTheRight ? k : start,
388 Side==
OnTheRight ? bs : m_vectors.rows()-start,
389 Side==
OnTheRight ? m_vectors.cols()-start : bs);
390 typename internal::conditional<Side==OnTheRight, Transpose<SubVectorsType>, SubVectorsType&>::type sub_vecs(sub_vecs1);
392 Index dstStart = dst.rows()-
rows()+m_shift+k;
394 Block<Dest,Dynamic,Dynamic> sub_dst(dst,
396 inputIsIdentity ? dstStart : 0,
398 inputIsIdentity ? dstRows : dst.cols());
399 apply_block_householder_on_the_left(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_reverse);
404 workspace.resize(dst.cols());
405 for(
Index k = 0; k < m_length; ++k)
407 Index actual_k = m_reverse ? k : m_length-k-1;
408 Index dstStart =
rows()-m_shift-actual_k;
409 dst.bottomRightCorner(dstStart, inputIsIdentity ? dstStart : dst.cols())
410 .applyHouseholderOnTheLeft(
essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
422 template<
typename OtherDerived>
426 res(other.template cast<
typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
427 applyThisOnTheLeft(res, internal::is_identity<OtherDerived>::value && res.rows()==res.cols());
431 template<
typename _VectorsType,
typename _CoeffsType,
int _S
ide>
friend struct internal::hseq_side_dependent_impl;
468 Index length()
const {
return m_length; }
471 Index shift()
const {
return m_shift; }
474 template <
typename VectorsType2,
typename CoeffsType2,
int S
ide2>
friend class HouseholderSequence;
494 bool reverseFlag()
const {
return m_reverse; }
496 typename VectorsType::Nested m_vectors;
497 typename CoeffsType::Nested m_coeffs;
501 enum { BlockSize = 48 };
512 template<
typename OtherDerived,
typename VectorsType,
typename CoeffsType,
int S
ide>
516 res(other.template cast<
typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
517 h.applyThisOnTheRight(res);
525 template<
typename VectorsType,
typename CoeffsType>
537 template<
typename VectorsType,
typename CoeffsType>
545 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H