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@gmail.com>
5// Copyright (C) 2009-2011 Gael Guennebaud <gael.guennebaud@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
14namespace Eigen {
15
16template<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
42namespace internal {
43
44template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
45struct permut_matrix_product_retval;
46template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
47struct permut_sparsematrix_product_retval;
48enum PermPermProduct_t {PermPermProduct};
49
50} // end namespace internal
51
52template<typename Derived>
53class 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 indices().size(); }
109
110    /** \returns the number of columns */
111    inline Index cols() const { return 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 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 size)
143    {
144      indices().resize(size);
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 size)
157    {
158      resize(size);
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
272namespace internal {
273template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
274struct 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
282template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
283class 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>& indices) : m_indices(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
386namespace internal {
387template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
388struct 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
396template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
397class 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* indices)
410      : m_indices(indices)
411    {}
412
413    inline Map(const Index* indices, Index size)
414      : m_indices(indices,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
460struct PermutationStorage {};
461
462template<typename _IndicesType> class TranspositionsWrapper;
463namespace internal {
464template<typename _IndicesType>
465struct 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
482template<typename _IndicesType>
483class 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& indices)
494      : m_indices(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  */
508template<typename Derived, typename PermutationDerived>
509inline const internal::permut_matrix_product_retval<PermutationDerived, Derived, OnTheRight>
510operator*(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  */
520template<typename Derived, typename PermutationDerived>
521inline const internal::permut_matrix_product_retval
522               <PermutationDerived, Derived, OnTheLeft>
523operator*(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
531namespace internal {
532
533template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
534struct traits<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
535{
536  typedef typename MatrixType::PlainObject ReturnType;
537};
538
539template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
540struct 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
545    permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
546      : m_permutation(perm), m_matrix(matrix)
547    {}
548
549    inline int rows() const { return m_matrix.rows(); }
550    inline int cols() const { return m_matrix.cols(); }
551
552    template<typename Dest> inline void evalTo(Dest& dst) const
553    {
554      const int n = Side==OnTheLeft ? rows() : cols();
555
556      if(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix))
557      {
558        // apply the permutation inplace
559        Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
560        mask.fill(false);
561        int r = 0;
562        while(r < m_permutation.size())
563        {
564          // search for the next seed
565          while(r<m_permutation.size() && mask[r]) r++;
566          if(r>=m_permutation.size())
567            break;
568          // we got one, let's follow it until we are back to the seed
569          int k0 = r++;
570          int kPrev = k0;
571          mask.coeffRef(k0) = true;
572          for(int k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
573          {
574                  Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
575            .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
576                       (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
577
578            mask.coeffRef(k) = true;
579            kPrev = k;
580          }
581        }
582      }
583      else
584      {
585        for(int i = 0; i < n; ++i)
586        {
587          Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
588               (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
589
590          =
591
592          Block<const MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
593               (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
594        }
595      }
596    }
597
598  protected:
599    const PermutationType& m_permutation;
600    typename MatrixType::Nested m_matrix;
601};
602
603/* Template partial specialization for transposed/inverse permutations */
604
605template<typename Derived>
606struct traits<Transpose<PermutationBase<Derived> > >
607 : traits<Derived>
608{};
609
610} // end namespace internal
611
612template<typename Derived>
613class Transpose<PermutationBase<Derived> >
614  : public EigenBase<Transpose<PermutationBase<Derived> > >
615{
616    typedef Derived PermutationType;
617    typedef typename PermutationType::IndicesType IndicesType;
618    typedef typename PermutationType::PlainPermutationType PlainPermutationType;
619  public:
620
621    #ifndef EIGEN_PARSED_BY_DOXYGEN
622    typedef internal::traits<PermutationType> Traits;
623    typedef typename Derived::DenseMatrixType DenseMatrixType;
624    enum {
625      Flags = Traits::Flags,
626      CoeffReadCost = Traits::CoeffReadCost,
627      RowsAtCompileTime = Traits::RowsAtCompileTime,
628      ColsAtCompileTime = Traits::ColsAtCompileTime,
629      MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
630      MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
631    };
632    typedef typename Traits::Scalar Scalar;
633    #endif
634
635    Transpose(const PermutationType& p) : m_permutation(p) {}
636
637    inline int rows() const { return m_permutation.rows(); }
638    inline int cols() const { return m_permutation.cols(); }
639
640    #ifndef EIGEN_PARSED_BY_DOXYGEN
641    template<typename DenseDerived>
642    void evalTo(MatrixBase<DenseDerived>& other) const
643    {
644      other.setZero();
645      for (int i=0; i<rows();++i)
646        other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1);
647    }
648    #endif
649
650    /** \return the equivalent permutation matrix */
651    PlainPermutationType eval() const { return *this; }
652
653    DenseMatrixType toDenseMatrix() const { return *this; }
654
655    /** \returns the matrix with the inverse permutation applied to the columns.
656      */
657    template<typename OtherDerived> friend
658    inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>
659    operator*(const MatrixBase<OtherDerived>& matrix, const Transpose& trPerm)
660    {
661      return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>(trPerm.m_permutation, matrix.derived());
662    }
663
664    /** \returns the matrix with the inverse permutation applied to the rows.
665      */
666    template<typename OtherDerived>
667    inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>
668    operator*(const MatrixBase<OtherDerived>& matrix) const
669    {
670      return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>(m_permutation, matrix.derived());
671    }
672
673    const PermutationType& nestedPermutation() const { return m_permutation; }
674
675  protected:
676    const PermutationType& m_permutation;
677};
678
679template<typename Derived>
680const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
681{
682  return derived();
683}
684
685} // end namespace Eigen
686
687#endif // EIGEN_PERMUTATIONMATRIX_H
688