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 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,
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);
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_dense(dst,m_vectors))
249 dst.diagonal().setOnes();
250 dst.template triangularView<StrictlyUpper>().setZero();
251 for(
Index k = vecs-1; k >= 0; --k)
255 dst.bottomRightCorner(cornerSize, cornerSize)
256 .applyHouseholderOnTheRight(
essentialVector(k), m_coeffs.coeff(k), workspace.data());
258 dst.bottomRightCorner(cornerSize, cornerSize)
259 .applyHouseholderOnTheLeft(
essentialVector(k), m_coeffs.coeff(k), workspace.data());
262 dst.col(k).tail(
rows()-k-1).setZero();
266 dst.col(k).tail(
rows()-k-1).setZero();
271 for(
Index k = vecs-1; k >= 0; --k)
275 dst.bottomRightCorner(cornerSize, cornerSize)
276 .applyHouseholderOnTheRight(
essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
278 dst.bottomRightCorner(cornerSize, cornerSize)
279 .applyHouseholderOnTheLeft(
essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
285 template<
typename Dest>
inline void applyThisOnTheRight(Dest& dst)
const
287 Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
288 applyThisOnTheRight(dst, workspace);
292 template<
typename Dest,
typename Workspace>
293 inline void applyThisOnTheRight(Dest& dst, Workspace& workspace)
const
295 workspace.resize(dst.rows());
296 for(
Index k = 0; k < m_length; ++k)
298 Index actual_k = m_trans ? m_length-k-1 : k;
299 dst.rightCols(
rows()-m_shift-actual_k)
300 .applyHouseholderOnTheRight(
essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
305 template<
typename Dest>
inline void applyThisOnTheLeft(Dest& dst)
const
307 Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace;
308 applyThisOnTheLeft(dst, workspace);
312 template<
typename Dest,
typename Workspace>
313 inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace)
const
315 const Index BlockSize = 48;
317 if(m_length>=BlockSize && dst.cols()>1)
319 for(
Index i = 0; i < m_length; i+=BlockSize)
321 Index end = m_trans ? (std::min)(m_length,i+BlockSize) : m_length-i;
322 Index k = m_trans ? i : (std::max)(
Index(0),end-BlockSize);
324 Index start = k + m_shift;
326 typedef Block<typename internal::remove_all<VectorsType>::type,
Dynamic,
Dynamic> SubVectorsType;
327 SubVectorsType sub_vecs1(m_vectors.const_cast_derived(), Side==
OnTheRight ? k : start,
329 Side==
OnTheRight ? bs : m_vectors.rows()-start,
330 Side==
OnTheRight ? m_vectors.cols()-start : bs);
331 typename internal::conditional<Side==OnTheRight, Transpose<SubVectorsType>, SubVectorsType&>::type sub_vecs(sub_vecs1);
332 Block<Dest,Dynamic,Dynamic> sub_dst(dst,dst.rows()-
rows()+m_shift+k,0,
rows()-m_shift-k,dst.cols());
333 apply_block_householder_on_the_left(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_trans);
338 workspace.resize(dst.cols());
339 for(
Index k = 0; k < m_length; ++k)
341 Index actual_k = m_trans ? k : m_length-k-1;
342 dst.bottomRows(
rows()-m_shift-actual_k)
343 .applyHouseholderOnTheLeft(
essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
355 template<
typename OtherDerived>
359 res(other.template cast<
typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
360 applyThisOnTheLeft(res);
364 template<
typename _VectorsType,
typename _CoeffsType,
int _S
ide>
friend struct internal::hseq_side_dependent_impl;
402 template <
typename VectorsType2,
typename CoeffsType2,
int S
ide2>
friend class HouseholderSequence;
420 bool trans()
const {
return m_trans; }
422 typename VectorsType::Nested m_vectors;
423 typename CoeffsType::Nested m_coeffs;
437 template<
typename OtherDerived,
typename VectorsType,
typename CoeffsType,
int S
ide>
441 res(other.template cast<
typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
442 h.applyThisOnTheRight(res);
450 template<
typename VectorsType,
typename CoeffsType>
462 template<
typename VectorsType,
typename CoeffsType>
470 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H