Eigen  3.3.3
HouseholderSequence.h
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
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends