11#ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H
12#define EIGEN_HOUSEHOLDER_SEQUENCE_H
14#include "./InternalHeaderCheck.h"
61template<
typename VectorsType,
typename CoeffsType,
int S
ide>
62struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
64 typedef typename VectorsType::Scalar Scalar;
65 typedef typename VectorsType::StorageIndex StorageIndex;
66 typedef typename VectorsType::StorageKind StorageKind;
68 RowsAtCompileTime = Side==
OnTheLeft ? traits<VectorsType>::RowsAtCompileTime
69 : traits<VectorsType>::ColsAtCompileTime,
70 ColsAtCompileTime = RowsAtCompileTime,
71 MaxRowsAtCompileTime = Side==
OnTheLeft ? traits<VectorsType>::MaxRowsAtCompileTime
72 : traits<VectorsType>::MaxColsAtCompileTime,
73 MaxColsAtCompileTime = MaxRowsAtCompileTime,
78struct HouseholderSequenceShape {};
80template<
typename VectorsType,
typename CoeffsType,
int S
ide>
81struct evaluator_traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
82 :
public evaluator_traits_base<HouseholderSequence<VectorsType,CoeffsType,Side> >
84 typedef HouseholderSequenceShape Shape;
87template<
typename VectorsType,
typename CoeffsType,
int S
ide>
88struct hseq_side_dependent_impl
90 typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType;
91 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType;
92 static EIGEN_DEVICE_FUNC
inline const EssentialVectorType essentialVector(
const HouseholderSequenceType& h,
Index k)
94 Index start = k+1+h.m_shift;
95 return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
99template<
typename VectorsType,
typename CoeffsType>
100struct hseq_side_dependent_impl<VectorsType, CoeffsType,
OnTheRight>
102 typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType;
103 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType;
104 static inline const EssentialVectorType essentialVector(
const HouseholderSequenceType& h,
Index k)
106 Index start = k+1+h.m_shift;
107 return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose();
111template<
typename OtherScalarType,
typename MatrixType>
struct matrix_type_times_scalar_type
113 typedef typename ScalarBinaryOpTraits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
115 typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
116 0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
122 :
public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
128 RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
129 ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
130 MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
131 MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
133 typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
136 typename internal::conditional<NumTraits<Scalar>::IsComplex,
137 typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type,
139 typename internal::conditional<NumTraits<Scalar>::IsComplex,
140 typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
147 typename internal::conditional<NumTraits<Scalar>::IsComplex,
148 typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
154 typename internal::conditional<NumTraits<Scalar>::IsComplex,
155 typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type,
186 : m_vectors(v), m_coeffs(h), m_reverse(false), m_length(v.diagonalSize()),
194 : m_vectors(other.m_vectors),
195 m_coeffs(other.m_coeffs),
196 m_reverse(other.m_reverse),
197 m_length(other.m_length),
198 m_shift(other.m_shift)
206 EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
213 EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
233 eigen_assert(k >= 0 && k < m_length);
234 return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*
this, k);
241 .setReverseFlag(!m_reverse)
250 .setReverseFlag(m_reverse)
260 inline typename internal::conditional<Cond,ConjugateReturnType,ConstHouseholderSequence>::type
263 typedef typename internal::conditional<Cond,ConjugateReturnType,ConstHouseholderSequence>::type ReturnType;
264 return ReturnType(m_vectors.template conjugateIf<Cond>(), m_coeffs.template conjugateIf<Cond>());
271 .setReverseFlag(!m_reverse)
280 template<
typename DestType>
281 inline EIGEN_DEVICE_FUNC
282 void evalTo(DestType& dst)
const
284 Matrix<Scalar, DestType::RowsAtCompileTime, 1,
286 evalTo(dst, workspace);
290 template<
typename Dest,
typename Workspace>
292 void evalTo(Dest& dst, Workspace& workspace)
const
294 workspace.resize(
rows());
295 Index vecs = m_length;
296 if(internal::is_same_dense(dst,m_vectors))
299 dst.diagonal().setOnes();
300 dst.template triangularView<StrictlyUpper>().setZero();
301 for(
Index k = vecs-1; k >= 0; --k)
305 dst.bottomRightCorner(cornerSize, cornerSize)
306 .applyHouseholderOnTheRight(
essentialVector(k), m_coeffs.coeff(k), workspace.data());
308 dst.bottomRightCorner(cornerSize, cornerSize)
309 .applyHouseholderOnTheLeft(
essentialVector(k), m_coeffs.coeff(k), workspace.data());
312 dst.col(k).tail(
rows()-k-1).setZero();
316 dst.col(k).tail(
rows()-k-1).setZero();
318 else if(m_length>BlockSize)
322 applyThisOnTheLeft(dst,workspace,
true);
324 applyThisOnTheLeft(dst,workspace,
true);
329 for(
Index k = vecs-1; k >= 0; --k)
333 dst.bottomRightCorner(cornerSize, cornerSize)
334 .applyHouseholderOnTheRight(
essentialVector(k), m_coeffs.coeff(k), workspace.data());
336 dst.bottomRightCorner(cornerSize, cornerSize)
337 .applyHouseholderOnTheLeft(
essentialVector(k), m_coeffs.coeff(k), workspace.data());
343 template<
typename Dest>
inline void applyThisOnTheRight(Dest& dst)
const
345 Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
346 applyThisOnTheRight(dst, workspace);
350 template<
typename Dest,
typename Workspace>
351 inline void applyThisOnTheRight(Dest& dst, Workspace& workspace)
const
353 workspace.resize(dst.rows());
354 for(
Index k = 0; k < m_length; ++k)
356 Index actual_k = m_reverse ? m_length-k-1 : k;
357 dst.rightCols(
rows()-m_shift-actual_k)
358 .applyHouseholderOnTheRight(
essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
363 template<
typename Dest>
inline void applyThisOnTheLeft(Dest& dst,
bool inputIsIdentity =
false)
const
365 Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace;
366 applyThisOnTheLeft(dst, workspace, inputIsIdentity);
370 template<
typename Dest,
typename Workspace>
371 inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace,
bool inputIsIdentity =
false)
const
373 if(inputIsIdentity && m_reverse)
374 inputIsIdentity =
false;
376 if(m_length>=BlockSize && dst.cols()>1)
379 Index blockSize = m_length<
Index(2*BlockSize) ? (m_length+1)/2 :
Index(BlockSize);
380 for(
Index i = 0; i < m_length; i+=blockSize)
382 Index end = m_reverse ? (std::min)(m_length,i+blockSize) : m_length-i;
383 Index k = m_reverse ? i : (std::max)(
Index(0),
end-blockSize);
385 Index start = k + m_shift;
387 typedef Block<typename internal::remove_all<VectorsType>::type,
Dynamic,
Dynamic> SubVectorsType;
388 SubVectorsType sub_vecs1(m_vectors.const_cast_derived(), Side==
OnTheRight ? k : start,
390 Side==
OnTheRight ? bs : m_vectors.rows()-start,
391 Side==
OnTheRight ? m_vectors.cols()-start : bs);
392 typename internal::conditional<Side==OnTheRight, Transpose<SubVectorsType>, SubVectorsType&>::type sub_vecs(sub_vecs1);
394 Index dstStart = dst.rows()-
rows()+m_shift+k;
396 Block<Dest,Dynamic,Dynamic> sub_dst(dst,
398 inputIsIdentity ? dstStart : 0,
400 inputIsIdentity ? dstRows : dst.cols());
401 apply_block_householder_on_the_left(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_reverse);
406 workspace.resize(dst.cols());
407 for(
Index k = 0; k < m_length; ++k)
409 Index actual_k = m_reverse ? k : m_length-k-1;
410 Index dstStart =
rows()-m_shift-actual_k;
411 dst.bottomRightCorner(dstStart, inputIsIdentity ? dstStart : dst.cols())
412 .applyHouseholderOnTheLeft(
essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
424 template<
typename OtherDerived>
428 res(other.template cast<
typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
429 applyThisOnTheLeft(res, internal::is_identity<OtherDerived>::value && res.rows()==res.cols());
433 template<
typename VectorsType_,
typename CoeffsType_,
int S
ide_>
friend struct internal::hseq_side_dependent_impl;
476 template <
typename VectorsType2,
typename CoeffsType2,
int S
ide2>
friend class HouseholderSequence;
496 bool reverseFlag()
const {
return m_reverse; }
498 typename VectorsType::Nested m_vectors;
499 typename CoeffsType::Nested m_coeffs;
503 enum { BlockSize = 48 };
514template<
typename OtherDerived,
typename VectorsType,
typename CoeffsType,
int S
ide>
518 res(other.template cast<
typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
519 h.applyThisOnTheRight(res);
527template<
typename VectorsType,
typename CoeffsType>
539template<
typename VectorsType,
typename CoeffsType>
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:107
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:67
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:123
AdjointReturnType inverse() const
Inverse of the Householder sequence (equals the adjoint).
Definition: HouseholderSequence.h:277
HouseholderSequence & setLength(Index length)
Sets the length of the Householder sequence.
Definition: HouseholderSequence.h:445
Index shift() const
Returns the shift of the Householder sequence.
Definition: HouseholderSequence.h:473
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Number of columns of transformation viewed as a matrix.
Definition: HouseholderSequence.h:214
HouseholderSequence & setShift(Index shift)
Sets the shift of the Householder sequence.
Definition: HouseholderSequence.h:463
internal::conditional< Cond, ConjugateReturnType, ConstHouseholderSequence >::type conjugateIf() const
Definition: HouseholderSequence.h:261
HouseholderSequence(const HouseholderSequence &other)
Copy constructor.
Definition: HouseholderSequence.h:193
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:425
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Number of rows of transformation viewed as a matrix.
Definition: HouseholderSequence.h:207
Index length() const
Returns the length of the Householder sequence.
Definition: HouseholderSequence.h:470
ConjugateReturnType conjugate() const
Complex conjugate of the Householder sequence.
Definition: HouseholderSequence.h:247
const EssentialVectorType essentialVector(Index k) const
Essential part of a Householder vector.
Definition: HouseholderSequence.h:231
TransposeReturnType transpose() const
Transpose of the Householder sequence.
Definition: HouseholderSequence.h:238
AdjointReturnType adjoint() const
Adjoint (conjugate transpose) of the Householder sequence.
Definition: HouseholderSequence.h:268
HouseholderSequence(const VectorsType &v, const CoeffsType &h)
Constructor.
Definition: HouseholderSequence.h:185
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:182
static const lastp1_t end
Definition: IndexedViewHelper.h:183
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:528
HouseholderSequence< VectorsType, CoeffsType, OnTheRight > rightHouseholderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:540
@ ColMajor
Definition: Constants.h:321
@ AutoAlign
Definition: Constants.h:325
@ OnTheLeft
Definition: Constants.h:334
@ OnTheRight
Definition: Constants.h:336
Namespace containing all symbols from the Eigen library.
Definition: B01_Experimental.dox:1
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Product< MatrixDerived, PermutationDerived, AliasFreeProduct > operator*(const MatrixBase< MatrixDerived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:517
const int Dynamic
Definition: Constants.h:24
Definition: EigenBase.h:32
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41