1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // Copyright (C) 2010 Benoit Jacob <jacob.benoit.1 (at) gmail.com> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H 12 #define EIGEN_HOUSEHOLDER_SEQUENCE_H 13 14 namespace Eigen { 15 16 /** \ingroup Householder_Module 17 * \householder_module 18 * \class HouseholderSequence 19 * \brief Sequence of Householder reflections acting on subspaces with decreasing size 20 * \tparam VectorsType type of matrix containing the Householder vectors 21 * \tparam CoeffsType type of vector containing the Householder coefficients 22 * \tparam Side either OnTheLeft (the default) or OnTheRight 23 * 24 * This class represents a product sequence of Householder reflections where the first Householder reflection 25 * acts on the whole space, the second Householder reflection leaves the one-dimensional subspace spanned by 26 * the first unit vector invariant, the third Householder reflection leaves the two-dimensional subspace 27 * spanned by the first two unit vectors invariant, and so on up to the last reflection which leaves all but 28 * one dimensions invariant and acts only on the last dimension. Such sequences of Householder reflections 29 * are used in several algorithms to zero out certain parts of a matrix. Indeed, the methods 30 * HessenbergDecomposition::matrixQ(), Tridiagonalization::matrixQ(), HouseholderQR::householderQ(), 31 * and ColPivHouseholderQR::householderQ() all return a %HouseholderSequence. 32 * 33 * More precisely, the class %HouseholderSequence represents an \f$ n \times n \f$ matrix \f$ H \f$ of the 34 * 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 35 * v_i^* \f$. The i-th Householder coefficient \f$ h_i \f$ is a scalar and the i-th Householder vector \f$ 36 * v_i \f$ is a vector of the form 37 * \f[ 38 * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ]. 39 * \f] 40 * The last \f$ n-i \f$ entries of \f$ v_i \f$ are called the essential part of the Householder vector. 41 * 42 * Typical usages are listed below, where H is a HouseholderSequence: 43 * \code 44 * A.applyOnTheRight(H); // A = A * H 45 * A.applyOnTheLeft(H); // A = H * A 46 * A.applyOnTheRight(H.adjoint()); // A = A * H^* 47 * A.applyOnTheLeft(H.adjoint()); // A = H^* * A 48 * MatrixXd Q = H; // conversion to a dense matrix 49 * \endcode 50 * In addition to the adjoint, you can also apply the inverse (=adjoint), the transpose, and the conjugate operators. 51 * 52 * See the documentation for HouseholderSequence(const VectorsType&, const CoeffsType&) for an example. 53 * 54 * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() 55 */ 56 57 namespace internal { 58 59 template<typename VectorsType, typename CoeffsType, int Side> 60 struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> > 61 { 62 typedef typename VectorsType::Scalar Scalar; 63 typedef typename VectorsType::Index Index; 64 typedef typename VectorsType::StorageKind StorageKind; 65 enum { 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, 72 Flags = 0 73 }; 74 }; 75 76 template<typename VectorsType, typename CoeffsType, int Side> 77 struct hseq_side_dependent_impl 78 { 79 typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType; 80 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType; 81 typedef typename VectorsType::Index Index; 82 static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k) 83 { 84 Index start = k+1+h.m_shift; 85 return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1); 86 } 87 }; 88 89 template<typename VectorsType, typename CoeffsType> 90 struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight> 91 { 92 typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType; 93 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType; 94 typedef typename VectorsType::Index Index; 95 static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k) 96 { 97 Index start = k+1+h.m_shift; 98 return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose(); 99 } 100 }; 101 102 template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type 103 { 104 typedef typename scalar_product_traits<OtherScalarType, typename MatrixType::Scalar>::ReturnType 105 ResultScalar; 106 typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime, 107 0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type; 108 }; 109 110 } // end namespace internal 111 112 template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence 113 : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> > 114 { 115 typedef typename internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType EssentialVectorType; 116 117 public: 118 enum { 119 RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime, 120 ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime, 121 MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime, 122 MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime 123 }; 124 typedef typename internal::traits<HouseholderSequence>::Scalar Scalar; 125 typedef typename VectorsType::Index Index; 126 127 typedef HouseholderSequence< 128 typename internal::conditional<NumTraits<Scalar>::IsComplex, 129 typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type, 130 VectorsType>::type, 131 typename internal::conditional<NumTraits<Scalar>::IsComplex, 132 typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type, 133 CoeffsType>::type, 134 Side 135 > ConjugateReturnType; 136 137 /** \brief Constructor. 138 * \param[in] v %Matrix containing the essential parts of the Householder vectors 139 * \param[in] h Vector containing the Householder coefficients 140 * 141 * Constructs the Householder sequence with coefficients given by \p h and vectors given by \p v. The 142 * i-th Householder coefficient \f$ h_i \f$ is given by \p h(i) and the essential part of the i-th 143 * Householder vector \f$ v_i \f$ is given by \p v(k,i) with \p k > \p i (the subdiagonal part of the 144 * i-th column). If \p v has fewer columns than rows, then the Householder sequence contains as many 145 * Householder reflections as there are columns. 146 * 147 * \note The %HouseholderSequence object stores \p v and \p h by reference. 148 * 149 * Example: \include HouseholderSequence_HouseholderSequence.cpp 150 * Output: \verbinclude HouseholderSequence_HouseholderSequence.out 151 * 152 * \sa setLength(), setShift() 153 */ 154 HouseholderSequence(const VectorsType& v, const CoeffsType& h) 155 : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()), 156 m_shift(0) 157 { 158 } 159 160 /** \brief Copy constructor. */ 161 HouseholderSequence(const HouseholderSequence& other) 162 : m_vectors(other.m_vectors), 163 m_coeffs(other.m_coeffs), 164 m_trans(other.m_trans), 165 m_length(other.m_length), 166 m_shift(other.m_shift) 167 { 168 } 169 170 /** \brief Number of rows of transformation viewed as a matrix. 171 * \returns Number of rows 172 * \details This equals the dimension of the space that the transformation acts on. 173 */ 174 Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); } 175 176 /** \brief Number of columns of transformation viewed as a matrix. 177 * \returns Number of columns 178 * \details This equals the dimension of the space that the transformation acts on. 179 */ 180 Index cols() const { return rows(); } 181 182 /** \brief Essential part of a Householder vector. 183 * \param[in] k Index of Householder reflection 184 * \returns Vector containing non-trivial entries of k-th Householder vector 185 * 186 * This function returns the essential part of the Householder vector \f$ v_i \f$. This is a vector of 187 * length \f$ n-i \f$ containing the last \f$ n-i \f$ entries of the vector 188 * \f[ 189 * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ]. 190 * \f] 191 * The index \f$ i \f$ equals \p k + shift(), corresponding to the k-th column of the matrix \p v 192 * passed to the constructor. 193 * 194 * \sa setShift(), shift() 195 */ 196 const EssentialVectorType essentialVector(Index k) const 197 { 198 eigen_assert(k >= 0 && k < m_length); 199 return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k); 200 } 201 202 /** \brief %Transpose of the Householder sequence. */ 203 HouseholderSequence transpose() const 204 { 205 return HouseholderSequence(*this).setTrans(!m_trans); 206 } 207 208 /** \brief Complex conjugate of the Householder sequence. */ 209 ConjugateReturnType conjugate() const 210 { 211 return ConjugateReturnType(m_vectors.conjugate(), m_coeffs.conjugate()) 212 .setTrans(m_trans) 213 .setLength(m_length) 214 .setShift(m_shift); 215 } 216 217 /** \brief Adjoint (conjugate transpose) of the Householder sequence. */ 218 ConjugateReturnType adjoint() const 219 { 220 return conjugate().setTrans(!m_trans); 221 } 222 223 /** \brief Inverse of the Householder sequence (equals the adjoint). */ 224 ConjugateReturnType inverse() const { return adjoint(); } 225 226 /** \internal */ 227 template<typename DestType> inline void evalTo(DestType& dst) const 228 { 229 Matrix<Scalar, DestType::RowsAtCompileTime, 1, 230 AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows()); 231 evalTo(dst, workspace); 232 } 233 234 /** \internal */ 235 template<typename Dest, typename Workspace> 236 void evalTo(Dest& dst, Workspace& workspace) const 237 { 238 workspace.resize(rows()); 239 Index vecs = m_length; 240 if( internal::is_same<typename internal::remove_all<VectorsType>::type,Dest>::value 241 && internal::extract_data(dst) == internal::extract_data(m_vectors)) 242 { 243 // in-place 244 dst.diagonal().setOnes(); 245 dst.template triangularView<StrictlyUpper>().setZero(); 246 for(Index k = vecs-1; k >= 0; --k) 247 { 248 Index cornerSize = rows() - k - m_shift; 249 if(m_trans) 250 dst.bottomRightCorner(cornerSize, cornerSize) 251 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data()); 252 else 253 dst.bottomRightCorner(cornerSize, cornerSize) 254 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data()); 255 256 // clear the off diagonal vector 257 dst.col(k).tail(rows()-k-1).setZero(); 258 } 259 // clear the remaining columns if needed 260 for(Index k = 0; k<cols()-vecs ; ++k) 261 dst.col(k).tail(rows()-k-1).setZero(); 262 } 263 else 264 { 265 dst.setIdentity(rows(), rows()); 266 for(Index k = vecs-1; k >= 0; --k) 267 { 268 Index cornerSize = rows() - k - m_shift; 269 if(m_trans) 270 dst.bottomRightCorner(cornerSize, cornerSize) 271 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0)); 272 else 273 dst.bottomRightCorner(cornerSize, cornerSize) 274 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0)); 275 } 276 } 277 } 278 279 /** \internal */ 280 template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const 281 { 282 Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows()); 283 applyThisOnTheRight(dst, workspace); 284 } 285 286 /** \internal */ 287 template<typename Dest, typename Workspace> 288 inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const 289 { 290 workspace.resize(dst.rows()); 291 for(Index k = 0; k < m_length; ++k) 292 { 293 Index actual_k = m_trans ? m_length-k-1 : k; 294 dst.rightCols(rows()-m_shift-actual_k) 295 .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data()); 296 } 297 } 298 299 /** \internal */ 300 template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const 301 { 302 Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace(dst.cols()); 303 applyThisOnTheLeft(dst, workspace); 304 } 305 306 /** \internal */ 307 template<typename Dest, typename Workspace> 308 inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const 309 { 310 workspace.resize(dst.cols()); 311 for(Index k = 0; k < m_length; ++k) 312 { 313 Index actual_k = m_trans ? k : m_length-k-1; 314 dst.bottomRows(rows()-m_shift-actual_k) 315 .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data()); 316 } 317 } 318 319 /** \brief Computes the product of a Householder sequence with a matrix. 320 * \param[in] other %Matrix being multiplied. 321 * \returns Expression object representing the product. 322 * 323 * This function computes \f$ HM \f$ where \f$ H \f$ is the Householder sequence represented by \p *this 324 * and \f$ M \f$ is the matrix \p other. 325 */ 326 template<typename OtherDerived> 327 typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const 328 { 329 typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type 330 res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>()); 331 applyThisOnTheLeft(res); 332 return res; 333 } 334 335 template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl; 336 337 /** \brief Sets the length of the Householder sequence. 338 * \param [in] length New value for the length. 339 * 340 * By default, the length \f$ n \f$ of the Householder sequence \f$ H = H_0 H_1 \ldots H_{n-1} \f$ is set 341 * to the number of columns of the matrix \p v passed to the constructor, or the number of rows if that 342 * is smaller. After this function is called, the length equals \p length. 343 * 344 * \sa length() 345 */ 346 HouseholderSequence& setLength(Index length) 347 { 348 m_length = length; 349 return *this; 350 } 351 352 /** \brief Sets the shift of the Householder sequence. 353 * \param [in] shift New value for the shift. 354 * 355 * By default, a %HouseholderSequence object represents \f$ H = H_0 H_1 \ldots H_{n-1} \f$ and the i-th 356 * column of the matrix \p v passed to the constructor corresponds to the i-th Householder 357 * reflection. After this function is called, the object represents \f$ H = H_{\mathrm{shift}} 358 * H_{\mathrm{shift}+1} \ldots H_{n-1} \f$ and the i-th column of \p v corresponds to the (shift+i)-th 359 * Householder reflection. 360 * 361 * \sa shift() 362 */ 363 HouseholderSequence& setShift(Index shift) 364 { 365 m_shift = shift; 366 return *this; 367 } 368 369 Index length() const { return m_length; } /**< \brief Returns the length of the Householder sequence. */ 370 Index shift() const { return m_shift; } /**< \brief Returns the shift of the Householder sequence. */ 371 372 /* Necessary for .adjoint() and .conjugate() */ 373 template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence; 374 375 protected: 376 377 /** \brief Sets the transpose flag. 378 * \param [in] trans New value of the transpose flag. 379 * 380 * By default, the transpose flag is not set. If the transpose flag is set, then this object represents 381 * \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$. 382 * 383 * \sa trans() 384 */ 385 HouseholderSequence& setTrans(bool trans) 386 { 387 m_trans = trans; 388 return *this; 389 } 390 391 bool trans() const { return m_trans; } /**< \brief Returns the transpose flag. */ 392 393 typename VectorsType::Nested m_vectors; 394 typename CoeffsType::Nested m_coeffs; 395 bool m_trans; 396 Index m_length; 397 Index m_shift; 398 }; 399 400 /** \brief Computes the product of a matrix with a Householder sequence. 401 * \param[in] other %Matrix being multiplied. 402 * \param[in] h %HouseholderSequence being multiplied. 403 * \returns Expression object representing the product. 404 * 405 * This function computes \f$ MH \f$ where \f$ M \f$ is the matrix \p other and \f$ H \f$ is the 406 * Householder sequence represented by \p h. 407 */ 408 template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side> 409 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence<VectorsType,CoeffsType,Side>& h) 410 { 411 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type 412 res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>()); 413 h.applyThisOnTheRight(res); 414 return res; 415 } 416 417 /** \ingroup Householder_Module \householder_module 418 * \brief Convenience function for constructing a Householder sequence. 419 * \returns A HouseholderSequence constructed from the specified arguments. 420 */ 421 template<typename VectorsType, typename CoeffsType> 422 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h) 423 { 424 return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h); 425 } 426 427 /** \ingroup Householder_Module \householder_module 428 * \brief Convenience function for constructing a Householder sequence. 429 * \returns A HouseholderSequence constructed from the specified arguments. 430 * \details This function differs from householderSequence() in that the template argument \p OnTheSide of 431 * the constructed HouseholderSequence is set to OnTheRight, instead of the default OnTheLeft. 432 */ 433 template<typename VectorsType, typename CoeffsType> 434 HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h) 435 { 436 return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h); 437 } 438 439 } // end namespace Eigen 440 441 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H 442
