src/Core/PermutationMatrix.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
0005 // Copyright (C) 2009-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_PERMUTATIONMATRIX_H
0012 #define EIGEN_PERMUTATIONMATRIX_H
0013
0014 namespace Eigen {
0015
0016 namespace internal {
0017
0018 enum PermPermProduct_t {PermPermProduct};
0019
0020 } // end namespace internal
0021
0022 /** \class PermutationBase
0023   * \ingroup Core_Module
0024   *
0025   * \brief Base class for permutations
0026   *
0027   * \tparam Derived the derived class
0028   *
0029   * This class is the base class for all expressions representing a permutation matrix,
0030   * internally stored as a vector of integers.
0031   * The convention followed here is that if \f$ \sigma \f$ is a permutation, the corresponding permutation matrix
0032   * \f$ P_\sigma \f$ is such that if \f$ (e_1,\ldots,e_p) \f$ is the canonical basis, we have:
0033   *  \f[ P_\sigma(e_i) = e_{\sigma(i)}. \f]
0034   * This convention ensures that for any two permutations \f$ \sigma, \tau \f$, we have:
0035   *  \f[ P_{\sigma\circ\tau} = P_\sigma P_\tau. \f]
0036   *
0037   * Permutation matrices are square and invertible.
0038   *
0039   * Notice that in addition to the member functions and operators listed here, there also are non-member
0040   * operator* to multiply any kind of permutation object with any kind of matrix expression (MatrixBase)
0041   * on either side.
0042   *
0043   * \sa class PermutationMatrix, class PermutationWrapper
0044   */
0045 template<typename Derived>
0046 class PermutationBase : public EigenBase<Derived>
0047 {
0048     typedef internal::traits<Derived> Traits;
0049     typedef EigenBase<Derived> Base;
0050   public:
0051
0052     #ifndef EIGEN_PARSED_BY_DOXYGEN
0053     typedef typename Traits::IndicesType IndicesType;
0054     enum {
0055       Flags = Traits::Flags,
0056       RowsAtCompileTime = Traits::RowsAtCompileTime,
0057       ColsAtCompileTime = Traits::ColsAtCompileTime,
0058       MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
0059       MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
0060     };
0061     typedef typename Traits::StorageIndex StorageIndex;
0062     typedef Matrix<StorageIndex,RowsAtCompileTime,ColsAtCompileTime,0,MaxRowsAtCompileTime,MaxColsAtCompileTime>
0063             DenseMatrixType;
0064     typedef PermutationMatrix<IndicesType::SizeAtCompileTime,IndicesType::MaxSizeAtCompileTime,StorageIndex>
0065             PlainPermutationType;
0066     typedef PlainPermutationType PlainObject;
0067     using Base::derived;
0068     typedef Inverse<Derived> InverseReturnType;
0069     typedef void Scalar;
0070     #endif
0071
0072     /** Copies the other permutation into *this */
0073     template<typename OtherDerived>
0074     Derived& operator=(const PermutationBase<OtherDerived>& other)
0075     {
0076       indices() = other.indices();
0077       return derived();
0078     }
0079
0080     /** Assignment from the Transpositions \a tr */
0081     template<typename OtherDerived>
0082     Derived& operator=(const TranspositionsBase<OtherDerived>& tr)
0083     {
0084       setIdentity(tr.size());
0085       for(Index k=size()-1; k>=0; --k)
0086         applyTranspositionOnTheRight(k,tr.coeff(k));
0087       return derived();
0088     }
0089
0090     /** \returns the number of rows */
0091     inline EIGEN_DEVICE_FUNC Index rows() const { return Index(indices().size()); }
0092
0093     /** \returns the number of columns */
0094     inline EIGEN_DEVICE_FUNC Index cols() const { return Index(indices().size()); }
0095
0096     /** \returns the size of a side of the respective square matrix, i.e., the number of indices */
0097     inline EIGEN_DEVICE_FUNC Index size() const { return Index(indices().size()); }
0098
0099     #ifndef EIGEN_PARSED_BY_DOXYGEN
0100     template<typename DenseDerived>
0101     void evalTo(MatrixBase<DenseDerived>& other) const
0102     {
0103       other.setZero();
0104       for (Index i=0; i<rows(); ++i)
0105         other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
0106     }
0107     #endif
0108
0109     /** \returns a Matrix object initialized from this permutation matrix. Notice that it
0110       * is inefficient to return this Matrix object by value. For efficiency, favor using
0111       * the Matrix constructor taking EigenBase objects.
0112       */
0113     DenseMatrixType toDenseMatrix() const
0114     {
0115       return derived();
0116     }
0117
0118     /** const version of indices(). */
0119     const IndicesType& indices() const { return derived().indices(); }
0120     /** \returns a reference to the stored array representing the permutation. */
0121     IndicesType& indices() { return derived().indices(); }
0122
0123     /** Resizes to given size.
0124       */
0125     inline void resize(Index newSize)
0126     {
0127       indices().resize(newSize);
0128     }
0129
0130     /** Sets *this to be the identity permutation matrix */
0131     void setIdentity()
0132     {
0133       StorageIndex n = StorageIndex(size());
0134       for(StorageIndex i = 0; i < n; ++i)
0135         indices().coeffRef(i) = i;
0136     }
0137
0138     /** Sets *this to be the identity permutation matrix of given size.
0139       */
0140     void setIdentity(Index newSize)
0141     {
0142       resize(newSize);
0143       setIdentity();
0144     }
0145
0146     /** Multiplies *this by the transposition \f$(ij)\f$ on the left.
0147       *
0148       * \returns a reference to *this.
0149       *
0150       * \warning This is much slower than applyTranspositionOnTheRight(Index,Index):
0151       * this has linear complexity and requires a lot of branching.
0152       *
0153       * \sa applyTranspositionOnTheRight(Index,Index)
0154       */
0155     Derived& applyTranspositionOnTheLeft(Index i, Index j)
0156     {
0157       eigen_assert(i>=0 && j>=0 && i<size() && j<size());
0158       for(Index k = 0; k < size(); ++k)
0159       {
0160         if(indices().coeff(k) == i) indices().coeffRef(k) = StorageIndex(j);
0161         else if(indices().coeff(k) == j) indices().coeffRef(k) = StorageIndex(i);
0162       }
0163       return derived();
0164     }
0165
0166     /** Multiplies *this by the transposition \f$(ij)\f$ on the right.
0167       *
0168       * \returns a reference to *this.
0169       *
0170       * This is a fast operation, it only consists in swapping two indices.
0171       *
0172       * \sa applyTranspositionOnTheLeft(Index,Index)
0173       */
0174     Derived& applyTranspositionOnTheRight(Index i, Index j)
0175     {
0176       eigen_assert(i>=0 && j>=0 && i<size() && j<size());
0177       std::swap(indices().coeffRef(i), indices().coeffRef(j));
0178       return derived();
0179     }
0180
0181     /** \returns the inverse permutation matrix.
0182       *
0183       * \note \blank \note_try_to_help_rvo
0184       */
0185     inline InverseReturnType inverse() const
0186     { return InverseReturnType(derived()); }
0187     /** \returns the tranpose permutation matrix.
0188       *
0189       * \note \blank \note_try_to_help_rvo
0190       */
0191     inline InverseReturnType transpose() const
0192     { return InverseReturnType(derived()); }
0193
0194     /**** multiplication helpers to hopefully get RVO ****/
0195
0196
0197 #ifndef EIGEN_PARSED_BY_DOXYGEN
0198   protected:
0199     template<typename OtherDerived>
0200     void assignTranspose(const PermutationBase<OtherDerived>& other)
0201     {
0202       for (Index i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i;
0203     }
0204     template<typename Lhs,typename Rhs>
0205     void assignProduct(const Lhs& lhs, const Rhs& rhs)
0206     {
0207       eigen_assert(lhs.cols() == rhs.rows());
0208       for (Index i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
0209     }
0210 #endif
0211
0212   public:
0213
0214     /** \returns the product permutation matrix.
0215       *
0216       * \note \blank \note_try_to_help_rvo
0217       */
0218     template<typename Other>
0219     inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
0220     { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
0221
0222     /** \returns the product of a permutation with another inverse permutation.
0223       *
0224       * \note \blank \note_try_to_help_rvo
0225       */
0226     template<typename Other>
0227     inline PlainPermutationType operator*(const InverseImpl<Other,PermutationStorage>& other) const
0228     { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
0229
0230     /** \returns the product of an inverse permutation with another permutation.
0231       *
0232       * \note \blank \note_try_to_help_rvo
0233       */
0234     template<typename Other> friend
0235     inline PlainPermutationType operator*(const InverseImpl<Other, PermutationStorage>& other, const PermutationBase& perm)
0236     { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
0237
0238     /** \returns the determinant of the permutation matrix, which is either 1 or -1 depending on the parity of the permutation.
0239       *
0240       * This function is O(\c n) procedure allocating a buffer of \c n booleans.
0241       */
0242     Index determinant() const
0243     {
0244       Index res = 1;
0245       Index n = size();
0246       Matrix<bool,RowsAtCompileTime,1,0,MaxRowsAtCompileTime> mask(n);
0247       mask.fill(false);
0248       Index r = 0;
0249       while(r < n)
0250       {
0251         // search for the next seed
0252         while(r<n && mask[r]) r++;
0253         if(r>=n)
0254           break;
0255         // we got one, let's follow it until we are back to the seed
0256         Index k0 = r++;
0257         mask.coeffRef(k0) = true;
0258         for(Index k=indices().coeff(k0); k!=k0; k=indices().coeff(k))
0259         {
0260           mask.coeffRef(k) = true;
0261           res = -res;
0262         }
0263       }
0264       return res;
0265     }
0266
0267   protected:
0268
0269 };
0270
0271 namespace internal {
0272 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex>
0273 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> >
0274  : traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
0275 {
0276   typedef PermutationStorage StorageKind;
0277   typedef Matrix<_StorageIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
0278   typedef _StorageIndex StorageIndex;
0279   typedef void Scalar;
0280 };
0281 }
0282
0283 /** \class PermutationMatrix
0284   * \ingroup Core_Module
0285   *
0286   * \brief Permutation matrix
0287   *
0288   * \tparam SizeAtCompileTime the number of rows/cols, or Dynamic
0289   * \tparam 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.
0290   * \tparam _StorageIndex the integer type of the indices
0291   *
0292   * This class represents a permutation matrix, internally stored as a vector of integers.
0293   *
0294   * \sa class PermutationBase, class PermutationWrapper, class DiagonalMatrix
0295   */
0296 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex>
0297 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> >
0298 {
0299     typedef PermutationBase<PermutationMatrix> Base;
0300     typedef internal::traits<PermutationMatrix> Traits;
0301   public:
0302
0303     typedef const PermutationMatrix& Nested;
0304
0305     #ifndef EIGEN_PARSED_BY_DOXYGEN
0306     typedef typename Traits::IndicesType IndicesType;
0307     typedef typename Traits::StorageIndex StorageIndex;
0308     #endif
0309
0310     inline PermutationMatrix()
0311     {}
0312
0313     /** Constructs an uninitialized permutation matrix of given size.
0314       */
0315     explicit inline PermutationMatrix(Index size) : m_indices(size)
0316     {
0317       eigen_internal_assert(size <= NumTraits<StorageIndex>::highest());
0318     }
0319
0320     /** Copy constructor. */
0321     template<typename OtherDerived>
0322     inline PermutationMatrix(const PermutationBase<OtherDerived>& other)
0323       : m_indices(other.indices()) {}
0324
0325     /** Generic constructor from expression of the indices. The indices
0326       * array has the meaning that the permutations sends each integer i to indices[i].
0327       *
0328       * \warning It is your responsibility to check that the indices array that you passes actually
0329       * describes a permutation, i.e., each value between 0 and n-1 occurs exactly once, where n is the
0330       * array's size.
0331       */
0332     template<typename Other>
0333     explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
0334     {}
0335
0336     /** Convert the Transpositions \a tr to a permutation matrix */
0337     template<typename Other>
0338     explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
0339       : m_indices(tr.size())
0340     {
0341       *this = tr;
0342     }
0343
0344     /** Copies the other permutation into *this */
0345     template<typename Other>
0346     PermutationMatrix& operator=(const PermutationBase<Other>& other)
0347     {
0348       m_indices = other.indices();
0349       return *this;
0350     }
0351
0352     /** Assignment from the Transpositions \a tr */
0353     template<typename Other>
0354     PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
0355     {
0356       return Base::operator=(tr.derived());
0357     }
0358
0359     /** const version of indices(). */
0360     const IndicesType& indices() const { return m_indices; }
0361     /** \returns a reference to the stored array representing the permutation. */
0362     IndicesType& indices() { return m_indices; }
0363
0364
0365     /**** multiplication helpers to hopefully get RVO ****/
0366
0367 #ifndef EIGEN_PARSED_BY_DOXYGEN
0368     template<typename Other>
0369     PermutationMatrix(const InverseImpl<Other,PermutationStorage>& other)
0370       : m_indices(other.derived().nestedExpression().size())
0371     {
0372       eigen_internal_assert(m_indices.size() <= NumTraits<StorageIndex>::highest());
0373       StorageIndex end = StorageIndex(m_indices.size());
0374       for (StorageIndex i=0; i<end;++i)
0375         m_indices.coeffRef(other.derived().nestedExpression().indices().coeff(i)) = i;
0376     }
0377     template<typename Lhs,typename Rhs>
0378     PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
0379       : m_indices(lhs.indices().size())
0380     {
0381       Base::assignProduct(lhs,rhs);
0382     }
0383 #endif
0384
0385   protected:
0386
0387     IndicesType m_indices;
0388 };
0389
0390
0391 namespace internal {
0392 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess>
0393 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> >
0394  : traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
0395 {
0396   typedef PermutationStorage StorageKind;
0397   typedef Map<const Matrix<_StorageIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType;
0398   typedef _StorageIndex StorageIndex;
0399   typedef void Scalar;
0400 };
0401 }
0402
0403 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess>
0404 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess>
0405   : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> >
0406 {
0407     typedef PermutationBase<Map> Base;
0408     typedef internal::traits<Map> Traits;
0409   public:
0410
0411     #ifndef EIGEN_PARSED_BY_DOXYGEN
0412     typedef typename Traits::IndicesType IndicesType;
0413     typedef typename IndicesType::Scalar StorageIndex;
0414     #endif
0415
0416     inline Map(const StorageIndex* indicesPtr)
0417       : m_indices(indicesPtr)
0418     {}
0419
0420     inline Map(const StorageIndex* indicesPtr, Index size)
0421       : m_indices(indicesPtr,size)
0422     {}
0423
0424     /** Copies the other permutation into *this */
0425     template<typename Other>
0426     Map& operator=(const PermutationBase<Other>& other)
0427     { return Base::operator=(other.derived()); }
0428
0429     /** Assignment from the Transpositions \a tr */
0430     template<typename Other>
0431     Map& operator=(const TranspositionsBase<Other>& tr)
0432     { return Base::operator=(tr.derived()); }
0433
0434     #ifndef EIGEN_PARSED_BY_DOXYGEN
0435     /** This is a special case of the templated operator=. Its purpose is to
0436       * prevent a default operator= from hiding the templated operator=.
0437       */
0438     Map& operator=(const Map& other)
0439     {
0440       m_indices = other.m_indices;
0441       return *this;
0442     }
0443     #endif
0444
0445     /** const version of indices(). */
0446     const IndicesType& indices() const { return m_indices; }
0447     /** \returns a reference to the stored array representing the permutation. */
0448     IndicesType& indices() { return m_indices; }
0449
0450   protected:
0451
0452     IndicesType m_indices;
0453 };
0454
0455 template<typename _IndicesType> class TranspositionsWrapper;
0456 namespace internal {
0457 template<typename _IndicesType>
0458 struct traits<PermutationWrapper<_IndicesType> >
0459 {
0460   typedef PermutationStorage StorageKind;
0461   typedef void Scalar;
0462   typedef typename _IndicesType::Scalar StorageIndex;
0463   typedef _IndicesType IndicesType;
0464   enum {
0465     RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
0466     ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
0467     MaxRowsAtCompileTime = IndicesType::MaxSizeAtCompileTime,
0468     MaxColsAtCompileTime = IndicesType::MaxSizeAtCompileTime,
0469     Flags = 0
0470   };
0471 };
0472 }
0473
0474 /** \class PermutationWrapper
0475   * \ingroup Core_Module
0476   *
0477   * \brief Class to view a vector of integers as a permutation matrix
0478   *
0479   * \tparam _IndicesType the type of the vector of integer (can be any compatible expression)
0480   *
0481   * This class allows to view any vector expression of integers as a permutation matrix.
0482   *
0483   * \sa class PermutationBase, class PermutationMatrix
0484   */
0485 template<typename _IndicesType>
0486 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
0487 {
0488     typedef PermutationBase<PermutationWrapper> Base;
0489     typedef internal::traits<PermutationWrapper> Traits;
0490   public:
0491
0492     #ifndef EIGEN_PARSED_BY_DOXYGEN
0493     typedef typename Traits::IndicesType IndicesType;
0494     #endif
0495
0496     inline PermutationWrapper(const IndicesType& indices)
0497       : m_indices(indices)
0498     {}
0499
0500     /** const version of indices(). */
0501     const typename internal::remove_all<typename IndicesType::Nested>::type&
0502     indices() const { return m_indices; }
0503
0504   protected:
0505
0506     typename IndicesType::Nested m_indices;
0507 };
0508
0509
0510 /** \returns the matrix with the permutation applied to the columns.
0511   */
0512 template<typename MatrixDerived, typename PermutationDerived>
0513 EIGEN_DEVICE_FUNC
0514 const Product<MatrixDerived, PermutationDerived, AliasFreeProduct>
0515 operator*(const MatrixBase<MatrixDerived> &matrix,
0516           const PermutationBase<PermutationDerived>& permutation)
0517 {
0518   return Product<MatrixDerived, PermutationDerived, AliasFreeProduct>
0519             (matrix.derived(), permutation.derived());
0520 }
0521
0522 /** \returns the matrix with the permutation applied to the rows.
0523   */
0524 template<typename PermutationDerived, typename MatrixDerived>
0525 EIGEN_DEVICE_FUNC
0526 const Product<PermutationDerived, MatrixDerived, AliasFreeProduct>
0527 operator*(const PermutationBase<PermutationDerived> &permutation,
0528           const MatrixBase<MatrixDerived>& matrix)
0529 {
0530   return Product<PermutationDerived, MatrixDerived, AliasFreeProduct>
0531             (permutation.derived(), matrix.derived());
0532 }
0533
0534
0535 template<typename PermutationType>
0536 class InverseImpl<PermutationType, PermutationStorage>
0537   : public EigenBase<Inverse<PermutationType> >
0538 {
0539     typedef typename PermutationType::PlainPermutationType PlainPermutationType;
0540     typedef internal::traits<PermutationType> PermTraits;
0541   protected:
0542     InverseImpl() {}
0543   public:
0544     typedef Inverse<PermutationType> InverseType;
0545     using EigenBase<Inverse<PermutationType> >::derived;
0546
0547     #ifndef EIGEN_PARSED_BY_DOXYGEN
0548     typedef typename PermutationType::DenseMatrixType DenseMatrixType;
0549     enum {
0550       RowsAtCompileTime = PermTraits::RowsAtCompileTime,
0551       ColsAtCompileTime = PermTraits::ColsAtCompileTime,
0552       MaxRowsAtCompileTime = PermTraits::MaxRowsAtCompileTime,
0553       MaxColsAtCompileTime = PermTraits::MaxColsAtCompileTime
0554     };
0555     #endif
0556
0557     #ifndef EIGEN_PARSED_BY_DOXYGEN
0558     template<typename DenseDerived>
0559     void evalTo(MatrixBase<DenseDerived>& other) const
0560     {
0561       other.setZero();
0562       for (Index i=0; i<derived().rows();++i)
0563         other.coeffRef(i, derived().nestedExpression().indices().coeff(i)) = typename DenseDerived::Scalar(1);
0564     }
0565     #endif
0566
0567     /** \return the equivalent permutation matrix */
0568     PlainPermutationType eval() const { return derived(); }
0569
0570     DenseMatrixType toDenseMatrix() const { return derived(); }
0571
0572     /** \returns the matrix with the inverse permutation applied to the columns.
0573       */
0574     template<typename OtherDerived> friend
0575     const Product<OtherDerived, InverseType, AliasFreeProduct>
0576     operator*(const MatrixBase<OtherDerived>& matrix, const InverseType& trPerm)
0577     {
0578       return Product<OtherDerived, InverseType, AliasFreeProduct>(matrix.derived(), trPerm.derived());
0579     }
0580
0581     /** \returns the matrix with the inverse permutation applied to the rows.
0582       */
0583     template<typename OtherDerived>
0584     const Product<InverseType, OtherDerived, AliasFreeProduct>
0585     operator*(const MatrixBase<OtherDerived>& matrix) const
0586     {
0587       return Product<InverseType, OtherDerived, AliasFreeProduct>(derived(), matrix.derived());
0588     }
0589 };
0590
0591 template<typename Derived>
0592 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
0593 {
0594   return derived();
0595 }
0596
0597 namespace internal {
0598
0599 template<> struct AssignmentKind<DenseShape,PermutationShape> { typedef EigenBase2EigenBase Kind; };
0600
0601 } // end namespace internal
0602
0603 } // end namespace Eigen
0604
0605 #endif // EIGEN_PERMUTATIONMATRIX_H