00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H
00027 #define EIGEN_HOUSEHOLDER_SEQUENCE_H
00028
00029 namespace Eigen {
00030
00072 namespace internal {
00073
00074 template<typename VectorsType, typename CoeffsType, int Side>
00075 struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
00076 {
00077 typedef typename VectorsType::Scalar Scalar;
00078 typedef typename VectorsType::Index Index;
00079 typedef typename VectorsType::StorageKind StorageKind;
00080 enum {
00081 RowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::RowsAtCompileTime
00082 : traits<VectorsType>::ColsAtCompileTime,
00083 ColsAtCompileTime = RowsAtCompileTime,
00084 MaxRowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::MaxRowsAtCompileTime
00085 : traits<VectorsType>::MaxColsAtCompileTime,
00086 MaxColsAtCompileTime = MaxRowsAtCompileTime,
00087 Flags = 0
00088 };
00089 };
00090
00091 template<typename VectorsType, typename CoeffsType, int Side>
00092 struct hseq_side_dependent_impl
00093 {
00094 typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType;
00095 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType;
00096 typedef typename VectorsType::Index Index;
00097 static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
00098 {
00099 Index start = k+1+h.m_shift;
00100 return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
00101 }
00102 };
00103
00104 template<typename VectorsType, typename CoeffsType>
00105 struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight>
00106 {
00107 typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType;
00108 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType;
00109 typedef typename VectorsType::Index Index;
00110 static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
00111 {
00112 Index start = k+1+h.m_shift;
00113 return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose();
00114 }
00115 };
00116
00117 template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type
00118 {
00119 typedef typename scalar_product_traits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
00120 ResultScalar;
00121 typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
00122 0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
00123 };
00124
00125 }
00126
00127 template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence
00128 : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
00129 {
00130 enum {
00131 RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
00132 ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
00133 MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
00134 MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
00135 };
00136 typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
00137 typedef typename VectorsType::Index Index;
00138
00139 typedef typename internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType
00140 EssentialVectorType;
00141
00142 public:
00143
00144 typedef HouseholderSequence<
00145 VectorsType,
00146 typename internal::conditional<NumTraits<Scalar>::IsComplex,
00147 typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
00148 CoeffsType>::type,
00149 Side
00150 > ConjugateReturnType;
00151
00169 HouseholderSequence(const VectorsType& v, const CoeffsType& h)
00170 : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()),
00171 m_shift(0)
00172 {
00173 }
00174
00176 HouseholderSequence(const HouseholderSequence& other)
00177 : m_vectors(other.m_vectors),
00178 m_coeffs(other.m_coeffs),
00179 m_trans(other.m_trans),
00180 m_length(other.m_length),
00181 m_shift(other.m_shift)
00182 {
00183 }
00184
00189 Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
00190
00195 Index cols() const { return rows(); }
00196
00211 const EssentialVectorType essentialVector(Index k) const
00212 {
00213 eigen_assert(k >= 0 && k < m_length);
00214 return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
00215 }
00216
00218 HouseholderSequence transpose() const
00219 {
00220 return HouseholderSequence(*this).setTrans(!m_trans);
00221 }
00222
00224 ConjugateReturnType conjugate() const
00225 {
00226 return ConjugateReturnType(m_vectors, m_coeffs.conjugate())
00227 .setTrans(m_trans)
00228 .setLength(m_length)
00229 .setShift(m_shift);
00230 }
00231
00233 ConjugateReturnType adjoint() const
00234 {
00235 return conjugate().setTrans(!m_trans);
00236 }
00237
00239 ConjugateReturnType inverse() const { return adjoint(); }
00240
00242 template<typename DestType> inline void evalTo(DestType& dst) const
00243 {
00244 Matrix<Scalar, DestType::RowsAtCompileTime, 1,
00245 AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows());
00246 evalTo(dst, workspace);
00247 }
00248
00250 template<typename Dest, typename Workspace>
00251 void evalTo(Dest& dst, Workspace& workspace) const
00252 {
00253 workspace.resize(rows());
00254 Index vecs = m_length;
00255 if( internal::is_same<typename internal::remove_all<VectorsType>::type,Dest>::value
00256 && internal::extract_data(dst) == internal::extract_data(m_vectors))
00257 {
00258
00259 dst.diagonal().setOnes();
00260 dst.template triangularView<StrictlyUpper>().setZero();
00261 for(Index k = vecs-1; k >= 0; --k)
00262 {
00263 Index cornerSize = rows() - k - m_shift;
00264 if(m_trans)
00265 dst.bottomRightCorner(cornerSize, cornerSize)
00266 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
00267 else
00268 dst.bottomRightCorner(cornerSize, cornerSize)
00269 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
00270
00271
00272 dst.col(k).tail(rows()-k-1).setZero();
00273 }
00274
00275 for(Index k = 0; k<cols()-vecs ; ++k)
00276 dst.col(k).tail(rows()-k-1).setZero();
00277 }
00278 else
00279 {
00280 dst.setIdentity(rows(), rows());
00281 for(Index k = vecs-1; k >= 0; --k)
00282 {
00283 Index cornerSize = rows() - k - m_shift;
00284 if(m_trans)
00285 dst.bottomRightCorner(cornerSize, cornerSize)
00286 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
00287 else
00288 dst.bottomRightCorner(cornerSize, cornerSize)
00289 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
00290 }
00291 }
00292 }
00293
00295 template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
00296 {
00297 Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
00298 applyThisOnTheRight(dst, workspace);
00299 }
00300
00302 template<typename Dest, typename Workspace>
00303 inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const
00304 {
00305 workspace.resize(dst.rows());
00306 for(Index k = 0; k < m_length; ++k)
00307 {
00308 Index actual_k = m_trans ? m_length-k-1 : k;
00309 dst.rightCols(rows()-m_shift-actual_k)
00310 .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
00311 }
00312 }
00313
00315 template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
00316 {
00317 Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace(dst.cols());
00318 applyThisOnTheLeft(dst, workspace);
00319 }
00320
00322 template<typename Dest, typename Workspace>
00323 inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const
00324 {
00325 workspace.resize(dst.cols());
00326 for(Index k = 0; k < m_length; ++k)
00327 {
00328 Index actual_k = m_trans ? k : m_length-k-1;
00329 dst.bottomRows(rows()-m_shift-actual_k)
00330 .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
00331 }
00332 }
00333
00341 template<typename OtherDerived>
00342 typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const
00343 {
00344 typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type
00345 res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
00346 applyThisOnTheLeft(res);
00347 return res;
00348 }
00349
00350 template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
00351
00361 HouseholderSequence& setLength(Index length)
00362 {
00363 m_length = length;
00364 return *this;
00365 }
00366
00378 HouseholderSequence& setShift(Index shift)
00379 {
00380 m_shift = shift;
00381 return *this;
00382 }
00383
00384 Index length() const { return m_length; }
00385 Index shift() const { return m_shift; }
00387
00388 template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence;
00389
00390 protected:
00391
00400 HouseholderSequence& setTrans(bool trans)
00401 {
00402 m_trans = trans;
00403 return *this;
00404 }
00405
00406 bool trans() const { return m_trans; }
00408 typename VectorsType::Nested m_vectors;
00409 typename CoeffsType::Nested m_coeffs;
00410 bool m_trans;
00411 Index m_length;
00412 Index m_shift;
00413 };
00414
00423 template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
00424 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence<VectorsType,CoeffsType,Side>& h)
00425 {
00426 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type
00427 res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
00428 h.applyThisOnTheRight(res);
00429 return res;
00430 }
00431
00436 template<typename VectorsType, typename CoeffsType>
00437 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
00438 {
00439 return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h);
00440 }
00441
00448 template<typename VectorsType, typename CoeffsType>
00449 HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h)
00450 {
00451 return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h);
00452 }
00453
00454 }
00455
00456 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H