PermutationMatrix.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 // Copyright (C) 2009-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
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_PERMUTATIONMATRIX_H
00012 #define EIGEN_PERMUTATIONMATRIX_H
00013 
00014 namespace Eigen { 
00015 
00016 template<int RowCol,typename IndicesType,typename MatrixType, typename StorageKind> class PermutedImpl;
00017 
00018 /** \class PermutationBase
00019   * \ingroup Core_Module
00020   *
00021   * \brief Base class for permutations
00022   *
00023   * \param Derived the derived class
00024   *
00025   * This class is the base class for all expressions representing a permutation matrix,
00026   * internally stored as a vector of integers.
00027   * The convention followed here is that if \f$ \sigma \f$ is a permutation, the corresponding permutation matrix
00028   * \f$ P_\sigma \f$ is such that if \f$ (e_1,\ldots,e_p) \f$ is the canonical basis, we have:
00029   *  \f[ P_\sigma(e_i) = e_{\sigma(i)}. \f]
00030   * This convention ensures that for any two permutations \f$ \sigma, \tau \f$, we have:
00031   *  \f[ P_{\sigma\circ\tau} = P_\sigma P_\tau. \f]
00032   *
00033   * Permutation matrices are square and invertible.
00034   *
00035   * Notice that in addition to the member functions and operators listed here, there also are non-member
00036   * operator* to multiply any kind of permutation object with any kind of matrix expression (MatrixBase)
00037   * on either side.
00038   *
00039   * \sa class PermutationMatrix, class PermutationWrapper
00040   */
00041 
00042 namespace internal {
00043 
00044 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
00045 struct permut_matrix_product_retval;
00046 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
00047 struct permut_sparsematrix_product_retval;
00048 enum PermPermProduct_t {PermPermProduct};
00049 
00050 } // end namespace internal
00051 
00052 template<typename Derived>
00053 class PermutationBase : public EigenBase<Derived>
00054 {
00055     typedef internal::traits<Derived> Traits;
00056     typedef EigenBase<Derived> Base;
00057   public:
00058 
00059     #ifndef EIGEN_PARSED_BY_DOXYGEN
00060     typedef typename Traits::IndicesType IndicesType;
00061     enum {
00062       Flags = Traits::Flags,
00063       CoeffReadCost = Traits::CoeffReadCost,
00064       RowsAtCompileTime = Traits::RowsAtCompileTime,
00065       ColsAtCompileTime = Traits::ColsAtCompileTime,
00066       MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
00067       MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
00068     };
00069     typedef typename Traits::Scalar Scalar;
00070     typedef typename Traits::Index Index;
00071     typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime,0,MaxRowsAtCompileTime,MaxColsAtCompileTime>
00072             DenseMatrixType;
00073     typedef PermutationMatrix<IndicesType::SizeAtCompileTime,IndicesType::MaxSizeAtCompileTime,Index>
00074             PlainPermutationType;
00075     using Base::derived ;
00076     #endif
00077 
00078     /** Copies the other permutation into *this */
00079     template<typename OtherDerived>
00080     Derived& operator=(const PermutationBase<OtherDerived>& other)
00081     {
00082       indices() = other.indices();
00083       return derived ();
00084     }
00085 
00086     /** Assignment from the Transpositions \a tr */
00087     template<typename OtherDerived>
00088     Derived& operator=(const TranspositionsBase<OtherDerived>& tr)
00089     {
00090       setIdentity(tr.size());
00091       for(Index k=size ()-1; k>=0; --k)
00092         applyTranspositionOnTheRight(k,tr.coeff(k));
00093       return derived ();
00094     }
00095 
00096     #ifndef EIGEN_PARSED_BY_DOXYGEN
00097     /** This is a special case of the templated operator=. Its purpose is to
00098       * prevent a default operator= from hiding the templated operator=.
00099       */
00100     Derived& operator=(const PermutationBase& other)
00101     {
00102       indices() = other.indices();
00103       return derived ();
00104     }
00105     #endif
00106 
00107     /** \returns the number of rows */
00108     inline Index rows () const { return Index(indices().size ()); }
00109 
00110     /** \returns the number of columns */
00111     inline Index cols () const { return Index(indices().size ()); }
00112 
00113     /** \returns the size of a side of the respective square matrix, i.e., the number of indices */
00114     inline Index size () const { return Index(indices().size ()); }
00115 
00116     #ifndef EIGEN_PARSED_BY_DOXYGEN
00117     template<typename DenseDerived>
00118     void evalTo(MatrixBase<DenseDerived>& other) const
00119     {
00120       other.setZero();
00121       for (int i=0; i<rows ();++i)
00122         other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
00123     }
00124     #endif
00125 
00126     /** \returns a Matrix object initialized from this permutation matrix. Notice that it
00127       * is inefficient to return this Matrix object by value. For efficiency, favor using
00128       * the Matrix constructor taking EigenBase objects.
00129       */
00130     DenseMatrixType toDenseMatrix () const
00131     {
00132       return derived ();
00133     }
00134 
00135     /** const version of indices(). */
00136     const IndicesType& indices() const { return derived ().indices(); }
00137     /** \returns a reference to the stored array representing the permutation. */
00138     IndicesType& indices () { return derived ().indices(); }
00139 
00140     /** Resizes to given size.
00141       */
00142     inline void resize(Index newSize)
00143     {
00144       indices().resize(newSize);
00145     }
00146 
00147     /** Sets *this to be the identity permutation matrix */
00148     void setIdentity()
00149     {
00150       for(Index i = 0; i < size (); ++i)
00151         indices().coeffRef(i) = i;
00152     }
00153 
00154     /** Sets *this to be the identity permutation matrix of given size.
00155       */
00156     void setIdentity(Index newSize)
00157     {
00158       resize(newSize);
00159       setIdentity();
00160     }
00161 
00162     /** Multiplies *this by the transposition \f$(ij)\f$ on the left.
00163       *
00164       * \returns a reference to *this.
00165       *
00166       * \warning This is much slower than applyTranspositionOnTheRight(int,int):
00167       * this has linear complexity and requires a lot of branching.
00168       *
00169       * \sa applyTranspositionOnTheRight(int,int)
00170       */
00171     Derived& applyTranspositionOnTheLeft(Index i, Index j)
00172     {
00173       eigen_assert(i>=0 && j>=0 && i<size () && j<size ());
00174       for(Index k = 0; k < size (); ++k)
00175       {
00176         if(indices().coeff(k) == i) indices().coeffRef(k) = j;
00177         else if(indices().coeff(k) == j) indices().coeffRef(k) = i;
00178       }
00179       return derived ();
00180     }
00181 
00182     /** Multiplies *this by the transposition \f$(ij)\f$ on the right.
00183       *
00184       * \returns a reference to *this.
00185       *
00186       * This is a fast operation, it only consists in swapping two indices.
00187       *
00188       * \sa applyTranspositionOnTheLeft(int,int)
00189       */
00190     Derived& applyTranspositionOnTheRight(Index i, Index j)
00191     {
00192       eigen_assert(i>=0 && j>=0 && i<size () && j<size ());
00193       std::swap(indices().coeffRef(i), indices().coeffRef(j));
00194       return derived ();
00195     }
00196 
00197     /** \returns the inverse permutation matrix.
00198       *
00199       * \note \note_try_to_help_rvo
00200       */
00201     inline Transpose<PermutationBase> inverse () const
00202     { return derived (); }
00203     /** \returns the tranpose permutation matrix.
00204       *
00205       * \note \note_try_to_help_rvo
00206       */
00207     inline Transpose<PermutationBase> transpose () const
00208     { return derived (); }
00209 
00210     /**** multiplication helpers to hopefully get RVO ****/
00211 
00212   
00213 #ifndef EIGEN_PARSED_BY_DOXYGEN
00214   protected:
00215     template<typename OtherDerived>
00216     void assignTranspose(const PermutationBase<OtherDerived>& other)
00217     {
00218       for (int i=0; i<rows ();++i) indices().coeffRef(other.indices().coeff(i)) = i;
00219     }
00220     template<typename Lhs,typename Rhs>
00221     void assignProduct(const Lhs& lhs, const Rhs& rhs)
00222     {
00223       eigen_assert(lhs.cols() == rhs.rows());
00224       for (int i=0; i<rows ();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
00225     }
00226 #endif
00227 
00228   public:
00229 
00230     /** \returns the product permutation matrix.
00231       *
00232       * \note \note_try_to_help_rvo
00233       */
00234     template<typename Other>
00235     inline PlainPermutationType operator* (const PermutationBase<Other>& other) const
00236     { return PlainPermutationType(internal::PermPermProduct, derived (), other.derived ()); }
00237 
00238     /** \returns the product of a permutation with another inverse permutation.
00239       *
00240       * \note \note_try_to_help_rvo
00241       */
00242     template<typename Other>
00243     inline PlainPermutationType operator* (const Transpose<PermutationBase<Other> >& other) const
00244     { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
00245 
00246     /** \returns the product of an inverse permutation with another permutation.
00247       *
00248       * \note \note_try_to_help_rvo
00249       */
00250     template<typename Other> friend
00251     inline PlainPermutationType operator* (const Transpose<PermutationBase<Other> >& other, const PermutationBase& perm)
00252     { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
00253     
00254     /** \returns the determinant of the permutation matrix, which is either 1 or -1 depending on the parity of the permutation.
00255       *
00256       * This function is O(\c n) procedure allocating a buffer of \c n booleans.
00257       */
00258     Index determinant () const
00259     {
00260       Index res = 1;
00261       Index n = size ();
00262       Matrix<bool,RowsAtCompileTime,1,0,MaxRowsAtCompileTime> mask(n);
00263       mask.fill(false);
00264       Index r = 0;
00265       while(r < n)
00266       {
00267         // search for the next seed
00268         while(r<n && mask[r]) r++;
00269         if(r>=n)
00270           break;
00271         // we got one, let's follow it until we are back to the seed
00272         Index k0 = r++;
00273         mask.coeffRef(k0) = true;
00274         for(Index k=indices().coeff(k0); k!=k0; k=indices().coeff(k))
00275         {
00276           mask.coeffRef(k) = true;
00277           res = -res;
00278         }
00279       }
00280       return res;
00281     }
00282 
00283   protected:
00284 
00285 };
00286 
00287 /** \class PermutationMatrix
00288   * \ingroup Core_Module
00289   *
00290   * \brief Permutation matrix
00291   *
00292   * \param SizeAtCompileTime the number of rows/cols, or Dynamic
00293   * \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.
00294   * \param IndexType the interger type of the indices
00295   *
00296   * This class represents a permutation matrix, internally stored as a vector of integers.
00297   *
00298   * \sa class PermutationBase, class PermutationWrapper, class DiagonalMatrix
00299   */
00300 
00301 namespace internal {
00302 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
00303 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
00304  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
00305 {
00306   typedef IndexType Index;
00307   typedef Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
00308 };
00309 }
00310 
00311 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
00312 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
00313 {
00314     typedef PermutationBase<PermutationMatrix> Base;
00315     typedef internal::traits<PermutationMatrix> Traits;
00316   public:
00317 
00318     #ifndef EIGEN_PARSED_BY_DOXYGEN
00319     typedef typename Traits::IndicesType IndicesType;
00320     #endif
00321 
00322     inline PermutationMatrix()
00323     {}
00324 
00325     /** Constructs an uninitialized permutation matrix of given size.
00326       */
00327     inline PermutationMatrix(int size ) : m_indices(size)
00328     {}
00329 
00330     /** Copy constructor. */
00331     template<typename OtherDerived>
00332     inline PermutationMatrix(const PermutationBase<OtherDerived>& other)
00333       : m_indices(other.indices()) {}
00334 
00335     #ifndef EIGEN_PARSED_BY_DOXYGEN
00336     /** Standard copy constructor. Defined only to prevent a default copy constructor
00337       * from hiding the other templated constructor */
00338     inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
00339     #endif
00340 
00341     /** Generic constructor from expression of the indices. The indices
00342       * array has the meaning that the permutations sends each integer i to indices[i].
00343       *
00344       * \warning It is your responsibility to check that the indices array that you passes actually
00345       * describes a permutation, i.e., each value between 0 and n-1 occurs exactly once, where n is the
00346       * array's size.
00347       */
00348     template<typename Other>
00349     explicit inline PermutationMatrix(const MatrixBase<Other>& a_indices) : m_indices(a_indices)
00350     {}
00351 
00352     /** Convert the Transpositions \a tr to a permutation matrix */
00353     template<typename Other>
00354     explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
00355       : m_indices(tr.size ())
00356     {
00357       *this = tr;
00358     }
00359 
00360     /** Copies the other permutation into *this */
00361     template<typename Other>
00362     PermutationMatrix& operator=(const PermutationBase<Other>& other)
00363     {
00364       m_indices = other.indices();
00365       return *this;
00366     }
00367 
00368     /** Assignment from the Transpositions \a tr */
00369     template<typename Other>
00370     PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
00371     {
00372       return Base::operator=(tr.derived());
00373     }
00374 
00375     #ifndef EIGEN_PARSED_BY_DOXYGEN
00376     /** This is a special case of the templated operator=. Its purpose is to
00377       * prevent a default operator= from hiding the templated operator=.
00378       */
00379     PermutationMatrix& operator=(const PermutationMatrix& other)
00380     {
00381       m_indices = other.m_indices;
00382       return *this;
00383     }
00384     #endif
00385 
00386     /** const version of indices(). */
00387     const IndicesType& indices() const { return m_indices; }
00388     /** \returns a reference to the stored array representing the permutation. */
00389     IndicesType& indices () { return m_indices; }
00390 
00391 
00392     /**** multiplication helpers to hopefully get RVO ****/
00393 
00394 #ifndef EIGEN_PARSED_BY_DOXYGEN
00395     template<typename Other>
00396     PermutationMatrix(const Transpose<PermutationBase<Other> >& other)
00397       : m_indices(other.nestedPermutation().size ())
00398     {
00399       for (int i=0; i<m_indices.size();++i) m_indices.coeffRef(other.nestedPermutation().indices().coeff(i)) = i;
00400     }
00401     template<typename Lhs,typename Rhs>
00402     PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
00403       : m_indices(lhs.indices().size ())
00404     {
00405       Base::assignProduct(lhs,rhs);
00406     }
00407 #endif
00408 
00409   protected:
00410 
00411     IndicesType m_indices;
00412 };
00413 
00414 
00415 namespace internal {
00416 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
00417 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
00418  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
00419 {
00420   typedef IndexType Index;
00421   typedef Map<const Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType;
00422 };
00423 }
00424 
00425 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
00426 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess>
00427   : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
00428 {
00429     typedef PermutationBase<Map> Base;
00430     typedef internal::traits<Map> Traits;
00431   public:
00432 
00433     #ifndef EIGEN_PARSED_BY_DOXYGEN
00434     typedef typename Traits::IndicesType IndicesType;
00435     typedef typename IndicesType::Scalar Index;
00436     #endif
00437 
00438     inline Map(const Index* indicesPtr)
00439       : m_indices(indicesPtr)
00440     {}
00441 
00442     inline Map(const Index* indicesPtr, Index size)
00443       : m_indices(indicesPtr,size)
00444     {}
00445 
00446     /** Copies the other permutation into *this */
00447     template<typename Other>
00448     Map& operator=(const PermutationBase<Other>& other)
00449     { return Base::operator=(other.derived()); }
00450 
00451     /** Assignment from the Transpositions \a tr */
00452     template<typename Other>
00453     Map& operator=(const TranspositionsBase<Other>& tr)
00454     { return Base::operator=(tr.derived()); }
00455 
00456     #ifndef EIGEN_PARSED_BY_DOXYGEN
00457     /** This is a special case of the templated operator=. Its purpose is to
00458       * prevent a default operator= from hiding the templated operator=.
00459       */
00460     Map& operator=(const Map& other)
00461     {
00462       m_indices = other.m_indices;
00463       return *this;
00464     }
00465     #endif
00466 
00467     /** const version of indices(). */
00468     const IndicesType& indices() const { return m_indices; }
00469     /** \returns a reference to the stored array representing the permutation. */
00470     IndicesType& indices() { return m_indices; }
00471 
00472   protected:
00473 
00474     IndicesType m_indices;
00475 };
00476 
00477 /** \class PermutationWrapper
00478   * \ingroup Core_Module
00479   *
00480   * \brief Class to view a vector of integers as a permutation matrix
00481   *
00482   * \param _IndicesType the type of the vector of integer (can be any compatible expression)
00483   *
00484   * This class allows to view any vector expression of integers as a permutation matrix.
00485   *
00486   * \sa class PermutationBase, class PermutationMatrix
00487   */
00488 
00489 struct PermutationStorage {};
00490 
00491 template<typename _IndicesType> class TranspositionsWrapper;
00492 namespace internal {
00493 template<typename _IndicesType>
00494 struct traits<PermutationWrapper<_IndicesType> >
00495 {
00496   typedef PermutationStorage StorageKind;
00497   typedef typename _IndicesType::Scalar Scalar;
00498   typedef typename _IndicesType::Scalar Index;
00499   typedef _IndicesType IndicesType;
00500   enum {
00501     RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
00502     ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
00503     MaxRowsAtCompileTime = IndicesType::MaxRowsAtCompileTime,
00504     MaxColsAtCompileTime = IndicesType::MaxColsAtCompileTime,
00505     Flags = 0,
00506     CoeffReadCost = _IndicesType::CoeffReadCost
00507   };
00508 };
00509 }
00510 
00511 template<typename _IndicesType>
00512 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
00513 {
00514     typedef PermutationBase<PermutationWrapper> Base;
00515     typedef internal::traits<PermutationWrapper> Traits;
00516   public:
00517 
00518     #ifndef EIGEN_PARSED_BY_DOXYGEN
00519     typedef typename Traits::IndicesType IndicesType;
00520     #endif
00521 
00522     inline PermutationWrapper(const IndicesType& a_indices)
00523       : m_indices(a_indices)
00524     {}
00525 
00526     /** const version of indices(). */
00527     const typename internal::remove_all<typename IndicesType::Nested>::type&
00528     indices() const { return m_indices; }
00529 
00530   protected:
00531 
00532     typename IndicesType::Nested m_indices;
00533 };
00534 
00535 /** \returns the matrix with the permutation applied to the columns.
00536   */
00537 template<typename Derived, typename PermutationDerived>
00538 inline const internal::permut_matrix_product_retval<PermutationDerived, Derived, OnTheRight>
00539 operator* (const MatrixBase<Derived>& matrix,
00540           const PermutationBase<PermutationDerived> &permutation)
00541 {
00542   return internal::permut_matrix_product_retval
00543            <PermutationDerived, Derived, OnTheRight>
00544            (permutation.derived (), matrix.derived());
00545 }
00546 
00547 /** \returns the matrix with the permutation applied to the rows.
00548   */
00549 template<typename Derived, typename PermutationDerived>
00550 inline const internal::permut_matrix_product_retval
00551                <PermutationDerived, Derived, OnTheLeft>
00552 operator* (const PermutationBase<PermutationDerived> &permutation,
00553           const MatrixBase<Derived>& matrix)
00554 {
00555   return internal::permut_matrix_product_retval
00556            <PermutationDerived, Derived, OnTheLeft>
00557            (permutation.derived (), matrix.derived());
00558 }
00559 
00560 namespace internal {
00561 
00562 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
00563 struct traits<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
00564 {
00565   typedef typename MatrixType::PlainObject ReturnType;
00566 };
00567 
00568 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
00569 struct permut_matrix_product_retval
00570  : public ReturnByValue<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
00571 {
00572     typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
00573     typedef typename MatrixType::Index Index;
00574 
00575     permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
00576       : m_permutation(perm), m_matrix(matrix)
00577     {}
00578 
00579     inline Index rows() const { return m_matrix.rows(); }
00580     inline Index cols() const { return m_matrix.cols(); }
00581 
00582     template<typename Dest> inline void evalTo(Dest& dst) const
00583     {
00584       const Index n = Side==OnTheLeft ? rows() : cols();
00585       // FIXME we need an is_same for expression that is not sensitive to constness. For instance
00586       // is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true.
00587       if(    is_same<MatrixTypeNestedCleaned,Dest>::value
00588           && blas_traits<MatrixTypeNestedCleaned>::HasUsableDirectAccess
00589           && blas_traits<Dest>::HasUsableDirectAccess
00590           && extract_data(dst) == extract_data(m_matrix))
00591       {
00592         // apply the permutation inplace
00593         Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
00594         mask.fill(false);
00595         Index r = 0;
00596         while(r < m_permutation.size())
00597         {
00598           // search for the next seed
00599           while(r<m_permutation.size() && mask[r]) r++;
00600           if(r>=m_permutation.size())
00601             break;
00602           // we got one, let's follow it until we are back to the seed
00603           Index k0 = r++;
00604           Index kPrev = k0;
00605           mask.coeffRef(k0) = true;
00606           for(Index k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
00607           {
00608                   Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
00609             .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
00610                        (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
00611 
00612             mask.coeffRef(k) = true;
00613             kPrev = k;
00614           }
00615         }
00616       }
00617       else
00618       {
00619         for(int i = 0; i < n; ++i)
00620         {
00621           Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
00622                (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
00623 
00624           =
00625 
00626           Block<const MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
00627                (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
00628         }
00629       }
00630     }
00631 
00632   protected:
00633     const PermutationType& m_permutation;
00634     typename MatrixType::Nested m_matrix;
00635 };
00636 
00637 /* Template partial specialization for transposed/inverse permutations */
00638 
00639 template<typename Derived>
00640 struct traits<Transpose<PermutationBase<Derived> > >
00641  : traits<Derived>
00642 {};
00643 
00644 } // end namespace internal
00645 
00646 template<typename Derived>
00647 class Transpose<PermutationBase<Derived> >
00648   : public EigenBase<Transpose<PermutationBase<Derived> > >
00649 {
00650     typedef Derived PermutationType;
00651     typedef typename PermutationType::IndicesType IndicesType;
00652     typedef typename PermutationType::PlainPermutationType PlainPermutationType;
00653   public:
00654 
00655     #ifndef EIGEN_PARSED_BY_DOXYGEN
00656     typedef internal::traits<PermutationType> Traits;
00657     typedef typename Derived::DenseMatrixType DenseMatrixType;
00658     enum {
00659       Flags = Traits::Flags,
00660       CoeffReadCost = Traits::CoeffReadCost,
00661       RowsAtCompileTime = Traits::RowsAtCompileTime,
00662       ColsAtCompileTime = Traits::ColsAtCompileTime,
00663       MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
00664       MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
00665     };
00666     typedef typename Traits::Scalar Scalar;
00667     #endif
00668 
00669     Transpose(const PermutationType& p) : m_permutation(p) {}
00670 
00671     inline int rows() const { return m_permutation.rows(); }
00672     inline int cols() const { return m_permutation.cols(); }
00673 
00674     #ifndef EIGEN_PARSED_BY_DOXYGEN
00675     template<typename DenseDerived>
00676     void evalTo(MatrixBase<DenseDerived>& other) const
00677     {
00678       other.setZero();
00679       for (int i=0; i<rows();++i)
00680         other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1);
00681     }
00682     #endif
00683 
00684     /** \return the equivalent permutation matrix */
00685     PlainPermutationType eval() const { return *this; }
00686 
00687     DenseMatrixType toDenseMatrix() const { return *this; }
00688 
00689     /** \returns the matrix with the inverse permutation applied to the columns.
00690       */
00691     template<typename OtherDerived> friend
00692     inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>
00693     operator* (const MatrixBase<OtherDerived>& matrix, const Transpose& trPerm)
00694     {
00695       return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>(trPerm.m_permutation, matrix.derived());
00696     }
00697 
00698     /** \returns the matrix with the inverse permutation applied to the rows.
00699       */
00700     template<typename OtherDerived>
00701     inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>
00702     operator* (const MatrixBase<OtherDerived>& matrix) const
00703     {
00704       return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>(m_permutation, matrix.derived());
00705     }
00706 
00707     const PermutationType& nestedPermutation() const { return m_permutation; }
00708 
00709   protected:
00710     const PermutationType& m_permutation;
00711 };
00712 
00713 template<typename Derived>
00714 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
00715 {
00716   return derived();
00717 }
00718 
00719 } // end namespace Eigen
00720 
00721 #endif // EIGEN_PERMUTATIONMATRIX_H