HouseholderSequence.h
Go to the documentation of this file.
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


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:37:24