// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2009 Benoit Jacob // Copyright (C) 2009-2011 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #ifndef EIGEN_PERMUTATIONMATRIX_H #define EIGEN_PERMUTATIONMATRIX_H namespace Eigen { template class PermutedImpl; /** \class PermutationBase * \ingroup Core_Module * * \brief Base class for permutations * * \param Derived the derived class * * This class is the base class for all expressions representing a permutation matrix, * internally stored as a vector of integers. * The convention followed here is that if \f$\sigma \f$ is a permutation, the corresponding permutation matrix * \f$P_\sigma \f$ is such that if \f$(e_1,\ldots,e_p) \f$ is the canonical basis, we have: * \f[ P_\sigma(e_i) = e_{\sigma(i)}. \f] * This convention ensures that for any two permutations \f$\sigma, \tau \f$, we have: * \f[ P_{\sigma\circ\tau} = P_\sigma P_\tau. \f] * * Permutation matrices are square and invertible. * * Notice that in addition to the member functions and operators listed here, there also are non-member * operator* to multiply any kind of permutation object with any kind of matrix expression (MatrixBase) * on either side. * * \sa class PermutationMatrix, class PermutationWrapper */ namespace internal { template struct permut_matrix_product_retval; template struct permut_sparsematrix_product_retval; enum PermPermProduct_t {PermPermProduct}; } // end namespace internal template class PermutationBase : public EigenBase { typedef internal::traits Traits; typedef EigenBase Base; public: #ifndef EIGEN_PARSED_BY_DOXYGEN typedef typename Traits::IndicesType IndicesType; enum { Flags = Traits::Flags, CoeffReadCost = Traits::CoeffReadCost, RowsAtCompileTime = Traits::RowsAtCompileTime, ColsAtCompileTime = Traits::ColsAtCompileTime, MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime, MaxColsAtCompileTime = Traits::MaxColsAtCompileTime }; typedef typename Traits::Scalar Scalar; typedef typename Traits::Index Index; typedef Matrix DenseMatrixType; typedef PermutationMatrix PlainPermutationType; using Base::derived; #endif /** Copies the other permutation into *this */ template Derived& operator=(const PermutationBase& other) { indices() = other.indices(); return derived(); } /** Assignment from the Transpositions \a tr */ template Derived& operator=(const TranspositionsBase& tr) { setIdentity(tr.size()); for(Index k=size()-1; k>=0; --k) applyTranspositionOnTheRight(k,tr.coeff(k)); return derived(); } #ifndef EIGEN_PARSED_BY_DOXYGEN /** This is a special case of the templated operator=. Its purpose is to * prevent a default operator= from hiding the templated operator=. */ Derived& operator=(const PermutationBase& other) { indices() = other.indices(); return derived(); } #endif /** \returns the number of rows */ inline Index rows() const { return Index(indices().size()); } /** \returns the number of columns */ inline Index cols() const { return Index(indices().size()); } /** \returns the size of a side of the respective square matrix, i.e., the number of indices */ inline Index size() const { return Index(indices().size()); } #ifndef EIGEN_PARSED_BY_DOXYGEN template void evalTo(MatrixBase& other) const { other.setZero(); for (int i=0; i=0 && j>=0 && i=0 && j>=0 && i inverse() const { return derived(); } /** \returns the tranpose permutation matrix. * * \note \note_try_to_help_rvo */ inline Transpose transpose() const { return derived(); } /**** multiplication helpers to hopefully get RVO ****/ #ifndef EIGEN_PARSED_BY_DOXYGEN protected: template void assignTranspose(const PermutationBase& other) { for (int i=0; i void assignProduct(const Lhs& lhs, const Rhs& rhs) { eigen_assert(lhs.cols() == rhs.rows()); for (int i=0; i inline PlainPermutationType operator*(const PermutationBase& other) const { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); } /** \returns the product of a permutation with another inverse permutation. * * \note \note_try_to_help_rvo */ template inline PlainPermutationType operator*(const Transpose >& other) const { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); } /** \returns the product of an inverse permutation with another permutation. * * \note \note_try_to_help_rvo */ template friend inline PlainPermutationType operator*(const Transpose >& other, const PermutationBase& perm) { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); } /** \returns the determinant of the permutation matrix, which is either 1 or -1 depending on the parity of the permutation. * * This function is O(\c n) procedure allocating a buffer of \c n booleans. */ Index determinant() const { Index res = 1; Index n = size(); Matrix mask(n); mask.fill(false); Index r = 0; while(r < n) { // search for the next seed while(r=n) break; // we got one, let's follow it until we are back to the seed Index k0 = r++; mask.coeffRef(k0) = true; for(Index k=indices().coeff(k0); k!=k0; k=indices().coeff(k)) { mask.coeffRef(k) = true; res = -res; } } return res; } protected: }; /** \class PermutationMatrix * \ingroup Core_Module * * \brief Permutation matrix * * \param SizeAtCompileTime the number of rows/cols, or Dynamic * \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. * \param IndexType the interger type of the indices * * This class represents a permutation matrix, internally stored as a vector of integers. * * \sa class PermutationBase, class PermutationWrapper, class DiagonalMatrix */ namespace internal { template struct traits > : traits > { typedef IndexType Index; typedef Matrix IndicesType; }; } template class PermutationMatrix : public PermutationBase > { typedef PermutationBase Base; typedef internal::traits Traits; public: #ifndef EIGEN_PARSED_BY_DOXYGEN typedef typename Traits::IndicesType IndicesType; #endif inline PermutationMatrix() {} /** Constructs an uninitialized permutation matrix of given size. */ inline PermutationMatrix(int size) : m_indices(size) {} /** Copy constructor. */ template inline PermutationMatrix(const PermutationBase& other) : m_indices(other.indices()) {} #ifndef EIGEN_PARSED_BY_DOXYGEN /** Standard copy constructor. Defined only to prevent a default copy constructor * from hiding the other templated constructor */ inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {} #endif /** Generic constructor from expression of the indices. The indices * array has the meaning that the permutations sends each integer i to indices[i]. * * \warning It is your responsibility to check that the indices array that you passes actually * describes a permutation, i.e., each value between 0 and n-1 occurs exactly once, where n is the * array's size. */ template explicit inline PermutationMatrix(const MatrixBase& a_indices) : m_indices(a_indices) {} /** Convert the Transpositions \a tr to a permutation matrix */ template explicit PermutationMatrix(const TranspositionsBase& tr) : m_indices(tr.size()) { *this = tr; } /** Copies the other permutation into *this */ template PermutationMatrix& operator=(const PermutationBase& other) { m_indices = other.indices(); return *this; } /** Assignment from the Transpositions \a tr */ template PermutationMatrix& operator=(const TranspositionsBase& tr) { return Base::operator=(tr.derived()); } #ifndef EIGEN_PARSED_BY_DOXYGEN /** This is a special case of the templated operator=. Its purpose is to * prevent a default operator= from hiding the templated operator=. */ PermutationMatrix& operator=(const PermutationMatrix& other) { m_indices = other.m_indices; return *this; } #endif /** const version of indices(). */ const IndicesType& indices() const { return m_indices; } /** \returns a reference to the stored array representing the permutation. */ IndicesType& indices() { return m_indices; } /**** multiplication helpers to hopefully get RVO ****/ #ifndef EIGEN_PARSED_BY_DOXYGEN template PermutationMatrix(const Transpose >& other) : m_indices(other.nestedPermutation().size()) { for (int i=0; i PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs) : m_indices(lhs.indices().size()) { Base::assignProduct(lhs,rhs); } #endif protected: IndicesType m_indices; }; namespace internal { template struct traits,_PacketAccess> > : traits > { typedef IndexType Index; typedef Map, _PacketAccess> IndicesType; }; } template class Map,_PacketAccess> : public PermutationBase,_PacketAccess> > { typedef PermutationBase Base; typedef internal::traits Traits; public: #ifndef EIGEN_PARSED_BY_DOXYGEN typedef typename Traits::IndicesType IndicesType; typedef typename IndicesType::Scalar Index; #endif inline Map(const Index* indicesPtr) : m_indices(indicesPtr) {} inline Map(const Index* indicesPtr, Index size) : m_indices(indicesPtr,size) {} /** Copies the other permutation into *this */ template Map& operator=(const PermutationBase& other) { return Base::operator=(other.derived()); } /** Assignment from the Transpositions \a tr */ template Map& operator=(const TranspositionsBase& tr) { return Base::operator=(tr.derived()); } #ifndef EIGEN_PARSED_BY_DOXYGEN /** This is a special case of the templated operator=. Its purpose is to * prevent a default operator= from hiding the templated operator=. */ Map& operator=(const Map& other) { m_indices = other.m_indices; return *this; } #endif /** const version of indices(). */ const IndicesType& indices() const { return m_indices; } /** \returns a reference to the stored array representing the permutation. */ IndicesType& indices() { return m_indices; } protected: IndicesType m_indices; }; /** \class PermutationWrapper * \ingroup Core_Module * * \brief Class to view a vector of integers as a permutation matrix * * \param _IndicesType the type of the vector of integer (can be any compatible expression) * * This class allows to view any vector expression of integers as a permutation matrix. * * \sa class PermutationBase, class PermutationMatrix */ struct PermutationStorage {}; template class TranspositionsWrapper; namespace internal { template struct traits > { typedef PermutationStorage StorageKind; typedef typename _IndicesType::Scalar Scalar; typedef typename _IndicesType::Scalar Index; typedef _IndicesType IndicesType; enum { RowsAtCompileTime = _IndicesType::SizeAtCompileTime, ColsAtCompileTime = _IndicesType::SizeAtCompileTime, MaxRowsAtCompileTime = IndicesType::MaxRowsAtCompileTime, MaxColsAtCompileTime = IndicesType::MaxColsAtCompileTime, Flags = 0, CoeffReadCost = _IndicesType::CoeffReadCost }; }; } template class PermutationWrapper : public PermutationBase > { typedef PermutationBase Base; typedef internal::traits Traits; public: #ifndef EIGEN_PARSED_BY_DOXYGEN typedef typename Traits::IndicesType IndicesType; #endif inline PermutationWrapper(const IndicesType& a_indices) : m_indices(a_indices) {} /** const version of indices(). */ const typename internal::remove_all::type& indices() const { return m_indices; } protected: typename IndicesType::Nested m_indices; }; /** \returns the matrix with the permutation applied to the columns. */ template inline const internal::permut_matrix_product_retval operator*(const MatrixBase& matrix, const PermutationBase &permutation) { return internal::permut_matrix_product_retval (permutation.derived(), matrix.derived()); } /** \returns the matrix with the permutation applied to the rows. */ template inline const internal::permut_matrix_product_retval operator*(const PermutationBase &permutation, const MatrixBase& matrix) { return internal::permut_matrix_product_retval (permutation.derived(), matrix.derived()); } namespace internal { template struct traits > { typedef typename MatrixType::PlainObject ReturnType; }; template struct permut_matrix_product_retval : public ReturnByValue > { typedef typename remove_all::type MatrixTypeNestedCleaned; typedef typename MatrixType::Index Index; permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix) : m_permutation(perm), m_matrix(matrix) {} inline Index rows() const { return m_matrix.rows(); } inline Index cols() const { return m_matrix.cols(); } template inline void evalTo(Dest& dst) const { const Index n = Side==OnTheLeft ? rows() : cols(); // FIXME we need an is_same for expression that is not sensitive to constness. For instance // is_same_xpr, Block >::value should be true. if( is_same::value && blas_traits::HasUsableDirectAccess && blas_traits::HasUsableDirectAccess && extract_data(dst) == extract_data(m_matrix)) { // apply the permutation inplace Matrix mask(m_permutation.size()); mask.fill(false); Index r = 0; while(r < m_permutation.size()) { // search for the next seed while(r=m_permutation.size()) break; // we got one, let's follow it until we are back to the seed Index k0 = r++; Index kPrev = k0; mask.coeffRef(k0) = true; for(Index k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k)) { Block(dst, k) .swap(Block (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev)); mask.coeffRef(k) = true; kPrev = k; } } } else { for(int i = 0; i < n; ++i) { Block (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i) = Block (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i); } } } protected: const PermutationType& m_permutation; typename MatrixType::Nested m_matrix; }; /* Template partial specialization for transposed/inverse permutations */ template struct traits > > : traits {}; } // end namespace internal template class Transpose > : public EigenBase > > { typedef Derived PermutationType; typedef typename PermutationType::IndicesType IndicesType; typedef typename PermutationType::PlainPermutationType PlainPermutationType; public: #ifndef EIGEN_PARSED_BY_DOXYGEN typedef internal::traits Traits; typedef typename Derived::DenseMatrixType DenseMatrixType; enum { Flags = Traits::Flags, CoeffReadCost = Traits::CoeffReadCost, RowsAtCompileTime = Traits::RowsAtCompileTime, ColsAtCompileTime = Traits::ColsAtCompileTime, MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime, MaxColsAtCompileTime = Traits::MaxColsAtCompileTime }; typedef typename Traits::Scalar Scalar; #endif Transpose(const PermutationType& p) : m_permutation(p) {} inline int rows() const { return m_permutation.rows(); } inline int cols() const { return m_permutation.cols(); } #ifndef EIGEN_PARSED_BY_DOXYGEN template void evalTo(MatrixBase& other) const { other.setZero(); for (int i=0; i friend inline const internal::permut_matrix_product_retval operator*(const MatrixBase& matrix, const Transpose& trPerm) { return internal::permut_matrix_product_retval(trPerm.m_permutation, matrix.derived()); } /** \returns the matrix with the inverse permutation applied to the rows. */ template inline const internal::permut_matrix_product_retval operator*(const MatrixBase& matrix) const { return internal::permut_matrix_product_retval(m_permutation, matrix.derived()); } const PermutationType& nestedPermutation() const { return m_permutation; } protected: const PermutationType& m_permutation; }; template const PermutationWrapper MatrixBase::asPermutation() const { return derived(); } } // end namespace Eigen #endif // EIGEN_PERMUTATIONMATRIX_H