Home | History | Annotate | Download | only in Core
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2010-2011 Gael Guennebaud <gael.guennebaud (at) inria.fr>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #ifndef EIGEN_TRANSPOSITIONS_H
     11 #define EIGEN_TRANSPOSITIONS_H
     12 
     13 namespace Eigen {
     14 
     15 /** \class Transpositions
     16   * \ingroup Core_Module
     17   *
     18   * \brief Represents a sequence of transpositions (row/column interchange)
     19   *
     20   * \param SizeAtCompileTime the number of transpositions, or Dynamic
     21   * \param MaxSizeAtCompileTime the maximum number of transpositions, or Dynamic. This optional parameter defaults to SizeAtCompileTime. Most of the time, you should not have to specify it.
     22   *
     23   * This class represents a permutation transformation as a sequence of \em n transpositions
     24   * \f$[T_{n-1} \ldots T_{i} \ldots T_{0}]\f$. It is internally stored as a vector of integers \c indices.
     25   * Each transposition \f$ T_{i} \f$ applied on the left of a matrix (\f$ T_{i} M\f$) interchanges
     26   * the rows \c i and \c indices[i] of the matrix \c M.
     27   * A transposition applied on the right (e.g., \f$ M T_{i}\f$) yields a column interchange.
     28   *
     29   * Compared to the class PermutationMatrix, such a sequence of transpositions is what is
     30   * computed during a decomposition with pivoting, and it is faster when applying the permutation in-place.
     31   *
     32   * To apply a sequence of transpositions to a matrix, simply use the operator * as in the following example:
     33   * \code
     34   * Transpositions tr;
     35   * MatrixXf mat;
     36   * mat = tr * mat;
     37   * \endcode
     38   * In this example, we detect that the matrix appears on both side, and so the transpositions
     39   * are applied in-place without any temporary or extra copy.
     40   *
     41   * \sa class PermutationMatrix
     42   */
     43 
     44 namespace internal {
     45 template<typename TranspositionType, typename MatrixType, int Side, bool Transposed=false> struct transposition_matrix_product_retval;
     46 }
     47 
     48 template<typename Derived>
     49 class TranspositionsBase
     50 {
     51     typedef internal::traits<Derived> Traits;
     52 
     53   public:
     54 
     55     typedef typename Traits::IndicesType IndicesType;
     56     typedef typename IndicesType::Scalar Index;
     57 
     58     Derived& derived() { return *static_cast<Derived*>(this); }
     59     const Derived& derived() const { return *static_cast<const Derived*>(this); }
     60 
     61     /** Copies the \a other transpositions into \c *this */
     62     template<typename OtherDerived>
     63     Derived& operator=(const TranspositionsBase<OtherDerived>& other)
     64     {
     65       indices() = other.indices();
     66       return derived();
     67     }
     68 
     69     #ifndef EIGEN_PARSED_BY_DOXYGEN
     70     /** This is a special case of the templated operator=. Its purpose is to
     71       * prevent a default operator= from hiding the templated operator=.
     72       */
     73     Derived& operator=(const TranspositionsBase& other)
     74     {
     75       indices() = other.indices();
     76       return derived();
     77     }
     78     #endif
     79 
     80     /** \returns the number of transpositions */
     81     inline Index size() const { return indices().size(); }
     82 
     83     /** Direct access to the underlying index vector */
     84     inline const Index& coeff(Index i) const { return indices().coeff(i); }
     85     /** Direct access to the underlying index vector */
     86     inline Index& coeffRef(Index i) { return indices().coeffRef(i); }
     87     /** Direct access to the underlying index vector */
     88     inline const Index& operator()(Index i) const { return indices()(i); }
     89     /** Direct access to the underlying index vector */
     90     inline Index& operator()(Index i) { return indices()(i); }
     91     /** Direct access to the underlying index vector */
     92     inline const Index& operator[](Index i) const { return indices()(i); }
     93     /** Direct access to the underlying index vector */
     94     inline Index& operator[](Index i) { return indices()(i); }
     95 
     96     /** const version of indices(). */
     97     const IndicesType& indices() const { return derived().indices(); }
     98     /** \returns a reference to the stored array representing the transpositions. */
     99     IndicesType& indices() { return derived().indices(); }
    100 
    101     /** Resizes to given size. */
    102     inline void resize(int newSize)
    103     {
    104       indices().resize(newSize);
    105     }
    106 
    107     /** Sets \c *this to represents an identity transformation */
    108     void setIdentity()
    109     {
    110       for(int i = 0; i < indices().size(); ++i)
    111         coeffRef(i) = i;
    112     }
    113 
    114     // FIXME: do we want such methods ?
    115     // might be usefull when the target matrix expression is complex, e.g.:
    116     // object.matrix().block(..,..,..,..) = trans * object.matrix().block(..,..,..,..);
    117     /*
    118     template<typename MatrixType>
    119     void applyForwardToRows(MatrixType& mat) const
    120     {
    121       for(Index k=0 ; k<size() ; ++k)
    122         if(m_indices(k)!=k)
    123           mat.row(k).swap(mat.row(m_indices(k)));
    124     }
    125 
    126     template<typename MatrixType>
    127     void applyBackwardToRows(MatrixType& mat) const
    128     {
    129       for(Index k=size()-1 ; k>=0 ; --k)
    130         if(m_indices(k)!=k)
    131           mat.row(k).swap(mat.row(m_indices(k)));
    132     }
    133     */
    134 
    135     /** \returns the inverse transformation */
    136     inline Transpose<TranspositionsBase> inverse() const
    137     { return Transpose<TranspositionsBase>(derived()); }
    138 
    139     /** \returns the tranpose transformation */
    140     inline Transpose<TranspositionsBase> transpose() const
    141     { return Transpose<TranspositionsBase>(derived()); }
    142 
    143   protected:
    144 };
    145 
    146 namespace internal {
    147 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
    148 struct traits<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType> >
    149 {
    150   typedef IndexType Index;
    151   typedef Matrix<Index, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
    152 };
    153 }
    154 
    155 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
    156 class Transpositions : public TranspositionsBase<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType> >
    157 {
    158     typedef internal::traits<Transpositions> Traits;
    159   public:
    160 
    161     typedef TranspositionsBase<Transpositions> Base;
    162     typedef typename Traits::IndicesType IndicesType;
    163     typedef typename IndicesType::Scalar Index;
    164 
    165     inline Transpositions() {}
    166 
    167     /** Copy constructor. */
    168     template<typename OtherDerived>
    169     inline Transpositions(const TranspositionsBase<OtherDerived>& other)
    170       : m_indices(other.indices()) {}
    171 
    172     #ifndef EIGEN_PARSED_BY_DOXYGEN
    173     /** Standard copy constructor. Defined only to prevent a default copy constructor
    174       * from hiding the other templated constructor */
    175     inline Transpositions(const Transpositions& other) : m_indices(other.indices()) {}
    176     #endif
    177 
    178     /** Generic constructor from expression of the transposition indices. */
    179     template<typename Other>
    180     explicit inline Transpositions(const MatrixBase<Other>& a_indices) : m_indices(a_indices)
    181     {}
    182 
    183     /** Copies the \a other transpositions into \c *this */
    184     template<typename OtherDerived>
    185     Transpositions& operator=(const TranspositionsBase<OtherDerived>& other)
    186     {
    187       return Base::operator=(other);
    188     }
    189 
    190     #ifndef EIGEN_PARSED_BY_DOXYGEN
    191     /** This is a special case of the templated operator=. Its purpose is to
    192       * prevent a default operator= from hiding the templated operator=.
    193       */
    194     Transpositions& operator=(const Transpositions& other)
    195     {
    196       m_indices = other.m_indices;
    197       return *this;
    198     }
    199     #endif
    200 
    201     /** Constructs an uninitialized permutation matrix of given size.
    202       */
    203     inline Transpositions(Index size) : m_indices(size)
    204     {}
    205 
    206     /** const version of indices(). */
    207     const IndicesType& indices() const { return m_indices; }
    208     /** \returns a reference to the stored array representing the transpositions. */
    209     IndicesType& indices() { return m_indices; }
    210 
    211   protected:
    212 
    213     IndicesType m_indices;
    214 };
    215 
    216 
    217 namespace internal {
    218 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
    219 struct traits<Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType>,_PacketAccess> >
    220 {
    221   typedef IndexType Index;
    222   typedef Map<const Matrix<Index,SizeAtCompileTime,1,0,MaxSizeAtCompileTime,1>, _PacketAccess> IndicesType;
    223 };
    224 }
    225 
    226 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int PacketAccess>
    227 class Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType>,PacketAccess>
    228  : public TranspositionsBase<Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType>,PacketAccess> >
    229 {
    230     typedef internal::traits<Map> Traits;
    231   public:
    232 
    233     typedef TranspositionsBase<Map> Base;
    234     typedef typename Traits::IndicesType IndicesType;
    235     typedef typename IndicesType::Scalar Index;
    236 
    237     inline Map(const Index* indicesPtr)
    238       : m_indices(indicesPtr)
    239     {}
    240 
    241     inline Map(const Index* indicesPtr, Index size)
    242       : m_indices(indicesPtr,size)
    243     {}
    244 
    245     /** Copies the \a other transpositions into \c *this */
    246     template<typename OtherDerived>
    247     Map& operator=(const TranspositionsBase<OtherDerived>& other)
    248     {
    249       return Base::operator=(other);
    250     }
    251 
    252     #ifndef EIGEN_PARSED_BY_DOXYGEN
    253     /** This is a special case of the templated operator=. Its purpose is to
    254       * prevent a default operator= from hiding the templated operator=.
    255       */
    256     Map& operator=(const Map& other)
    257     {
    258       m_indices = other.m_indices;
    259       return *this;
    260     }
    261     #endif
    262 
    263     /** const version of indices(). */
    264     const IndicesType& indices() const { return m_indices; }
    265 
    266     /** \returns a reference to the stored array representing the transpositions. */
    267     IndicesType& indices() { return m_indices; }
    268 
    269   protected:
    270 
    271     IndicesType m_indices;
    272 };
    273 
    274 namespace internal {
    275 template<typename _IndicesType>
    276 struct traits<TranspositionsWrapper<_IndicesType> >
    277 {
    278   typedef typename _IndicesType::Scalar Index;
    279   typedef _IndicesType IndicesType;
    280 };
    281 }
    282 
    283 template<typename _IndicesType>
    284 class TranspositionsWrapper
    285  : public TranspositionsBase<TranspositionsWrapper<_IndicesType> >
    286 {
    287     typedef internal::traits<TranspositionsWrapper> Traits;
    288   public:
    289 
    290     typedef TranspositionsBase<TranspositionsWrapper> Base;
    291     typedef typename Traits::IndicesType IndicesType;
    292     typedef typename IndicesType::Scalar Index;
    293 
    294     inline TranspositionsWrapper(IndicesType& a_indices)
    295       : m_indices(a_indices)
    296     {}
    297 
    298     /** Copies the \a other transpositions into \c *this */
    299     template<typename OtherDerived>
    300     TranspositionsWrapper& operator=(const TranspositionsBase<OtherDerived>& other)
    301     {
    302       return Base::operator=(other);
    303     }
    304 
    305     #ifndef EIGEN_PARSED_BY_DOXYGEN
    306     /** This is a special case of the templated operator=. Its purpose is to
    307       * prevent a default operator= from hiding the templated operator=.
    308       */
    309     TranspositionsWrapper& operator=(const TranspositionsWrapper& other)
    310     {
    311       m_indices = other.m_indices;
    312       return *this;
    313     }
    314     #endif
    315 
    316     /** const version of indices(). */
    317     const IndicesType& indices() const { return m_indices; }
    318 
    319     /** \returns a reference to the stored array representing the transpositions. */
    320     IndicesType& indices() { return m_indices; }
    321 
    322   protected:
    323 
    324     const typename IndicesType::Nested m_indices;
    325 };
    326 
    327 /** \returns the \a matrix with the \a transpositions applied to the columns.
    328   */
    329 template<typename Derived, typename TranspositionsDerived>
    330 inline const internal::transposition_matrix_product_retval<TranspositionsDerived, Derived, OnTheRight>
    331 operator*(const MatrixBase<Derived>& matrix,
    332           const TranspositionsBase<TranspositionsDerived> &transpositions)
    333 {
    334   return internal::transposition_matrix_product_retval
    335            <TranspositionsDerived, Derived, OnTheRight>
    336            (transpositions.derived(), matrix.derived());
    337 }
    338 
    339 /** \returns the \a matrix with the \a transpositions applied to the rows.
    340   */
    341 template<typename Derived, typename TranspositionDerived>
    342 inline const internal::transposition_matrix_product_retval
    343                <TranspositionDerived, Derived, OnTheLeft>
    344 operator*(const TranspositionsBase<TranspositionDerived> &transpositions,
    345           const MatrixBase<Derived>& matrix)
    346 {
    347   return internal::transposition_matrix_product_retval
    348            <TranspositionDerived, Derived, OnTheLeft>
    349            (transpositions.derived(), matrix.derived());
    350 }
    351 
    352 namespace internal {
    353 
    354 template<typename TranspositionType, typename MatrixType, int Side, bool Transposed>
    355 struct traits<transposition_matrix_product_retval<TranspositionType, MatrixType, Side, Transposed> >
    356 {
    357   typedef typename MatrixType::PlainObject ReturnType;
    358 };
    359 
    360 template<typename TranspositionType, typename MatrixType, int Side, bool Transposed>
    361 struct transposition_matrix_product_retval
    362  : public ReturnByValue<transposition_matrix_product_retval<TranspositionType, MatrixType, Side, Transposed> >
    363 {
    364     typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
    365     typedef typename TranspositionType::Index Index;
    366 
    367     transposition_matrix_product_retval(const TranspositionType& tr, const MatrixType& matrix)
    368       : m_transpositions(tr), m_matrix(matrix)
    369     {}
    370 
    371     inline int rows() const { return m_matrix.rows(); }
    372     inline int cols() const { return m_matrix.cols(); }
    373 
    374     template<typename Dest> inline void evalTo(Dest& dst) const
    375     {
    376       const int size = m_transpositions.size();
    377       Index j = 0;
    378 
    379       if(!(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix)))
    380         dst = m_matrix;
    381 
    382       for(int k=(Transposed?size-1:0) ; Transposed?k>=0:k<size ; Transposed?--k:++k)
    383         if((j=m_transpositions.coeff(k))!=k)
    384         {
    385           if(Side==OnTheLeft)
    386             dst.row(k).swap(dst.row(j));
    387           else if(Side==OnTheRight)
    388             dst.col(k).swap(dst.col(j));
    389         }
    390     }
    391 
    392   protected:
    393     const TranspositionType& m_transpositions;
    394     typename MatrixType::Nested m_matrix;
    395 };
    396 
    397 } // end namespace internal
    398 
    399 /* Template partial specialization for transposed/inverse transpositions */
    400 
    401 template<typename TranspositionsDerived>
    402 class Transpose<TranspositionsBase<TranspositionsDerived> >
    403 {
    404     typedef TranspositionsDerived TranspositionType;
    405     typedef typename TranspositionType::IndicesType IndicesType;
    406   public:
    407 
    408     Transpose(const TranspositionType& t) : m_transpositions(t) {}
    409 
    410     inline int size() const { return m_transpositions.size(); }
    411 
    412     /** \returns the \a matrix with the inverse transpositions applied to the columns.
    413       */
    414     template<typename Derived> friend
    415     inline const internal::transposition_matrix_product_retval<TranspositionType, Derived, OnTheRight, true>
    416     operator*(const MatrixBase<Derived>& matrix, const Transpose& trt)
    417     {
    418       return internal::transposition_matrix_product_retval<TranspositionType, Derived, OnTheRight, true>(trt.m_transpositions, matrix.derived());
    419     }
    420 
    421     /** \returns the \a matrix with the inverse transpositions applied to the rows.
    422       */
    423     template<typename Derived>
    424     inline const internal::transposition_matrix_product_retval<TranspositionType, Derived, OnTheLeft, true>
    425     operator*(const MatrixBase<Derived>& matrix) const
    426     {
    427       return internal::transposition_matrix_product_retval<TranspositionType, Derived, OnTheLeft, true>(m_transpositions, matrix.derived());
    428     }
    429 
    430   protected:
    431     const TranspositionType& m_transpositions;
    432 };
    433 
    434 } // end namespace Eigen
    435 
    436 #endif // EIGEN_TRANSPOSITIONS_H
    437