// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2009 Gael Guennebaud // Copyright (C) 2010 Benoit Jacob // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H #define EIGEN_HOUSEHOLDER_SEQUENCE_H namespace Eigen { /** \ingroup Householder_Module * \householder_module * \class HouseholderSequence * \brief Sequence of Householder reflections acting on subspaces with decreasing size * \tparam VectorsType type of matrix containing the Householder vectors * \tparam CoeffsType type of vector containing the Householder coefficients * \tparam Side either OnTheLeft (the default) or OnTheRight * * This class represents a product sequence of Householder reflections where the first Householder reflection * acts on the whole space, the second Householder reflection leaves the one-dimensional subspace spanned by * the first unit vector invariant, the third Householder reflection leaves the two-dimensional subspace * spanned by the first two unit vectors invariant, and so on up to the last reflection which leaves all but * one dimensions invariant and acts only on the last dimension. Such sequences of Householder reflections * are used in several algorithms to zero out certain parts of a matrix. Indeed, the methods * HessenbergDecomposition::matrixQ(), Tridiagonalization::matrixQ(), HouseholderQR::householderQ(), * and ColPivHouseholderQR::householderQ() all return a %HouseholderSequence. * * More precisely, the class %HouseholderSequence represents an \f$n \times n \f$ matrix \f$H \f$ of the * 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 * v_i^* \f$. The i-th Householder coefficient \f$h_i \f$ is a scalar and the i-th Householder vector \f$* v_i \f$ is a vector of the form * \f[ * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ]. * \f] * The last \f$n-i \f$ entries of \f$v_i \f$ are called the essential part of the Householder vector. * * Typical usages are listed below, where H is a HouseholderSequence: * \code * A.applyOnTheRight(H); // A = A * H * A.applyOnTheLeft(H); // A = H * A * A.applyOnTheRight(H.adjoint()); // A = A * H^* * A.applyOnTheLeft(H.adjoint()); // A = H^* * A * MatrixXd Q = H; // conversion to a dense matrix * \endcode * In addition to the adjoint, you can also apply the inverse (=adjoint), the transpose, and the conjugate operators. * * See the documentation for HouseholderSequence(const VectorsType&, const CoeffsType&) for an example. * * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() */ namespace internal { template struct traits > { typedef typename VectorsType::Scalar Scalar; typedef typename VectorsType::Index Index; typedef typename VectorsType::StorageKind StorageKind; enum { RowsAtCompileTime = Side==OnTheLeft ? traits::RowsAtCompileTime : traits::ColsAtCompileTime, ColsAtCompileTime = RowsAtCompileTime, MaxRowsAtCompileTime = Side==OnTheLeft ? traits::MaxRowsAtCompileTime : traits::MaxColsAtCompileTime, MaxColsAtCompileTime = MaxRowsAtCompileTime, Flags = 0 }; }; template struct hseq_side_dependent_impl { typedef Block EssentialVectorType; typedef HouseholderSequence HouseholderSequenceType; typedef typename VectorsType::Index Index; static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k) { Index start = k+1+h.m_shift; return Block(h.m_vectors, start, k, h.rows()-start, 1); } }; template struct hseq_side_dependent_impl { typedef Transpose > EssentialVectorType; typedef HouseholderSequence HouseholderSequenceType; typedef typename VectorsType::Index Index; static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k) { Index start = k+1+h.m_shift; return Block(h.m_vectors, k, start, 1, h.rows()-start).transpose(); } }; template struct matrix_type_times_scalar_type { typedef typename scalar_product_traits::ReturnType ResultScalar; typedef Matrix Type; }; } // end namespace internal template class HouseholderSequence : public EigenBase > { typedef typename internal::hseq_side_dependent_impl::EssentialVectorType EssentialVectorType; public: enum { RowsAtCompileTime = internal::traits::RowsAtCompileTime, ColsAtCompileTime = internal::traits::ColsAtCompileTime, MaxRowsAtCompileTime = internal::traits::MaxRowsAtCompileTime, MaxColsAtCompileTime = internal::traits::MaxColsAtCompileTime }; typedef typename internal::traits::Scalar Scalar; typedef typename VectorsType::Index Index; typedef HouseholderSequence< typename internal::conditional::IsComplex, typename internal::remove_all::type, VectorsType>::type, typename internal::conditional::IsComplex, typename internal::remove_all::type, CoeffsType>::type, Side > ConjugateReturnType; /** \brief Constructor. * \param[in] v %Matrix containing the essential parts of the Householder vectors * \param[in] h Vector containing the Householder coefficients * * Constructs the Householder sequence with coefficients given by \p h and vectors given by \p v. The * i-th Householder coefficient \f$h_i \f$ is given by \p h(i) and the essential part of the i-th * Householder vector \f$v_i \f$ is given by \p v(k,i) with \p k > \p i (the subdiagonal part of the * i-th column). If \p v has fewer columns than rows, then the Householder sequence contains as many * Householder reflections as there are columns. * * \note The %HouseholderSequence object stores \p v and \p h by reference. * * Example: \include HouseholderSequence_HouseholderSequence.cpp * Output: \verbinclude HouseholderSequence_HouseholderSequence.out * * \sa setLength(), setShift() */ HouseholderSequence(const VectorsType& v, const CoeffsType& h) : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()), m_shift(0) { } /** \brief Copy constructor. */ HouseholderSequence(const HouseholderSequence& other) : m_vectors(other.m_vectors), m_coeffs(other.m_coeffs), m_trans(other.m_trans), m_length(other.m_length), m_shift(other.m_shift) { } /** \brief Number of rows of transformation viewed as a matrix. * \returns Number of rows * \details This equals the dimension of the space that the transformation acts on. */ Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); } /** \brief Number of columns of transformation viewed as a matrix. * \returns Number of columns * \details This equals the dimension of the space that the transformation acts on. */ Index cols() const { return rows(); } /** \brief Essential part of a Householder vector. * \param[in] k Index of Householder reflection * \returns Vector containing non-trivial entries of k-th Householder vector * * This function returns the essential part of the Householder vector \f$v_i \f$. This is a vector of * length \f$n-i \f$ containing the last \f$n-i \f$ entries of the vector * \f[ * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ]. * \f] * The index \f$i \f$ equals \p k + shift(), corresponding to the k-th column of the matrix \p v * passed to the constructor. * * \sa setShift(), shift() */ const EssentialVectorType essentialVector(Index k) const { eigen_assert(k >= 0 && k < m_length); return internal::hseq_side_dependent_impl::essentialVector(*this, k); } /** \brief %Transpose of the Householder sequence. */ HouseholderSequence transpose() const { return HouseholderSequence(*this).setTrans(!m_trans); } /** \brief Complex conjugate of the Householder sequence. */ ConjugateReturnType conjugate() const { return ConjugateReturnType(m_vectors.conjugate(), m_coeffs.conjugate()) .setTrans(m_trans) .setLength(m_length) .setShift(m_shift); } /** \brief Adjoint (conjugate transpose) of the Householder sequence. */ ConjugateReturnType adjoint() const { return conjugate().setTrans(!m_trans); } /** \brief Inverse of the Householder sequence (equals the adjoint). */ ConjugateReturnType inverse() const { return adjoint(); } /** \internal */ template inline void evalTo(DestType& dst) const { Matrix workspace(rows()); evalTo(dst, workspace); } /** \internal */ template void evalTo(Dest& dst, Workspace& workspace) const { workspace.resize(rows()); Index vecs = m_length; if( internal::is_same::type,Dest>::value && internal::extract_data(dst) == internal::extract_data(m_vectors)) { // in-place dst.diagonal().setOnes(); dst.template triangularView().setZero(); for(Index k = vecs-1; k >= 0; --k) { Index cornerSize = rows() - k - m_shift; if(m_trans) dst.bottomRightCorner(cornerSize, cornerSize) .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data()); else dst.bottomRightCorner(cornerSize, cornerSize) .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data()); // clear the off diagonal vector dst.col(k).tail(rows()-k-1).setZero(); } // clear the remaining columns if needed for(Index k = 0; k= 0; --k) { Index cornerSize = rows() - k - m_shift; if(m_trans) dst.bottomRightCorner(cornerSize, cornerSize) .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0)); else dst.bottomRightCorner(cornerSize, cornerSize) .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0)); } } } /** \internal */ template inline void applyThisOnTheRight(Dest& dst) const { Matrix workspace(dst.rows()); applyThisOnTheRight(dst, workspace); } /** \internal */ template inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const { workspace.resize(dst.rows()); for(Index k = 0; k < m_length; ++k) { Index actual_k = m_trans ? m_length-k-1 : k; dst.rightCols(rows()-m_shift-actual_k) .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data()); } } /** \internal */ template inline void applyThisOnTheLeft(Dest& dst) const { Matrix workspace(dst.cols()); applyThisOnTheLeft(dst, workspace); } /** \internal */ template inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const { workspace.resize(dst.cols()); for(Index k = 0; k < m_length; ++k) { Index actual_k = m_trans ? k : m_length-k-1; dst.bottomRows(rows()-m_shift-actual_k) .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data()); } } /** \brief Computes the product of a Householder sequence with a matrix. * \param[in] other %Matrix being multiplied. * \returns Expression object representing the product. * * This function computes \f$HM \f$ where \f$H \f$ is the Householder sequence represented by \p *this * and \f$M \f$ is the matrix \p other. */ template typename internal::matrix_type_times_scalar_type::Type operator*(const MatrixBase& other) const { typename internal::matrix_type_times_scalar_type::Type res(other.template cast::ResultScalar>()); applyThisOnTheLeft(res); return res; } template friend struct internal::hseq_side_dependent_impl; /** \brief Sets the length of the Householder sequence. * \param [in] length New value for the length. * * By default, the length \f$n \f$ of the Householder sequence \f$H = H_0 H_1 \ldots H_{n-1} \f$ is set * to the number of columns of the matrix \p v passed to the constructor, or the number of rows if that * is smaller. After this function is called, the length equals \p length. * * \sa length() */ HouseholderSequence& setLength(Index length) { m_length = length; return *this; } /** \brief Sets the shift of the Householder sequence. * \param [in] shift New value for the shift. * * By default, a %HouseholderSequence object represents \f$H = H_0 H_1 \ldots H_{n-1} \f$ and the i-th * column of the matrix \p v passed to the constructor corresponds to the i-th Householder * reflection. After this function is called, the object represents \f$H = H_{\mathrm{shift}} * H_{\mathrm{shift}+1} \ldots H_{n-1} \f$ and the i-th column of \p v corresponds to the (shift+i)-th * Householder reflection. * * \sa shift() */ HouseholderSequence& setShift(Index shift) { m_shift = shift; return *this; } Index length() const { return m_length; } /**< \brief Returns the length of the Householder sequence. */ Index shift() const { return m_shift; } /**< \brief Returns the shift of the Householder sequence. */ /* Necessary for .adjoint() and .conjugate() */ template friend class HouseholderSequence; protected: /** \brief Sets the transpose flag. * \param [in] trans New value of the transpose flag. * * By default, the transpose flag is not set. If the transpose flag is set, then this object represents * \f$H^T = H_{n-1}^T \ldots H_1^T H_0^T \f$ instead of \f$H = H_0 H_1 \ldots H_{n-1} \f$. * * \sa trans() */ HouseholderSequence& setTrans(bool trans) { m_trans = trans; return *this; } bool trans() const { return m_trans; } /**< \brief Returns the transpose flag. */ typename VectorsType::Nested m_vectors; typename CoeffsType::Nested m_coeffs; bool m_trans; Index m_length; Index m_shift; }; /** \brief Computes the product of a matrix with a Householder sequence. * \param[in] other %Matrix being multiplied. * \param[in] h %HouseholderSequence being multiplied. * \returns Expression object representing the product. * * This function computes \f$MH \f$ where \f$M \f$ is the matrix \p other and \f$H \f$ is the * Householder sequence represented by \p h. */ template typename internal::matrix_type_times_scalar_type::Type operator*(const MatrixBase& other, const HouseholderSequence& h) { typename internal::matrix_type_times_scalar_type::Type res(other.template cast::ResultScalar>()); h.applyThisOnTheRight(res); return res; } /** \ingroup Householder_Module \householder_module * \brief Convenience function for constructing a Householder sequence. * \returns A HouseholderSequence constructed from the specified arguments. */ template HouseholderSequence householderSequence(const VectorsType& v, const CoeffsType& h) { return HouseholderSequence(v, h); } /** \ingroup Householder_Module \householder_module * \brief Convenience function for constructing a Householder sequence. * \returns A HouseholderSequence constructed from the specified arguments. * \details This function differs from householderSequence() in that the template argument \p OnTheSide of * the constructed HouseholderSequence is set to OnTheRight, instead of the default OnTheLeft. */ template HouseholderSequence rightHouseholderSequence(const VectorsType& v, const CoeffsType& h) { return HouseholderSequence(v, h); } } // end namespace Eigen #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H