$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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::Index Index; 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 template<typename VectorsType, typename CoeffsType, int Side> 00077 struct hseq_side_dependent_impl 00078 { 00079 typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType; 00080 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType; 00081 typedef typename VectorsType::Index Index; 00082 static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k) 00083 { 00084 Index start = k+1+h.m_shift; 00085 return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1); 00086 } 00087 }; 00088 00089 template<typename VectorsType, typename CoeffsType> 00090 struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight> 00091 { 00092 typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType; 00093 typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> 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,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose(); 00099 } 00100 }; 00101 00102 template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type 00103 { 00104 typedef typename scalar_product_traits<OtherScalarType, typename MatrixType::Scalar>::ReturnType 00105 ResultScalar; 00106 typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime, 00107 0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type; 00108 }; 00109 00110 } // end namespace internal 00111 00112 template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence 00113 : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> > 00114 { 00115 typedef typename internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType EssentialVectorType; 00116 00117 public: 00118 enum { 00119 RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime, 00120 ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime, 00121 MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime, 00122 MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime 00123 }; 00124 typedef typename internal::traits<HouseholderSequence>::Scalar Scalar; 00125 typedef typename VectorsType::Index Index; 00126 00127 typedef HouseholderSequence< 00128 typename internal::conditional<NumTraits<Scalar>::IsComplex, 00129 typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type, 00130 VectorsType>::type, 00131 typename internal::conditional<NumTraits<Scalar>::IsComplex, 00132 typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type, 00133 CoeffsType>::type, 00134 Side 00135 > ConjugateReturnType; 00136 00154 HouseholderSequence(const VectorsType& v, const CoeffsType& h) 00155 : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()), 00156 m_shift(0) 00157 { 00158 } 00159 00161 HouseholderSequence(const HouseholderSequence& other) 00162 : m_vectors(other.m_vectors), 00163 m_coeffs(other.m_coeffs), 00164 m_trans(other.m_trans), 00165 m_length(other.m_length), 00166 m_shift(other.m_shift) 00167 { 00168 } 00169 00174 Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); } 00175 00180 Index cols() const { return rows(); } 00181 00196 const EssentialVectorType essentialVector(Index k) const 00197 { 00198 eigen_assert(k >= 0 && k < m_length); 00199 return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k); 00200 } 00201 00203 HouseholderSequence transpose() const 00204 { 00205 return HouseholderSequence(*this).setTrans(!m_trans); 00206 } 00207 00209 ConjugateReturnType conjugate() const 00210 { 00211 return ConjugateReturnType(m_vectors.conjugate(), m_coeffs.conjugate()) 00212 .setTrans(m_trans) 00213 .setLength(m_length) 00214 .setShift(m_shift); 00215 } 00216 00218 ConjugateReturnType adjoint() const 00219 { 00220 return conjugate().setTrans(!m_trans); 00221 } 00222 00224 ConjugateReturnType inverse() const { return adjoint(); } 00225 00227 template<typename DestType> inline void evalTo(DestType& dst) const 00228 { 00229 Matrix<Scalar, DestType::RowsAtCompileTime, 1, 00230 AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows()); 00231 evalTo(dst, workspace); 00232 } 00233 00235 template<typename Dest, typename Workspace> 00236 void evalTo(Dest& dst, Workspace& workspace) const 00237 { 00238 workspace.resize(rows()); 00239 Index vecs = m_length; 00240 if( internal::is_same<typename internal::remove_all<VectorsType>::type,Dest>::value 00241 && internal::extract_data(dst) == internal::extract_data(m_vectors)) 00242 { 00243 // in-place 00244 dst.diagonal().setOnes(); 00245 dst.template triangularView<StrictlyUpper>().setZero(); 00246 for(Index k = vecs-1; k >= 0; --k) 00247 { 00248 Index cornerSize = rows() - k - m_shift; 00249 if(m_trans) 00250 dst.bottomRightCorner(cornerSize, cornerSize) 00251 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data()); 00252 else 00253 dst.bottomRightCorner(cornerSize, cornerSize) 00254 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data()); 00255 00256 // clear the off diagonal vector 00257 dst.col(k).tail(rows()-k-1).setZero(); 00258 } 00259 // clear the remaining columns if needed 00260 for(Index k = 0; k<cols()-vecs ; ++k) 00261 dst.col(k).tail(rows()-k-1).setZero(); 00262 } 00263 else 00264 { 00265 dst.setIdentity(rows(), rows()); 00266 for(Index k = vecs-1; k >= 0; --k) 00267 { 00268 Index cornerSize = rows() - k - m_shift; 00269 if(m_trans) 00270 dst.bottomRightCorner(cornerSize, cornerSize) 00271 .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0)); 00272 else 00273 dst.bottomRightCorner(cornerSize, cornerSize) 00274 .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0)); 00275 } 00276 } 00277 } 00278 00280 template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const 00281 { 00282 Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows()); 00283 applyThisOnTheRight(dst, workspace); 00284 } 00285 00287 template<typename Dest, typename Workspace> 00288 inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const 00289 { 00290 workspace.resize(dst.rows()); 00291 for(Index k = 0; k < m_length; ++k) 00292 { 00293 Index actual_k = m_trans ? m_length-k-1 : k; 00294 dst.rightCols(rows()-m_shift-actual_k) 00295 .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data()); 00296 } 00297 } 00298 00300 template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const 00301 { 00302 Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace(dst.cols()); 00303 applyThisOnTheLeft(dst, workspace); 00304 } 00305 00307 template<typename Dest, typename Workspace> 00308 inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const 00309 { 00310 workspace.resize(dst.cols()); 00311 for(Index k = 0; k < m_length; ++k) 00312 { 00313 Index actual_k = m_trans ? k : m_length-k-1; 00314 dst.bottomRows(rows()-m_shift-actual_k) 00315 .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data()); 00316 } 00317 } 00318 00326 template<typename OtherDerived> 00327 typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const 00328 { 00329 typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type 00330 res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>()); 00331 applyThisOnTheLeft(res); 00332 return res; 00333 } 00334 00335 template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl; 00336 00346 HouseholderSequence& setLength(Index length) 00347 { 00348 m_length = length; 00349 return *this; 00350 } 00351 00363 HouseholderSequence& setShift(Index shift) 00364 { 00365 m_shift = shift; 00366 return *this; 00367 } 00368 00369 Index length() const { return m_length; } 00370 Index shift() const { return m_shift; } 00372 /* Necessary for .adjoint() and .conjugate() */ 00373 template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence; 00374 00375 protected: 00376 00385 HouseholderSequence& setTrans(bool trans) 00386 { 00387 m_trans = trans; 00388 return *this; 00389 } 00390 00391 bool trans() const { return m_trans; } 00393 typename VectorsType::Nested m_vectors; 00394 typename CoeffsType::Nested m_coeffs; 00395 bool m_trans; 00396 Index m_length; 00397 Index m_shift; 00398 }; 00399 00408 template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side> 00409 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence<VectorsType,CoeffsType,Side>& h) 00410 { 00411 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type 00412 res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>()); 00413 h.applyThisOnTheRight(res); 00414 return res; 00415 } 00416 00421 template<typename VectorsType, typename CoeffsType> 00422 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h) 00423 { 00424 return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h); 00425 } 00426 00433 template<typename VectorsType, typename CoeffsType> 00434 HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h) 00435 { 00436 return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h); 00437 } 00438 00439 } // end namespace Eigen 00440 00441 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H