![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H 00012 #define EIGEN_HOUSEHOLDER_SEQUENCE_H 00013 00014 namespace Eigen { 00015 00057 namespace internal { 00058 00059 template<typename VectorsType, typename CoeffsType, int Side> 00060 struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> > 00061 { 00062 typedef typename VectorsType::Scalar Scalar; 00063 typedef typename VectorsType::StorageIndex StorageIndex; 00064 typedef typename VectorsType::StorageKind StorageKind; 00065 enum { 00066 RowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::RowsAtCompileTime 00067 : traits<VectorsType>::ColsAtCompileTime, 00068 ColsAtCompileTime = RowsAtCompileTime, 00069 MaxRowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::MaxRowsAtCompileTime 00070 : traits<VectorsType>::MaxColsAtCompileTime, 00071 MaxColsAtCompileTime = MaxRowsAtCompileTime, 00072 Flags = 0 00073 }; 00074 }; 00075 00076 struct HouseholderSequenceShape {}; 00077 00078 template<typename VectorsType, typename CoeffsType, int Side> 00079 struct evaluator_traits<HouseholderSequence<VectorsType,CoeffsType,Side> > 00080 : public evaluator_traits_base<HouseholderSequence<VectorsType,CoeffsType,Side> > 00081 { 00082 typedef HouseholderSequenceShape Shape; 00083 }; 00084 00085 template<typename VectorsType, typename CoeffsType, int Side> 00086 struct hseq_side_dependent_impl 00087 { 00088 typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType; 00089 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType; 00090 static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k) 00091 { 00092 Index start = k+1+h.m_shift; 00093 return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1); 00094 } 00095 }; 00096 00097 template<typename VectorsType, typename CoeffsType> 00098 struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight> 00099 { 00100 typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType; 00101 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType; 00102 static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k) 00103 { 00104 Index start = k+1+h.m_shift; 00105 return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose(); 00106 } 00107 }; 00108 00109 template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type 00110 { 00111 typedef typename ScalarBinaryOpTraits<OtherScalarType, typename MatrixType::Scalar>::ReturnType 00112 ResultScalar; 00113 typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime, 00114 0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type; 00115 }; 00116 00117 } // end namespace internal 00118 00119 template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence 00120 : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> > 00121 { 00122 typedef typename internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType EssentialVectorType; 00123 00124 public: 00125 enum { 00126 RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime, 00127 ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime, 00128 MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime, 00129 MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime 00130 }; 00131 typedef typename internal::traits<HouseholderSequence>::Scalar Scalar; 00132 00133 typedef HouseholderSequence< 00134 typename internal::conditional<NumTraits<Scalar>::IsComplex, 00135 typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type, 00136 VectorsType>::type, 00137 typename internal::conditional<NumTraits<Scalar>::IsComplex, 00138 typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type, 00139 CoeffsType>::type, 00140 Side 00141 > ConjugateReturnType; 00142 00160 HouseholderSequence(const VectorsType& v, const CoeffsType& h) 00161 : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()), 00162 m_shift(0) 00163 { 00164 } 00165 00167 HouseholderSequence(const HouseholderSequence& other) 00168 : m_vectors(other.m_vectors), 00169 m_coeffs(other.m_coeffs), 00170 m_trans(other.m_trans), 00171 m_length(other.m_length), 00172 m_shift(other.m_shift) 00173 { 00174 } 00175 00180 Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); } 00181 00186 Index cols() const { return rows(); } 00187 00202 const EssentialVectorType essentialVector(Index k) const 00203 { 00204 eigen_assert(k >= 0 && k < m_length); 00205 return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k); 00206 } 00207 00209 HouseholderSequence transpose() const 00210 { 00211 return HouseholderSequence(*this).setTrans(!m_trans); 00212 } 00213 00215 ConjugateReturnType conjugate() const 00216 { 00217 return ConjugateReturnType(m_vectors.conjugate(), m_coeffs.conjugate()) 00218 .setTrans(m_trans) 00219 .setLength(m_length) 00220 .setShift(m_shift); 00221 } 00222 00224 ConjugateReturnType adjoint() const 00225 { 00226 return conjugate().setTrans(!m_trans); 00227 } 00228 00230 ConjugateReturnType inverse() const { return adjoint(); } 00231 00233 template<typename DestType> inline void evalTo(DestType& dst) const 00234 { 00235 Matrix<Scalar, DestType::RowsAtCompileTime, 1, 00236 AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows()); 00237 evalTo(dst, workspace); 00238 } 00239 00241 template<typename Dest, typename Workspace> 00242 void evalTo(Dest& dst, Workspace& workspace) const 00243 { 00244 workspace.resize(rows()); 00245 Index vecs = m_length; 00246 if(internal::is_same_dense(dst,m_vectors)) 00247 { 00248 // in-place 00249 dst.diagonal().setOnes(); 00250 dst.template triangularView<StrictlyUpper>().setZero(); 00251 for(Index k = vecs-1; k >= 0; --k) 00252 { 00253 Index cornerSize = rows() - k - m_shift; 00254 if(m_trans) 00255 dst.bottomRightCorner(cornerSize, cornerSize) 00256 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data()); 00257 else 00258 dst.bottomRightCorner(cornerSize, cornerSize) 00259 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data()); 00260 00261 // clear the off diagonal vector 00262 dst.col(k).tail(rows()-k-1).setZero(); 00263 } 00264 // clear the remaining columns if needed 00265 for(Index k = 0; k<cols()-vecs ; ++k) 00266 dst.col(k).tail(rows()-k-1).setZero(); 00267 } 00268 else 00269 { 00270 dst.setIdentity(rows(), rows()); 00271 for(Index k = vecs-1; k >= 0; --k) 00272 { 00273 Index cornerSize = rows() - k - m_shift; 00274 if(m_trans) 00275 dst.bottomRightCorner(cornerSize, cornerSize) 00276 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0)); 00277 else 00278 dst.bottomRightCorner(cornerSize, cornerSize) 00279 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0)); 00280 } 00281 } 00282 } 00283 00285 template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const 00286 { 00287 Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows()); 00288 applyThisOnTheRight(dst, workspace); 00289 } 00290 00292 template<typename Dest, typename Workspace> 00293 inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const 00294 { 00295 workspace.resize(dst.rows()); 00296 for(Index k = 0; k < m_length; ++k) 00297 { 00298 Index actual_k = m_trans ? m_length-k-1 : k; 00299 dst.rightCols(rows()-m_shift-actual_k) 00300 .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data()); 00301 } 00302 } 00303 00305 template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const 00306 { 00307 Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace; 00308 applyThisOnTheLeft(dst, workspace); 00309 } 00310 00312 template<typename Dest, typename Workspace> 00313 inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const 00314 { 00315 const Index BlockSize = 48; 00316 // if the entries are large enough, then apply the reflectors by block 00317 if(m_length>=BlockSize && dst.cols()>1) 00318 { 00319 for(Index i = 0; i < m_length; i+=BlockSize) 00320 { 00321 Index end = m_trans ? (std::min)(m_length,i+BlockSize) : m_length-i; 00322 Index k = m_trans ? i : (std::max)(Index(0),end-BlockSize); 00323 Index bs = end-k; 00324 Index start = k + m_shift; 00325 00326 typedef Block<typename internal::remove_all<VectorsType>::type,Dynamic,Dynamic> SubVectorsType; 00327 SubVectorsType sub_vecs1(m_vectors.const_cast_derived(), Side==OnTheRight ? k : start, 00328 Side==OnTheRight ? start : k, 00329 Side==OnTheRight ? bs : m_vectors.rows()-start, 00330 Side==OnTheRight ? m_vectors.cols()-start : bs); 00331 typename internal::conditional<Side==OnTheRight, Transpose<SubVectorsType>, SubVectorsType&>::type sub_vecs(sub_vecs1); 00332 Block<Dest,Dynamic,Dynamic> sub_dst(dst,dst.rows()-rows()+m_shift+k,0, rows()-m_shift-k,dst.cols()); 00333 apply_block_householder_on_the_left(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_trans); 00334 } 00335 } 00336 else 00337 { 00338 workspace.resize(dst.cols()); 00339 for(Index k = 0; k < m_length; ++k) 00340 { 00341 Index actual_k = m_trans ? k : m_length-k-1; 00342 dst.bottomRows(rows()-m_shift-actual_k) 00343 .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data()); 00344 } 00345 } 00346 } 00347 00355 template<typename OtherDerived> 00356 typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const 00357 { 00358 typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type 00359 res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>()); 00360 applyThisOnTheLeft(res); 00361 return res; 00362 } 00363 00364 template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl; 00365 00375 HouseholderSequence& setLength(Index length) 00376 { 00377 m_length = length; 00378 return *this; 00379 } 00380 00392 HouseholderSequence& setShift(Index shift) 00393 { 00394 m_shift = shift; 00395 return *this; 00396 } 00397 00398 Index length() const { return m_length; } 00399 Index shift() const { return m_shift; } 00401 /* Necessary for .adjoint() and .conjugate() */ 00402 template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence; 00403 00404 protected: 00405 00414 HouseholderSequence& setTrans(bool trans) 00415 { 00416 m_trans = trans; 00417 return *this; 00418 } 00419 00420 bool trans() const { return m_trans; } 00422 typename VectorsType::Nested m_vectors; 00423 typename CoeffsType::Nested m_coeffs; 00424 bool m_trans; 00425 Index m_length; 00426 Index m_shift; 00427 }; 00428 00437 template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side> 00438 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence<VectorsType,CoeffsType,Side>& h) 00439 { 00440 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type 00441 res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>()); 00442 h.applyThisOnTheRight(res); 00443 return res; 00444 } 00445 00450 template<typename VectorsType, typename CoeffsType> 00451 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h) 00452 { 00453 return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h); 00454 } 00455 00462 template<typename VectorsType, typename CoeffsType> 00463 HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h) 00464 { 00465 return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h); 00466 } 00467 00468 } // end namespace Eigen 00469 00470 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H