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 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H
00027 #define EIGEN_HOUSEHOLDER_SEQUENCE_H
00028 
00029 namespace Eigen { 
00030 
00072 namespace internal {
00073 
00074 template<typename VectorsType, typename CoeffsType, int Side>
00075 struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
00076 {
00077   typedef typename VectorsType::Scalar Scalar;
00078   typedef typename VectorsType::Index Index;
00079   typedef typename VectorsType::StorageKind StorageKind;
00080   enum {
00081     RowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::RowsAtCompileTime
00082                                         : traits<VectorsType>::ColsAtCompileTime,
00083     ColsAtCompileTime = RowsAtCompileTime,
00084     MaxRowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::MaxRowsAtCompileTime
00085                                            : traits<VectorsType>::MaxColsAtCompileTime,
00086     MaxColsAtCompileTime = MaxRowsAtCompileTime,
00087     Flags = 0
00088   };
00089 };
00090 
00091 template<typename VectorsType, typename CoeffsType, int Side>
00092 struct hseq_side_dependent_impl
00093 {
00094   typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType;
00095   typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType;
00096   typedef typename VectorsType::Index Index;
00097   static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
00098   {
00099     Index start = k+1+h.m_shift;
00100     return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
00101   }
00102 };
00103 
00104 template<typename VectorsType, typename CoeffsType>
00105 struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight>
00106 {
00107   typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType;
00108   typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType;
00109   typedef typename VectorsType::Index Index;
00110   static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
00111   {
00112     Index start = k+1+h.m_shift;
00113     return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose();
00114   }
00115 };
00116 
00117 template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type
00118 {
00119   typedef typename scalar_product_traits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
00120     ResultScalar;
00121   typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
00122                  0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
00123 };
00124 
00125 } // end namespace internal
00126 
00127 template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence
00128   : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
00129 {
00130     enum {
00131       RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
00132       ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
00133       MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
00134       MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
00135     };
00136     typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
00137     typedef typename VectorsType::Index Index;
00138 
00139     typedef typename internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType
00140             EssentialVectorType;
00141 
00142   public:
00143 
00144     typedef HouseholderSequence<
00145       VectorsType,
00146       typename internal::conditional<NumTraits<Scalar>::IsComplex,
00147         typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
00148         CoeffsType>::type,
00149       Side
00150     > ConjugateReturnType;
00151 
00169     HouseholderSequence(const VectorsType& v, const CoeffsType& h)
00170       : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()),
00171         m_shift(0)
00172     {
00173     }
00174 
00176     HouseholderSequence(const HouseholderSequence& other)
00177       : m_vectors(other.m_vectors),
00178         m_coeffs(other.m_coeffs),
00179         m_trans(other.m_trans),
00180         m_length(other.m_length),
00181         m_shift(other.m_shift)
00182     {
00183     }
00184 
00189     Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
00190 
00195     Index cols() const { return rows(); }
00196 
00211     const EssentialVectorType essentialVector(Index k) const
00212     {
00213       eigen_assert(k >= 0 && k < m_length);
00214       return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
00215     }
00216 
00218     HouseholderSequence transpose() const
00219     {
00220       return HouseholderSequence(*this).setTrans(!m_trans);
00221     }
00222 
00224     ConjugateReturnType conjugate() const
00225     {
00226       return ConjugateReturnType(m_vectors, m_coeffs.conjugate())
00227              .setTrans(m_trans)
00228              .setLength(m_length)
00229              .setShift(m_shift);
00230     }
00231 
00233     ConjugateReturnType adjoint() const
00234     {
00235       return conjugate().setTrans(!m_trans);
00236     }
00237 
00239     ConjugateReturnType inverse() const { return adjoint(); }
00240 
00242     template<typename DestType> inline void evalTo(DestType& dst) const
00243     {
00244       Matrix<Scalar, DestType::RowsAtCompileTime, 1,
00245              AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows());
00246       evalTo(dst, workspace);
00247     }
00248 
00250     template<typename Dest, typename Workspace>
00251     void evalTo(Dest& dst, Workspace& workspace) const
00252     {
00253       workspace.resize(rows());
00254       Index vecs = m_length;
00255       if(    internal::is_same<typename internal::remove_all<VectorsType>::type,Dest>::value
00256           && internal::extract_data(dst) == internal::extract_data(m_vectors))
00257       {
00258         // in-place
00259         dst.diagonal().setOnes();
00260         dst.template triangularView<StrictlyUpper>().setZero();
00261         for(Index k = vecs-1; k >= 0; --k)
00262         {
00263           Index cornerSize = rows() - k - m_shift;
00264           if(m_trans)
00265             dst.bottomRightCorner(cornerSize, cornerSize)
00266                .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
00267           else
00268             dst.bottomRightCorner(cornerSize, cornerSize)
00269                .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
00270 
00271           // clear the off diagonal vector
00272           dst.col(k).tail(rows()-k-1).setZero();
00273         }
00274         // clear the remaining columns if needed
00275         for(Index k = 0; k<cols()-vecs ; ++k)
00276           dst.col(k).tail(rows()-k-1).setZero();
00277       }
00278       else
00279       {
00280         dst.setIdentity(rows(), rows());
00281         for(Index k = vecs-1; k >= 0; --k)
00282         {
00283           Index cornerSize = rows() - k - m_shift;
00284           if(m_trans)
00285             dst.bottomRightCorner(cornerSize, cornerSize)
00286                .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
00287           else
00288             dst.bottomRightCorner(cornerSize, cornerSize)
00289                .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
00290         }
00291       }
00292     }
00293 
00295     template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
00296     {
00297       Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
00298       applyThisOnTheRight(dst, workspace);
00299     }
00300 
00302     template<typename Dest, typename Workspace>
00303     inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const
00304     {
00305       workspace.resize(dst.rows());
00306       for(Index k = 0; k < m_length; ++k)
00307       {
00308         Index actual_k = m_trans ? m_length-k-1 : k;
00309         dst.rightCols(rows()-m_shift-actual_k)
00310            .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
00311       }
00312     }
00313 
00315     template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
00316     {
00317       Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace(dst.cols());
00318       applyThisOnTheLeft(dst, workspace);
00319     }
00320 
00322     template<typename Dest, typename Workspace>
00323     inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const
00324     {
00325       workspace.resize(dst.cols());
00326       for(Index k = 0; k < m_length; ++k)
00327       {
00328         Index actual_k = m_trans ? k : m_length-k-1;
00329         dst.bottomRows(rows()-m_shift-actual_k)
00330            .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
00331       }
00332     }
00333 
00341     template<typename OtherDerived>
00342     typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const
00343     {
00344       typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type
00345         res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
00346       applyThisOnTheLeft(res);
00347       return res;
00348     }
00349 
00350     template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
00351 
00361     HouseholderSequence& setLength(Index length)
00362     {
00363       m_length = length;
00364       return *this;
00365     }
00366 
00378     HouseholderSequence& setShift(Index shift)
00379     {
00380       m_shift = shift;
00381       return *this;
00382     }
00383 
00384     Index length() const { return m_length; }  
00385     Index shift() const { return m_shift; }    
00387     /* Necessary for .adjoint() and .conjugate() */
00388     template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence;
00389 
00390   protected:
00391 
00400     HouseholderSequence& setTrans(bool trans)
00401     {
00402       m_trans = trans;
00403       return *this;
00404     }
00405 
00406     bool trans() const { return m_trans; }     
00408     typename VectorsType::Nested m_vectors;
00409     typename CoeffsType::Nested m_coeffs;
00410     bool m_trans;
00411     Index m_length;
00412     Index m_shift;
00413 };
00414 
00423 template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
00424 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence<VectorsType,CoeffsType,Side>& h)
00425 {
00426   typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type
00427     res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
00428   h.applyThisOnTheRight(res);
00429   return res;
00430 }
00431 
00436 template<typename VectorsType, typename CoeffsType>
00437 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
00438 {
00439   return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h);
00440 }
00441 
00448 template<typename VectorsType, typename CoeffsType>
00449 HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h)
00450 {
00451   return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h);
00452 }
00453 
00454 } // end namespace Eigen
00455 
00456 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H