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
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070 namespace internal {
00071
00072 template<typename VectorsType, typename CoeffsType, int Side>
00073 struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
00074 {
00075 typedef typename VectorsType::Scalar Scalar;
00076 typedef typename VectorsType::Index Index;
00077 typedef typename VectorsType::StorageKind StorageKind;
00078 enum {
00079 RowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::RowsAtCompileTime
00080 : traits<VectorsType>::ColsAtCompileTime,
00081 ColsAtCompileTime = RowsAtCompileTime,
00082 MaxRowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::MaxRowsAtCompileTime
00083 : traits<VectorsType>::MaxColsAtCompileTime,
00084 MaxColsAtCompileTime = MaxRowsAtCompileTime,
00085 Flags = 0
00086 };
00087 };
00088
00089 template<typename VectorsType, typename CoeffsType, int Side>
00090 struct hseq_side_dependent_impl
00091 {
00092 typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType;
00093 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType;
00094 typedef typename VectorsType::Index Index;
00095 static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
00096 {
00097 Index start = k+1+h.m_shift;
00098 return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
00099 }
00100 };
00101
00102 template<typename VectorsType, typename CoeffsType>
00103 struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight>
00104 {
00105 typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType;
00106 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType;
00107 typedef typename VectorsType::Index Index;
00108 static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
00109 {
00110 Index start = k+1+h.m_shift;
00111 return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose();
00112 }
00113 };
00114
00115 template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type
00116 {
00117 typedef typename scalar_product_traits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
00118 ResultScalar;
00119 typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
00120 0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
00121 };
00122
00123 }
00124
00125 template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence
00126 : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
00127 {
00128 enum {
00129 RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
00130 ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
00131 MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
00132 MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
00133 };
00134 typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
00135 typedef typename VectorsType::Index Index;
00136
00137 typedef typename internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType
00138 EssentialVectorType;
00139
00140 public:
00141
00142 typedef HouseholderSequence<
00143 VectorsType,
00144 typename internal::conditional<NumTraits<Scalar>::IsComplex,
00145 typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
00146 CoeffsType>::type,
00147 Side
00148 > ConjugateReturnType;
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167 HouseholderSequence(const VectorsType& v, const CoeffsType& h)
00168 : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()),
00169 m_shift(0)
00170 {
00171 }
00172
00173
00174 HouseholderSequence(const HouseholderSequence& other)
00175 : m_vectors(other.m_vectors),
00176 m_coeffs(other.m_coeffs),
00177 m_trans(other.m_trans),
00178 m_length(other.m_length),
00179 m_shift(other.m_shift)
00180 {
00181 }
00182
00183
00184
00185
00186
00187 Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
00188
00189
00190
00191
00192
00193 Index cols() const { return rows(); }
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209 const EssentialVectorType essentialVector(Index k) const
00210 {
00211 eigen_assert(k >= 0 && k < m_length);
00212 return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
00213 }
00214
00215
00216 HouseholderSequence transpose() const
00217 {
00218 return HouseholderSequence(*this).setTrans(!m_trans);
00219 }
00220
00221
00222 ConjugateReturnType conjugate() const
00223 {
00224 return ConjugateReturnType(m_vectors, m_coeffs.conjugate())
00225 .setTrans(m_trans)
00226 .setLength(m_length)
00227 .setShift(m_shift);
00228 }
00229
00230
00231 ConjugateReturnType adjoint() const
00232 {
00233 return conjugate().setTrans(!m_trans);
00234 }
00235
00236
00237 ConjugateReturnType inverse() const { return adjoint(); }
00238
00239
00240 template<typename DestType> void evalTo(DestType& dst) const
00241 {
00242 Index vecs = m_length;
00243
00244 Matrix<Scalar, DestType::RowsAtCompileTime, 1,
00245 AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> temp(rows());
00246 if( internal::is_same<typename internal::remove_all<VectorsType>::type,DestType>::value
00247 && internal::extract_data(dst) == internal::extract_data(m_vectors))
00248 {
00249
00250 dst.diagonal().setOnes();
00251 dst.template triangularView<StrictlyUpper>().setZero();
00252 for(Index k = vecs-1; k >= 0; --k)
00253 {
00254 Index cornerSize = rows() - k - m_shift;
00255 if(m_trans)
00256 dst.bottomRightCorner(cornerSize, cornerSize)
00257 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
00258 else
00259 dst.bottomRightCorner(cornerSize, cornerSize)
00260 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
00261
00262
00263 dst.col(k).tail(rows()-k-1).setZero();
00264 }
00265
00266 for(Index k = 0; k<cols()-vecs ; ++k)
00267 dst.col(k).tail(rows()-k-1).setZero();
00268 }
00269 else
00270 {
00271 dst.setIdentity(rows(), rows());
00272 for(Index k = vecs-1; k >= 0; --k)
00273 {
00274 Index cornerSize = rows() - k - m_shift;
00275 if(m_trans)
00276 dst.bottomRightCorner(cornerSize, cornerSize)
00277 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
00278 else
00279 dst.bottomRightCorner(cornerSize, cornerSize)
00280 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
00281 }
00282 }
00283 }
00284
00285
00286 template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
00287 {
00288 Matrix<Scalar,1,Dest::RowsAtCompileTime> temp(dst.rows());
00289 for(Index k = 0; k < m_length; ++k)
00290 {
00291 Index actual_k = m_trans ? m_length-k-1 : k;
00292 dst.rightCols(rows()-m_shift-actual_k)
00293 .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
00294 }
00295 }
00296
00297
00298 template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
00299 {
00300 Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols());
00301 for(Index k = 0; k < m_length; ++k)
00302 {
00303 Index actual_k = m_trans ? k : m_length-k-1;
00304 dst.bottomRows(rows()-m_shift-actual_k)
00305 .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
00306 }
00307 }
00308
00309
00310
00311
00312
00313
00314
00315
00316 template<typename OtherDerived>
00317 typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const
00318 {
00319 typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type
00320 res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::ResultScalar>());
00321 applyThisOnTheLeft(res);
00322 return res;
00323 }
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333 template<typename OtherDerived> friend
00334 typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence& h)
00335 {
00336 typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type
00337 res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::ResultScalar>());
00338 h.applyThisOnTheRight(res);
00339 return res;
00340 }
00341
00342 template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352 HouseholderSequence& setTrans(bool trans)
00353 {
00354 m_trans = trans;
00355 return *this;
00356 }
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367 HouseholderSequence& setLength(Index length)
00368 {
00369 m_length = length;
00370 return *this;
00371 }
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384 HouseholderSequence& setShift(Index shift)
00385 {
00386 m_shift = shift;
00387 return *this;
00388 }
00389
00390 bool trans() const { return m_trans; }
00391 Index length() const { return m_length; }
00392 Index shift() const { return m_shift; }
00393
00394 protected:
00395 typename VectorsType::Nested m_vectors;
00396 typename CoeffsType::Nested m_coeffs;
00397 bool m_trans;
00398 Index m_length;
00399 Index m_shift;
00400 };
00401
00402
00403
00404
00405
00406 template<typename VectorsType, typename CoeffsType>
00407 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
00408 {
00409 return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h);
00410 }
00411
00412
00413
00414
00415
00416
00417
00418 template<typename VectorsType, typename CoeffsType>
00419 HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h)
00420 {
00421 return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h);
00422 }
00423
00424 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H