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) 2009 Benoit Jacob <jacob.benoit.1 (at) gmail.com>
      5 // Copyright (C) 2009-2011 Gael Guennebaud <gael.guennebaud (at) inria.fr>
      6 //
      7 // This Source Code Form is subject to the terms of the Mozilla
      8 // Public License v. 2.0. If a copy of the MPL was not distributed
      9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
     10 
     11 #ifndef EIGEN_PERMUTATIONMATRIX_H
     12 #define EIGEN_PERMUTATIONMATRIX_H
     13 
     14 namespace Eigen {
     15 
     16 template<int RowCol,typename IndicesType,typename MatrixType, typename StorageKind> class PermutedImpl;
     17 
     18 /** \class PermutationBase
     19   * \ingroup Core_Module
     20   *
     21   * \brief Base class for permutations
     22   *
     23   * \param Derived the derived class
     24   *
     25   * This class is the base class for all expressions representing a permutation matrix,
     26   * internally stored as a vector of integers.
     27   * The convention followed here is that if \f$ \sigma \f$ is a permutation, the corresponding permutation matrix
     28   * \f$ P_\sigma \f$ is such that if \f$ (e_1,\ldots,e_p) \f$ is the canonical basis, we have:
     29   *  \f[ P_\sigma(e_i) = e_{\sigma(i)}. \f]
     30   * This convention ensures that for any two permutations \f$ \sigma, \tau \f$, we have:
     31   *  \f[ P_{\sigma\circ\tau} = P_\sigma P_\tau. \f]
     32   *
     33   * Permutation matrices are square and invertible.
     34   *
     35   * Notice that in addition to the member functions and operators listed here, there also are non-member
     36   * operator* to multiply any kind of permutation object with any kind of matrix expression (MatrixBase)
     37   * on either side.
     38   *
     39   * \sa class PermutationMatrix, class PermutationWrapper
     40   */
     41 
     42 namespace internal {
     43 
     44 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
     45 struct permut_matrix_product_retval;
     46 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
     47 struct permut_sparsematrix_product_retval;
     48 enum PermPermProduct_t {PermPermProduct};
     49 
     50 } // end namespace internal
     51 
     52 template<typename Derived>
     53 class PermutationBase : public EigenBase<Derived>
     54 {
     55     typedef internal::traits<Derived> Traits;
     56     typedef EigenBase<Derived> Base;
     57   public:
     58 
     59     #ifndef EIGEN_PARSED_BY_DOXYGEN
     60     typedef typename Traits::IndicesType IndicesType;
     61     enum {
     62       Flags = Traits::Flags,
     63       CoeffReadCost = Traits::CoeffReadCost,
     64       RowsAtCompileTime = Traits::RowsAtCompileTime,
     65       ColsAtCompileTime = Traits::ColsAtCompileTime,
     66       MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
     67       MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
     68     };
     69     typedef typename Traits::Scalar Scalar;
     70     typedef typename Traits::Index Index;
     71     typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime,0,MaxRowsAtCompileTime,MaxColsAtCompileTime>
     72             DenseMatrixType;
     73     typedef PermutationMatrix<IndicesType::SizeAtCompileTime,IndicesType::MaxSizeAtCompileTime,Index>
     74             PlainPermutationType;
     75     using Base::derived;
     76     #endif
     77 
     78     /** Copies the other permutation into *this */
     79     template<typename OtherDerived>
     80     Derived& operator=(const PermutationBase<OtherDerived>& other)
     81     {
     82       indices() = other.indices();
     83       return derived();
     84     }
     85 
     86     /** Assignment from the Transpositions \a tr */
     87     template<typename OtherDerived>
     88     Derived& operator=(const TranspositionsBase<OtherDerived>& tr)
     89     {
     90       setIdentity(tr.size());
     91       for(Index k=size()-1; k>=0; --k)
     92         applyTranspositionOnTheRight(k,tr.coeff(k));
     93       return derived();
     94     }
     95 
     96     #ifndef EIGEN_PARSED_BY_DOXYGEN
     97     /** This is a special case of the templated operator=. Its purpose is to
     98       * prevent a default operator= from hiding the templated operator=.
     99       */
    100     Derived& operator=(const PermutationBase& other)
    101     {
    102       indices() = other.indices();
    103       return derived();
    104     }
    105     #endif
    106 
    107     /** \returns the number of rows */
    108     inline Index rows() const { return Index(indices().size()); }
    109 
    110     /** \returns the number of columns */
    111     inline Index cols() const { return Index(indices().size()); }
    112 
    113     /** \returns the size of a side of the respective square matrix, i.e., the number of indices */
    114     inline Index size() const { return Index(indices().size()); }
    115 
    116     #ifndef EIGEN_PARSED_BY_DOXYGEN
    117     template<typename DenseDerived>
    118     void evalTo(MatrixBase<DenseDerived>& other) const
    119     {
    120       other.setZero();
    121       for (int i=0; i<rows();++i)
    122         other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
    123     }
    124     #endif
    125 
    126     /** \returns a Matrix object initialized from this permutation matrix. Notice that it
    127       * is inefficient to return this Matrix object by value. For efficiency, favor using
    128       * the Matrix constructor taking EigenBase objects.
    129       */
    130     DenseMatrixType toDenseMatrix() const
    131     {
    132       return derived();
    133     }
    134 
    135     /** const version of indices(). */
    136     const IndicesType& indices() const { return derived().indices(); }
    137     /** \returns a reference to the stored array representing the permutation. */
    138     IndicesType& indices() { return derived().indices(); }
    139 
    140     /** Resizes to given size.
    141       */
    142     inline void resize(Index newSize)
    143     {
    144       indices().resize(newSize);
    145     }
    146 
    147     /** Sets *this to be the identity permutation matrix */
    148     void setIdentity()
    149     {
    150       for(Index i = 0; i < size(); ++i)
    151         indices().coeffRef(i) = i;
    152     }
    153 
    154     /** Sets *this to be the identity permutation matrix of given size.
    155       */
    156     void setIdentity(Index newSize)
    157     {
    158       resize(newSize);
    159       setIdentity();
    160     }
    161 
    162     /** Multiplies *this by the transposition \f$(ij)\f$ on the left.
    163       *
    164       * \returns a reference to *this.
    165       *
    166       * \warning This is much slower than applyTranspositionOnTheRight(int,int):
    167       * this has linear complexity and requires a lot of branching.
    168       *
    169       * \sa applyTranspositionOnTheRight(int,int)
    170       */
    171     Derived& applyTranspositionOnTheLeft(Index i, Index j)
    172     {
    173       eigen_assert(i>=0 && j>=0 && i<size() && j<size());
    174       for(Index k = 0; k < size(); ++k)
    175       {
    176         if(indices().coeff(k) == i) indices().coeffRef(k) = j;
    177         else if(indices().coeff(k) == j) indices().coeffRef(k) = i;
    178       }
    179       return derived();
    180     }
    181 
    182     /** Multiplies *this by the transposition \f$(ij)\f$ on the right.
    183       *
    184       * \returns a reference to *this.
    185       *
    186       * This is a fast operation, it only consists in swapping two indices.
    187       *
    188       * \sa applyTranspositionOnTheLeft(int,int)
    189       */
    190     Derived& applyTranspositionOnTheRight(Index i, Index j)
    191     {
    192       eigen_assert(i>=0 && j>=0 && i<size() && j<size());
    193       std::swap(indices().coeffRef(i), indices().coeffRef(j));
    194       return derived();
    195     }
    196 
    197     /** \returns the inverse permutation matrix.
    198       *
    199       * \note \note_try_to_help_rvo
    200       */
    201     inline Transpose<PermutationBase> inverse() const
    202     { return derived(); }
    203     /** \returns the tranpose permutation matrix.
    204       *
    205       * \note \note_try_to_help_rvo
    206       */
    207     inline Transpose<PermutationBase> transpose() const
    208     { return derived(); }
    209 
    210     /**** multiplication helpers to hopefully get RVO ****/
    211 
    212 
    213 #ifndef EIGEN_PARSED_BY_DOXYGEN
    214   protected:
    215     template<typename OtherDerived>
    216     void assignTranspose(const PermutationBase<OtherDerived>& other)
    217     {
    218       for (int i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i;
    219     }
    220     template<typename Lhs,typename Rhs>
    221     void assignProduct(const Lhs& lhs, const Rhs& rhs)
    222     {
    223       eigen_assert(lhs.cols() == rhs.rows());
    224       for (int i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
    225     }
    226 #endif
    227 
    228   public:
    229 
    230     /** \returns the product permutation matrix.
    231       *
    232       * \note \note_try_to_help_rvo
    233       */
    234     template<typename Other>
    235     inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
    236     { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
    237 
    238     /** \returns the product of a permutation with another inverse permutation.
    239       *
    240       * \note \note_try_to_help_rvo
    241       */
    242     template<typename Other>
    243     inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other) const
    244     { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
    245 
    246     /** \returns the product of an inverse permutation with another permutation.
    247       *
    248       * \note \note_try_to_help_rvo
    249       */
    250     template<typename Other> friend
    251     inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other, const PermutationBase& perm)
    252     { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
    253 
    254   protected:
    255 
    256 };
    257 
    258 /** \class PermutationMatrix
    259   * \ingroup Core_Module
    260   *
    261   * \brief Permutation matrix
    262   *
    263   * \param SizeAtCompileTime the number of rows/cols, or Dynamic
    264   * \param MaxSizeAtCompileTime the maximum number of rows/cols, or Dynamic. This optional parameter defaults to SizeAtCompileTime. Most of the time, you should not have to specify it.
    265   * \param IndexType the interger type of the indices
    266   *
    267   * This class represents a permutation matrix, internally stored as a vector of integers.
    268   *
    269   * \sa class PermutationBase, class PermutationWrapper, class DiagonalMatrix
    270   */
    271 
    272 namespace internal {
    273 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
    274 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
    275  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
    276 {
    277   typedef IndexType Index;
    278   typedef Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
    279 };
    280 }
    281 
    282 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
    283 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
    284 {
    285     typedef PermutationBase<PermutationMatrix> Base;
    286     typedef internal::traits<PermutationMatrix> Traits;
    287   public:
    288 
    289     #ifndef EIGEN_PARSED_BY_DOXYGEN
    290     typedef typename Traits::IndicesType IndicesType;
    291     #endif
    292 
    293     inline PermutationMatrix()
    294     {}
    295 
    296     /** Constructs an uninitialized permutation matrix of given size.
    297       */
    298     inline PermutationMatrix(int size) : m_indices(size)
    299     {}
    300 
    301     /** Copy constructor. */
    302     template<typename OtherDerived>
    303     inline PermutationMatrix(const PermutationBase<OtherDerived>& other)
    304       : m_indices(other.indices()) {}
    305 
    306     #ifndef EIGEN_PARSED_BY_DOXYGEN
    307     /** Standard copy constructor. Defined only to prevent a default copy constructor
    308       * from hiding the other templated constructor */
    309     inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
    310     #endif
    311 
    312     /** Generic constructor from expression of the indices. The indices
    313       * array has the meaning that the permutations sends each integer i to indices[i].
    314       *
    315       * \warning It is your responsibility to check that the indices array that you passes actually
    316       * describes a permutation, i.e., each value between 0 and n-1 occurs exactly once, where n is the
    317       * array's size.
    318       */
    319     template<typename Other>
    320     explicit inline PermutationMatrix(const MatrixBase<Other>& a_indices) : m_indices(a_indices)
    321     {}
    322 
    323     /** Convert the Transpositions \a tr to a permutation matrix */
    324     template<typename Other>
    325     explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
    326       : m_indices(tr.size())
    327     {
    328       *this = tr;
    329     }
    330 
    331     /** Copies the other permutation into *this */
    332     template<typename Other>
    333     PermutationMatrix& operator=(const PermutationBase<Other>& other)
    334     {
    335       m_indices = other.indices();
    336       return *this;
    337     }
    338 
    339     /** Assignment from the Transpositions \a tr */
    340     template<typename Other>
    341     PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
    342     {
    343       return Base::operator=(tr.derived());
    344     }
    345 
    346     #ifndef EIGEN_PARSED_BY_DOXYGEN
    347     /** This is a special case of the templated operator=. Its purpose is to
    348       * prevent a default operator= from hiding the templated operator=.
    349       */
    350     PermutationMatrix& operator=(const PermutationMatrix& other)
    351     {
    352       m_indices = other.m_indices;
    353       return *this;
    354     }
    355     #endif
    356 
    357     /** const version of indices(). */
    358     const IndicesType& indices() const { return m_indices; }
    359     /** \returns a reference to the stored array representing the permutation. */
    360     IndicesType& indices() { return m_indices; }
    361 
    362 
    363     /**** multiplication helpers to hopefully get RVO ****/
    364 
    365 #ifndef EIGEN_PARSED_BY_DOXYGEN
    366     template<typename Other>
    367     PermutationMatrix(const Transpose<PermutationBase<Other> >& other)
    368       : m_indices(other.nestedPermutation().size())
    369     {
    370       for (int i=0; i<m_indices.size();++i) m_indices.coeffRef(other.nestedPermutation().indices().coeff(i)) = i;
    371     }
    372     template<typename Lhs,typename Rhs>
    373     PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
    374       : m_indices(lhs.indices().size())
    375     {
    376       Base::assignProduct(lhs,rhs);
    377     }
    378 #endif
    379 
    380   protected:
    381 
    382     IndicesType m_indices;
    383 };
    384 
    385 
    386 namespace internal {
    387 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
    388 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
    389  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
    390 {
    391   typedef IndexType Index;
    392   typedef Map<const Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType;
    393 };
    394 }
    395 
    396 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
    397 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess>
    398   : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
    399 {
    400     typedef PermutationBase<Map> Base;
    401     typedef internal::traits<Map> Traits;
    402   public:
    403 
    404     #ifndef EIGEN_PARSED_BY_DOXYGEN
    405     typedef typename Traits::IndicesType IndicesType;
    406     typedef typename IndicesType::Scalar Index;
    407     #endif
    408 
    409     inline Map(const Index* indicesPtr)
    410       : m_indices(indicesPtr)
    411     {}
    412 
    413     inline Map(const Index* indicesPtr, Index size)
    414       : m_indices(indicesPtr,size)
    415     {}
    416 
    417     /** Copies the other permutation into *this */
    418     template<typename Other>
    419     Map& operator=(const PermutationBase<Other>& other)
    420     { return Base::operator=(other.derived()); }
    421 
    422     /** Assignment from the Transpositions \a tr */
    423     template<typename Other>
    424     Map& operator=(const TranspositionsBase<Other>& tr)
    425     { return Base::operator=(tr.derived()); }
    426 
    427     #ifndef EIGEN_PARSED_BY_DOXYGEN
    428     /** This is a special case of the templated operator=. Its purpose is to
    429       * prevent a default operator= from hiding the templated operator=.
    430       */
    431     Map& operator=(const Map& other)
    432     {
    433       m_indices = other.m_indices;
    434       return *this;
    435     }
    436     #endif
    437 
    438     /** const version of indices(). */
    439     const IndicesType& indices() const { return m_indices; }
    440     /** \returns a reference to the stored array representing the permutation. */
    441     IndicesType& indices() { return m_indices; }
    442 
    443   protected:
    444 
    445     IndicesType m_indices;
    446 };
    447 
    448 /** \class PermutationWrapper
    449   * \ingroup Core_Module
    450   *
    451   * \brief Class to view a vector of integers as a permutation matrix
    452   *
    453   * \param _IndicesType the type of the vector of integer (can be any compatible expression)
    454   *
    455   * This class allows to view any vector expression of integers as a permutation matrix.
    456   *
    457   * \sa class PermutationBase, class PermutationMatrix
    458   */
    459 
    460 struct PermutationStorage {};
    461 
    462 template<typename _IndicesType> class TranspositionsWrapper;
    463 namespace internal {
    464 template<typename _IndicesType>
    465 struct traits<PermutationWrapper<_IndicesType> >
    466 {
    467   typedef PermutationStorage StorageKind;
    468   typedef typename _IndicesType::Scalar Scalar;
    469   typedef typename _IndicesType::Scalar Index;
    470   typedef _IndicesType IndicesType;
    471   enum {
    472     RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
    473     ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
    474     MaxRowsAtCompileTime = IndicesType::MaxRowsAtCompileTime,
    475     MaxColsAtCompileTime = IndicesType::MaxColsAtCompileTime,
    476     Flags = 0,
    477     CoeffReadCost = _IndicesType::CoeffReadCost
    478   };
    479 };
    480 }
    481 
    482 template<typename _IndicesType>
    483 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
    484 {
    485     typedef PermutationBase<PermutationWrapper> Base;
    486     typedef internal::traits<PermutationWrapper> Traits;
    487   public:
    488 
    489     #ifndef EIGEN_PARSED_BY_DOXYGEN
    490     typedef typename Traits::IndicesType IndicesType;
    491     #endif
    492 
    493     inline PermutationWrapper(const IndicesType& a_indices)
    494       : m_indices(a_indices)
    495     {}
    496 
    497     /** const version of indices(). */
    498     const typename internal::remove_all<typename IndicesType::Nested>::type&
    499     indices() const { return m_indices; }
    500 
    501   protected:
    502 
    503     typename IndicesType::Nested m_indices;
    504 };
    505 
    506 /** \returns the matrix with the permutation applied to the columns.
    507   */
    508 template<typename Derived, typename PermutationDerived>
    509 inline const internal::permut_matrix_product_retval<PermutationDerived, Derived, OnTheRight>
    510 operator*(const MatrixBase<Derived>& matrix,
    511           const PermutationBase<PermutationDerived> &permutation)
    512 {
    513   return internal::permut_matrix_product_retval
    514            <PermutationDerived, Derived, OnTheRight>
    515            (permutation.derived(), matrix.derived());
    516 }
    517 
    518 /** \returns the matrix with the permutation applied to the rows.
    519   */
    520 template<typename Derived, typename PermutationDerived>
    521 inline const internal::permut_matrix_product_retval
    522                <PermutationDerived, Derived, OnTheLeft>
    523 operator*(const PermutationBase<PermutationDerived> &permutation,
    524           const MatrixBase<Derived>& matrix)
    525 {
    526   return internal::permut_matrix_product_retval
    527            <PermutationDerived, Derived, OnTheLeft>
    528            (permutation.derived(), matrix.derived());
    529 }
    530 
    531 namespace internal {
    532 
    533 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
    534 struct traits<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
    535 {
    536   typedef typename MatrixType::PlainObject ReturnType;
    537 };
    538 
    539 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
    540 struct permut_matrix_product_retval
    541  : public ReturnByValue<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
    542 {
    543     typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
    544     typedef typename MatrixType::Index Index;
    545 
    546     permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
    547       : m_permutation(perm), m_matrix(matrix)
    548     {}
    549 
    550     inline Index rows() const { return m_matrix.rows(); }
    551     inline Index cols() const { return m_matrix.cols(); }
    552 
    553     template<typename Dest> inline void evalTo(Dest& dst) const
    554     {
    555       const Index n = Side==OnTheLeft ? rows() : cols();
    556       // FIXME we need an is_same for expression that is not sensitive to constness. For instance
    557       // is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true.
    558       if(    is_same<MatrixTypeNestedCleaned,Dest>::value
    559           && blas_traits<MatrixTypeNestedCleaned>::HasUsableDirectAccess
    560           && blas_traits<Dest>::HasUsableDirectAccess
    561           && extract_data(dst) == extract_data(m_matrix))
    562       {
    563         // apply the permutation inplace
    564         Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
    565         mask.fill(false);
    566         Index r = 0;
    567         while(r < m_permutation.size())
    568         {
    569           // search for the next seed
    570           while(r<m_permutation.size() && mask[r]) r++;
    571           if(r>=m_permutation.size())
    572             break;
    573           // we got one, let's follow it until we are back to the seed
    574           Index k0 = r++;
    575           Index kPrev = k0;
    576           mask.coeffRef(k0) = true;
    577           for(Index k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
    578           {
    579                   Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
    580             .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
    581                        (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
    582 
    583             mask.coeffRef(k) = true;
    584             kPrev = k;
    585           }
    586         }
    587       }
    588       else
    589       {
    590         for(int i = 0; i < n; ++i)
    591         {
    592           Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
    593                (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
    594 
    595           =
    596 
    597           Block<const MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
    598                (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
    599         }
    600       }
    601     }
    602 
    603   protected:
    604     const PermutationType& m_permutation;
    605     typename MatrixType::Nested m_matrix;
    606 };
    607 
    608 /* Template partial specialization for transposed/inverse permutations */
    609 
    610 template<typename Derived>
    611 struct traits<Transpose<PermutationBase<Derived> > >
    612  : traits<Derived>
    613 {};
    614 
    615 } // end namespace internal
    616 
    617 template<typename Derived>
    618 class Transpose<PermutationBase<Derived> >
    619   : public EigenBase<Transpose<PermutationBase<Derived> > >
    620 {
    621     typedef Derived PermutationType;
    622     typedef typename PermutationType::IndicesType IndicesType;
    623     typedef typename PermutationType::PlainPermutationType PlainPermutationType;
    624   public:
    625 
    626     #ifndef EIGEN_PARSED_BY_DOXYGEN
    627     typedef internal::traits<PermutationType> Traits;
    628     typedef typename Derived::DenseMatrixType DenseMatrixType;
    629     enum {
    630       Flags = Traits::Flags,
    631       CoeffReadCost = Traits::CoeffReadCost,
    632       RowsAtCompileTime = Traits::RowsAtCompileTime,
    633       ColsAtCompileTime = Traits::ColsAtCompileTime,
    634       MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
    635       MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
    636     };
    637     typedef typename Traits::Scalar Scalar;
    638     #endif
    639 
    640     Transpose(const PermutationType& p) : m_permutation(p) {}
    641 
    642     inline int rows() const { return m_permutation.rows(); }
    643     inline int cols() const { return m_permutation.cols(); }
    644 
    645     #ifndef EIGEN_PARSED_BY_DOXYGEN
    646     template<typename DenseDerived>
    647     void evalTo(MatrixBase<DenseDerived>& other) const
    648     {
    649       other.setZero();
    650       for (int i=0; i<rows();++i)
    651         other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1);
    652     }
    653     #endif
    654 
    655     /** \return the equivalent permutation matrix */
    656     PlainPermutationType eval() const { return *this; }
    657 
    658     DenseMatrixType toDenseMatrix() const { return *this; }
    659 
    660     /** \returns the matrix with the inverse permutation applied to the columns.
    661       */
    662     template<typename OtherDerived> friend
    663     inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>
    664     operator*(const MatrixBase<OtherDerived>& matrix, const Transpose& trPerm)
    665     {
    666       return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>(trPerm.m_permutation, matrix.derived());
    667     }
    668 
    669     /** \returns the matrix with the inverse permutation applied to the rows.
    670       */
    671     template<typename OtherDerived>
    672     inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>
    673     operator*(const MatrixBase<OtherDerived>& matrix) const
    674     {
    675       return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>(m_permutation, matrix.derived());
    676     }
    677 
    678     const PermutationType& nestedPermutation() const { return m_permutation; }
    679 
    680   protected:
    681     const PermutationType& m_permutation;
    682 };
    683 
    684 template<typename Derived>
    685 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
    686 {
    687   return derived();
    688 }
    689 
    690 } // end namespace Eigen
    691 
    692 #endif // EIGEN_PERMUTATIONMATRIX_H
    693