src/Householder/HouseholderSequence.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 // Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H
0012 #define EIGEN_HOUSEHOLDER_SEQUENCE_H
0013
0014 namespace Eigen {
0015
0016 /** \ingroup Householder_Module
0017   * \householder_module
0018   * \class HouseholderSequence
0019   * \brief Sequence of Householder reflections acting on subspaces with decreasing size
0020   * \tparam VectorsType type of matrix containing the Householder vectors
0021   * \tparam CoeffsType  type of vector containing the Householder coefficients
0022   * \tparam Side        either OnTheLeft (the default) or OnTheRight
0023   *
0024   * This class represents a product sequence of Householder reflections where the first Householder reflection
0025   * acts on the whole space, the second Householder reflection leaves the one-dimensional subspace spanned by
0026   * the first unit vector invariant, the third Householder reflection leaves the two-dimensional subspace
0027   * spanned by the first two unit vectors invariant, and so on up to the last reflection which leaves all but
0028   * one dimensions invariant and acts only on the last dimension. Such sequences of Householder reflections
0029   * are used in several algorithms to zero out certain parts of a matrix. Indeed, the methods
0030   * HessenbergDecomposition::matrixQ(), Tridiagonalization::matrixQ(), HouseholderQR::householderQ(),
0031   * and ColPivHouseholderQR::householderQ() all return a %HouseholderSequence.
0032   *
0033   * More precisely, the class %HouseholderSequence represents an \f$ n \times n \f$ matrix \f$ H \f$ of the
0034   * form \f$ H = \prod_{i=0}^{n-1} H_i \f$ where the i-th Householder reflection is \f$ H_i = I - h_i v_i
0035   * v_i^* \f$. The i-th Householder coefficient \f$ h_i \f$ is a scalar and the i-th Householder vector \f$
0036   * v_i \f$ is a vector of the form
0037   * \f[
0038   * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ].
0039   * \f]
0040   * The last \f$ n-i \f$ entries of \f$ v_i \f$ are called the essential part of the Householder vector.
0041   *
0042   * Typical usages are listed below, where H is a HouseholderSequence:
0043   * \code
0044   * A.applyOnTheRight(H);             // A = A * H
0045   * A.applyOnTheLeft(H);              // A = H * A
0046   * A.applyOnTheRight(H.adjoint());   // A = A * H^*
0047   * A.applyOnTheLeft(H.adjoint());    // A = H^* * A
0048   * MatrixXd Q = H;                   // conversion to a dense matrix
0049   * \endcode
0050   * In addition to the adjoint, you can also apply the inverse (=adjoint), the transpose, and the conjugate operators.
0051   *
0052   * See the documentation for HouseholderSequence(const VectorsType&, const CoeffsType&) for an example.
0053   *
0054   * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
0055   */
0056
0057 namespace internal {
0058
0059 template<typename VectorsType, typename CoeffsType, int Side>
0060 struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
0061 {
0062   typedef typename VectorsType::Scalar Scalar;
0063   typedef typename VectorsType::StorageIndex StorageIndex;
0064   typedef typename VectorsType::StorageKind StorageKind;
0065   enum {
0066     RowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::RowsAtCompileTime
0067                                         : traits<VectorsType>::ColsAtCompileTime,
0068     ColsAtCompileTime = RowsAtCompileTime,
0069     MaxRowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::MaxRowsAtCompileTime
0070                                            : traits<VectorsType>::MaxColsAtCompileTime,
0071     MaxColsAtCompileTime = MaxRowsAtCompileTime,
0072     Flags = 0
0073   };
0074 };
0075
0076 struct HouseholderSequenceShape {};
0077
0078 template<typename VectorsType, typename CoeffsType, int Side>
0079 struct evaluator_traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
0080   : public evaluator_traits_base<HouseholderSequence<VectorsType,CoeffsType,Side> >
0081 {
0082   typedef HouseholderSequenceShape Shape;
0083 };
0084
0085 template<typename VectorsType, typename CoeffsType, int Side>
0086 struct hseq_side_dependent_impl
0087 {
0088   typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType;
0089   typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType;
0090   static EIGEN_DEVICE_FUNC inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
0091   {
0092     Index start = k+1+h.m_shift;
0093     return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
0094   }
0095 };
0096
0097 template<typename VectorsType, typename CoeffsType>
0098 struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight>
0099 {
0100   typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType;
0101   typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType;
0102   static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
0103   {
0104     Index start = k+1+h.m_shift;
0105     return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose();
0106   }
0107 };
0108
0109 template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type
0110 {
0111   typedef typename ScalarBinaryOpTraits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
0112     ResultScalar;
0113   typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
0114                  0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
0115 };
0116
0117 } // end namespace internal
0118
0119 template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence
0120   : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
0121 {
0122     typedef typename internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType EssentialVectorType;
0123
0124   public:
0125     enum {
0126       RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
0127       ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
0128       MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
0129       MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
0130     };
0131     typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
0132
0133     typedef HouseholderSequence<
0134       typename internal::conditional<NumTraits<Scalar>::IsComplex,
0135         typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type,
0136         VectorsType>::type,
0137       typename internal::conditional<NumTraits<Scalar>::IsComplex,
0138         typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
0139         CoeffsType>::type,
0140       Side
0141     > ConjugateReturnType;
0142
0143     typedef HouseholderSequence<
0144       VectorsType,
0145       typename internal::conditional<NumTraits<Scalar>::IsComplex,
0146         typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
0147         CoeffsType>::type,
0148       Side
0149     > AdjointReturnType;
0150
0151     typedef HouseholderSequence<
0152       typename internal::conditional<NumTraits<Scalar>::IsComplex,
0153         typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type,
0154         VectorsType>::type,
0155       CoeffsType,
0156       Side
0157     > TransposeReturnType;
0158
0159     typedef HouseholderSequence<
0160       typename internal::add_const<VectorsType>::type,
0161       typename internal::add_const<CoeffsType>::type,
0162       Side
0163     > ConstHouseholderSequence;
0164
0165     /** \brief Constructor.
0166       * \param[in]  v      %Matrix containing the essential parts of the Householder vectors
0167       * \param[in]  h      Vector containing the Householder coefficients
0168       *
0169       * Constructs the Householder sequence with coefficients given by \p h and vectors given by \p v. The
0170       * i-th Householder coefficient \f$ h_i \f$ is given by \p h(i) and the essential part of the i-th
0171       * Householder vector \f$ v_i \f$ is given by \p v(k,i) with \p k > \p i (the subdiagonal part of the
0172       * i-th column). If \p v has fewer columns than rows, then the Householder sequence contains as many
0173       * Householder reflections as there are columns.
0174       *
0175       * \note The %HouseholderSequence object stores \p v and \p h by reference.
0176       *
0177       * Example: \include HouseholderSequence_HouseholderSequence.cpp
0178       * Output: \verbinclude HouseholderSequence_HouseholderSequence.out
0179       *
0180       * \sa setLength(), setShift()
0181       */
0182     EIGEN_DEVICE_FUNC
0183     HouseholderSequence(const VectorsType& v, const CoeffsType& h)
0184       : m_vectors(v), m_coeffs(h), m_reverse(false), m_length(v.diagonalSize()),
0185         m_shift(0)
0186     {
0187     }
0188
0189     /** \brief Copy constructor. */
0190     EIGEN_DEVICE_FUNC
0191     HouseholderSequence(const HouseholderSequence& other)
0192       : m_vectors(other.m_vectors),
0193         m_coeffs(other.m_coeffs),
0194         m_reverse(other.m_reverse),
0195         m_length(other.m_length),
0196         m_shift(other.m_shift)
0197     {
0198     }
0199
0200     /** \brief Number of rows of transformation viewed as a matrix.
0201       * \returns Number of rows
0202       * \details This equals the dimension of the space that the transformation acts on.
0203       */
0204     EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
0205     Index rows() const EIGEN_NOEXCEPT { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
0206
0207     /** \brief Number of columns of transformation viewed as a matrix.
0208       * \returns Number of columns
0209       * \details This equals the dimension of the space that the transformation acts on.
0210       */
0211     EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
0212     Index cols() const EIGEN_NOEXCEPT { return rows(); }
0213
0214     /** \brief Essential part of a Householder vector.
0215       * \param[in]  k  Index of Householder reflection
0216       * \returns    Vector containing non-trivial entries of k-th Householder vector
0217       *
0218       * This function returns the essential part of the Householder vector \f$ v_i \f$. This is a vector of
0219       * length \f$ n-i \f$ containing the last \f$ n-i \f$ entries of the vector
0220       * \f[
0221       * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ].
0222       * \f]
0223       * The index \f$ i \f$ equals \p k + shift(), corresponding to the k-th column of the matrix \p v
0224       * passed to the constructor.
0225       *
0226       * \sa setShift(), shift()
0227       */
0228     EIGEN_DEVICE_FUNC
0229     const EssentialVectorType essentialVector(Index k) const
0230     {
0231       eigen_assert(k >= 0 && k < m_length);
0232       return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
0233     }
0234
0235     /** \brief %Transpose of the Householder sequence. */
0236     TransposeReturnType transpose() const
0237     {
0238       return TransposeReturnType(m_vectors.conjugate(), m_coeffs)
0239               .setReverseFlag(!m_reverse)
0240               .setLength(m_length)
0241               .setShift(m_shift);
0242     }
0243
0244     /** \brief Complex conjugate of the Householder sequence. */
0245     ConjugateReturnType conjugate() const
0246     {
0247       return ConjugateReturnType(m_vectors.conjugate(), m_coeffs.conjugate())
0248              .setReverseFlag(m_reverse)
0249              .setLength(m_length)
0250              .setShift(m_shift);
0251     }
0252
0253     /** \returns an expression of the complex conjugate of \c *this if Cond==true,
0254      *           returns \c *this otherwise.
0255      */
0256     template<bool Cond>
0257     EIGEN_DEVICE_FUNC
0258     inline typename internal::conditional<Cond,ConjugateReturnType,ConstHouseholderSequence>::type
0259     conjugateIf() const
0260     {
0261       typedef typename internal::conditional<Cond,ConjugateReturnType,ConstHouseholderSequence>::type ReturnType;
0262       return ReturnType(m_vectors.template conjugateIf<Cond>(), m_coeffs.template conjugateIf<Cond>());
0263     }
0264
0265     /** \brief Adjoint (conjugate transpose) of the Householder sequence. */
0266     AdjointReturnType adjoint() const
0267     {
0268       return AdjointReturnType(m_vectors, m_coeffs.conjugate())
0269               .setReverseFlag(!m_reverse)
0270               .setLength(m_length)
0271               .setShift(m_shift);
0272     }
0273
0274     /** \brief Inverse of the Householder sequence (equals the adjoint). */
0275     AdjointReturnType inverse() const { return adjoint(); }
0276
0277     /** \internal */
0278     template<typename DestType>
0279     inline EIGEN_DEVICE_FUNC
0280     void evalTo(DestType& dst) const
0281     {
0282       Matrix<Scalar, DestType::RowsAtCompileTime, 1,
0283              AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows());
0284       evalTo(dst, workspace);
0285     }
0286
0287     /** \internal */
0288     template<typename Dest, typename Workspace>
0289     EIGEN_DEVICE_FUNC
0290     void evalTo(Dest& dst, Workspace& workspace) const
0291     {
0292       workspace.resize(rows());
0293       Index vecs = m_length;
0294       if(internal::is_same_dense(dst,m_vectors))
0295       {
0296         // in-place
0297         dst.diagonal().setOnes();
0298         dst.template triangularView<StrictlyUpper>().setZero();
0299         for(Index k = vecs-1; k >= 0; --k)
0300         {
0301           Index cornerSize = rows() - k - m_shift;
0302           if(m_reverse)
0303             dst.bottomRightCorner(cornerSize, cornerSize)
0304                .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
0305           else
0306             dst.bottomRightCorner(cornerSize, cornerSize)
0307                .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
0308
0309           // clear the off diagonal vector
0310           dst.col(k).tail(rows()-k-1).setZero();
0311         }
0312         // clear the remaining columns if needed
0313         for(Index k = 0; k<cols()-vecs ; ++k)
0314           dst.col(k).tail(rows()-k-1).setZero();
0315       }
0316       else if(m_length>BlockSize)
0317       {
0318         dst.setIdentity(rows(), rows());
0319         if(m_reverse)
0320           applyThisOnTheLeft(dst,workspace,true);
0321         else
0322           applyThisOnTheLeft(dst,workspace,true);
0323       }
0324       else
0325       {
0326         dst.setIdentity(rows(), rows());
0327         for(Index k = vecs-1; k >= 0; --k)
0328         {
0329           Index cornerSize = rows() - k - m_shift;
0330           if(m_reverse)
0331             dst.bottomRightCorner(cornerSize, cornerSize)
0332                .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
0333           else
0334             dst.bottomRightCorner(cornerSize, cornerSize)
0335                .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
0336         }
0337       }
0338     }
0339
0340     /** \internal */
0341     template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
0342     {
0343       Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
0344       applyThisOnTheRight(dst, workspace);
0345     }
0346
0347     /** \internal */
0348     template<typename Dest, typename Workspace>
0349     inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const
0350     {
0351       workspace.resize(dst.rows());
0352       for(Index k = 0; k < m_length; ++k)
0353       {
0354         Index actual_k = m_reverse ? m_length-k-1 : k;
0355         dst.rightCols(rows()-m_shift-actual_k)
0356            .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
0357       }
0358     }
0359
0360     /** \internal */
0361     template<typename Dest> inline void applyThisOnTheLeft(Dest& dst, bool inputIsIdentity = false) const
0362     {
0363       Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace;
0364       applyThisOnTheLeft(dst, workspace, inputIsIdentity);
0365     }
0366
0367     /** \internal */
0368     template<typename Dest, typename Workspace>
0369     inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace, bool inputIsIdentity = false) const
0370     {
0371       if(inputIsIdentity && m_reverse)
0372         inputIsIdentity = false;
0373       // if the entries are large enough, then apply the reflectors by block
0374       if(m_length>=BlockSize && dst.cols()>1)
0375       {
0376         // Make sure we have at least 2 useful blocks, otherwise it is point-less:
0377         Index blockSize = m_length<Index(2*BlockSize) ? (m_length+1)/2 : Index(BlockSize);
0378         for(Index i = 0; i < m_length; i+=blockSize)
0379         {
0380           Index end = m_reverse ? (std::min)(m_length,i+blockSize) : m_length-i;
0381           Index k = m_reverse ? i : (std::max)(Index(0),end-blockSize);
0382           Index bs = end-k;
0383           Index start = k + m_shift;
0384
0385           typedef Block<typename internal::remove_all<VectorsType>::type,Dynamic,Dynamic> SubVectorsType;
0386           SubVectorsType sub_vecs1(m_vectors.const_cast_derived(), Side==OnTheRight ? k : start,
0387                                                                    Side==OnTheRight ? start : k,
0388                                                                    Side==OnTheRight ? bs : m_vectors.rows()-start,
0389                                                                    Side==OnTheRight ? m_vectors.cols()-start : bs);
0390           typename internal::conditional<Side==OnTheRight, Transpose<SubVectorsType>, SubVectorsType&>::type sub_vecs(sub_vecs1);
0391
0392           Index dstStart = dst.rows()-rows()+m_shift+k;
0393           Index dstRows  = rows()-m_shift-k;
0394           Block<Dest,Dynamic,Dynamic> sub_dst(dst,
0395                                               dstStart,
0396                                               inputIsIdentity ? dstStart : 0,
0397                                               dstRows,
0398                                               inputIsIdentity ? dstRows : dst.cols());
0399           apply_block_householder_on_the_left(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_reverse);
0400         }
0401       }
0402       else
0403       {
0404         workspace.resize(dst.cols());
0405         for(Index k = 0; k < m_length; ++k)
0406         {
0407           Index actual_k = m_reverse ? k : m_length-k-1;
0408           Index dstStart = rows()-m_shift-actual_k;
0409           dst.bottomRightCorner(dstStart, inputIsIdentity ? dstStart : dst.cols())
0410             .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
0411         }
0412       }
0413     }
0414
0415     /** \brief Computes the product of a Householder sequence with a matrix.
0416       * \param[in]  other  %Matrix being multiplied.
0417       * \returns    Expression object representing the product.
0418       *
0419       * This function computes \f$ HM \f$ where \f$ H \f$ is the Householder sequence represented by \p *this
0420       * and \f$ M \f$ is the matrix \p other.
0421       */
0422     template<typename OtherDerived>
0423     typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const
0424     {
0425       typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type
0426         res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
0427       applyThisOnTheLeft(res, internal::is_identity<OtherDerived>::value && res.rows()==res.cols());
0428       return res;
0429     }
0430
0431     template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
0432
0433     /** \brief Sets the length of the Householder sequence.
0434       * \param [in]  length  New value for the length.
0435       *
0436       * By default, the length \f$ n \f$ of the Householder sequence \f$ H = H_0 H_1 \ldots H_{n-1} \f$ is set
0437       * to the number of columns of the matrix \p v passed to the constructor, or the number of rows if that
0438       * is smaller. After this function is called, the length equals \p length.
0439       *
0440       * \sa length()
0441       */
0442     EIGEN_DEVICE_FUNC
0443     HouseholderSequence& setLength(Index length)
0444     {
0445       m_length = length;
0446       return *this;
0447     }
0448
0449     /** \brief Sets the shift of the Householder sequence.
0450       * \param [in]  shift  New value for the shift.
0451       *
0452       * By default, a %HouseholderSequence object represents \f$ H = H_0 H_1 \ldots H_{n-1} \f$ and the i-th
0453       * column of the matrix \p v passed to the constructor corresponds to the i-th Householder
0454       * reflection. After this function is called, the object represents \f$ H = H_{\mathrm{shift}}
0455       * H_{\mathrm{shift}+1} \ldots H_{n-1} \f$ and the i-th column of \p v corresponds to the (shift+i)-th
0456       * Householder reflection.
0457       *
0458       * \sa shift()
0459       */
0460     EIGEN_DEVICE_FUNC
0461     HouseholderSequence& setShift(Index shift)
0462     {
0463       m_shift = shift;
0464       return *this;
0465     }
0466
0467     EIGEN_DEVICE_FUNC
0468     Index length() const { return m_length; }  /**< \brief Returns the length of the Householder sequence. */
0469
0470     EIGEN_DEVICE_FUNC
0471     Index shift() const { return m_shift; }    /**< \brief Returns the shift of the Householder sequence. */
0472
0473     /* Necessary for .adjoint() and .conjugate() */
0474     template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence;
0475
0476   protected:
0477
0478     /** \internal
0479       * \brief Sets the reverse flag.
0480       * \param [in]  reverse  New value of the reverse flag.
0481       *
0482       * By default, the reverse flag is not set. If the reverse flag is set, then this object represents
0483       * \f$ H^r = H_{n-1} \ldots H_1 H_0 \f$ instead of \f$ H = H_0 H_1 \ldots H_{n-1} \f$.
0484       * \note For real valued HouseholderSequence this is equivalent to transposing \f$ H \f$.
0485       *
0486       * \sa reverseFlag(), transpose(), adjoint()
0487       */
0488     HouseholderSequence& setReverseFlag(bool reverse)
0489     {
0490       m_reverse = reverse;
0491       return *this;
0492     }
0493
0494     bool reverseFlag() const { return m_reverse; }     /**< \internal \brief Returns the reverse flag. */
0495
0496     typename VectorsType::Nested m_vectors;
0497     typename CoeffsType::Nested m_coeffs;
0498     bool m_reverse;
0499     Index m_length;
0500     Index m_shift;
0501     enum { BlockSize = 48 };
0502 };
0503
0504 /** \brief Computes the product of a matrix with a Householder sequence.
0505   * \param[in]  other  %Matrix being multiplied.
0506   * \param[in]  h      %HouseholderSequence being multiplied.
0507   * \returns    Expression object representing the product.
0508   *
0509   * This function computes \f$ MH \f$ where \f$ M \f$ is the matrix \p other and \f$ H \f$ is the
0510   * Householder sequence represented by \p h.
0511   */
0512 template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
0513 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence<VectorsType,CoeffsType,Side>& h)
0514 {
0515   typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type
0516     res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
0517   h.applyThisOnTheRight(res);
0518   return res;
0519 }
0520
0521 /** \ingroup Householder_Module \householder_module
0522   * \brief Convenience function for constructing a Householder sequence.
0523   * \returns A HouseholderSequence constructed from the specified arguments.
0524   */
0525 template<typename VectorsType, typename CoeffsType>
0526 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
0527 {
0528   return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h);
0529 }
0530
0531 /** \ingroup Householder_Module \householder_module
0532   * \brief Convenience function for constructing a Householder sequence.
0533   * \returns A HouseholderSequence constructed from the specified arguments.
0534   * \details This function differs from householderSequence() in that the template argument \p OnTheSide of
0535   * the constructed HouseholderSequence is set to OnTheRight, instead of the default OnTheLeft.
0536   */
0537 template<typename VectorsType, typename CoeffsType>
0538 HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h)
0539 {
0540   return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h);
0541 }
0542
0543 } // end namespace Eigen
0544
0545 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H