blob: 73a73006b45c414bb2b4de441d146ed0ccb24a7d [file] [log] [blame]
 // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2009 Benoit Jacob // Copyright (C) 2009-2015 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 #include "./InternalHeaderCheck.h" namespace Eigen { namespace internal { enum PermPermProduct_t {PermPermProduct}; } // end namespace internal /** \class PermutationBase * \ingroup Core_Module * * \brief Base class for permutations * * \tparam 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 */ 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, RowsAtCompileTime = Traits::RowsAtCompileTime, ColsAtCompileTime = Traits::ColsAtCompileTime, MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime, MaxColsAtCompileTime = Traits::MaxColsAtCompileTime }; typedef typename Traits::StorageIndex StorageIndex; typedef Matrix DenseMatrixType; typedef PermutationMatrix PlainPermutationType; typedef PlainPermutationType PlainObject; using Base::derived; typedef Inverse InverseReturnType; typedef void Scalar; #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(); } /** \returns the number of rows */ inline EIGEN_DEVICE_FUNC Index rows() const { return Index(indices().size()); } /** \returns the number of columns */ inline EIGEN_DEVICE_FUNC 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 EIGEN_DEVICE_FUNC Index size() const { return Index(indices().size()); } #ifndef EIGEN_PARSED_BY_DOXYGEN template void evalTo(MatrixBase& other) const { other.setZero(); for (Index i=0; i=0 && j>=0 && i=0 && j>=0 && i void assignTranspose(const PermutationBase& other) { for (Index i=0; i void assignProduct(const Lhs& lhs, const Rhs& rhs) { eigen_assert(lhs.cols() == rhs.rows()); for (Index 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 \blank \note_try_to_help_rvo */ template inline PlainPermutationType operator*(const InverseImpl& other) const { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); } /** \returns the product of an inverse permutation with another permutation. * * \note \blank \note_try_to_help_rvo */ template friend inline PlainPermutationType operator*(const InverseImpl& 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: }; namespace internal { template struct traits > : traits > { typedef PermutationStorage StorageKind; typedef Matrix IndicesType; typedef StorageIndex_ StorageIndex; typedef void Scalar; }; } /** \class PermutationMatrix * \ingroup Core_Module * * \brief Permutation matrix * * \tparam SizeAtCompileTime the number of rows/cols, or Dynamic * \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. * \tparam StorageIndex_ the integer type of the indices * * This class represents a permutation matrix, internally stored as a vector of integers. * * \sa class PermutationBase, class PermutationWrapper, class DiagonalMatrix */ template class PermutationMatrix : public PermutationBase > { typedef PermutationBase Base; typedef internal::traits Traits; public: typedef const PermutationMatrix& Nested; #ifndef EIGEN_PARSED_BY_DOXYGEN typedef typename Traits::IndicesType IndicesType; typedef typename Traits::StorageIndex StorageIndex; #endif inline PermutationMatrix() {} /** Constructs an uninitialized permutation matrix of given size. */ explicit inline PermutationMatrix(Index size) : m_indices(size) { eigen_internal_assert(size <= NumTraits::highest()); } /** Copy constructor. */ template inline PermutationMatrix(const PermutationBase& other) : m_indices(other.indices()) {} /** 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& indices) : m_indices(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()); } /** 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 InverseImpl& other) : m_indices(other.derived().nestedExpression().size()) { eigen_internal_assert(m_indices.size() <= NumTraits::highest()); StorageIndex end = StorageIndex(m_indices.size()); for (StorageIndex 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 PermutationStorage StorageKind; typedef Map, PacketAccess_> IndicesType; typedef StorageIndex_ StorageIndex; typedef void Scalar; }; } 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 StorageIndex; #endif inline Map(const StorageIndex* indicesPtr) : m_indices(indicesPtr) {} inline Map(const StorageIndex* 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; }; template class TranspositionsWrapper; namespace internal { template struct traits > { typedef PermutationStorage StorageKind; typedef void Scalar; typedef typename IndicesType_::Scalar StorageIndex; typedef IndicesType_ IndicesType; enum { RowsAtCompileTime = IndicesType_::SizeAtCompileTime, ColsAtCompileTime = IndicesType_::SizeAtCompileTime, MaxRowsAtCompileTime = IndicesType::MaxSizeAtCompileTime, MaxColsAtCompileTime = IndicesType::MaxSizeAtCompileTime, Flags = 0 }; }; } /** \class PermutationWrapper * \ingroup Core_Module * * \brief Class to view a vector of integers as a permutation matrix * * \tparam 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 */ 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& indices) : m_indices(indices) {} /** const version of indices(). */ const internal::remove_all_t& indices() const { return m_indices; } protected: typename IndicesType::Nested m_indices; }; /** \returns the matrix with the permutation applied to the columns. */ template EIGEN_DEVICE_FUNC const Product operator*(const MatrixBase &matrix, const PermutationBase& permutation) { return Product (matrix.derived(), permutation.derived()); } /** \returns the matrix with the permutation applied to the rows. */ template EIGEN_DEVICE_FUNC const Product operator*(const PermutationBase &permutation, const MatrixBase& matrix) { return Product (permutation.derived(), matrix.derived()); } template class InverseImpl : public EigenBase > { typedef typename PermutationType::PlainPermutationType PlainPermutationType; typedef internal::traits PermTraits; protected: InverseImpl() {} public: typedef Inverse InverseType; using EigenBase >::derived; #ifndef EIGEN_PARSED_BY_DOXYGEN typedef typename PermutationType::DenseMatrixType DenseMatrixType; enum { RowsAtCompileTime = PermTraits::RowsAtCompileTime, ColsAtCompileTime = PermTraits::ColsAtCompileTime, MaxRowsAtCompileTime = PermTraits::MaxRowsAtCompileTime, MaxColsAtCompileTime = PermTraits::MaxColsAtCompileTime }; #endif #ifndef EIGEN_PARSED_BY_DOXYGEN template void evalTo(MatrixBase& other) const { other.setZero(); for (Index i=0; i friend const Product operator*(const MatrixBase& matrix, const InverseType& trPerm) { return Product(matrix.derived(), trPerm.derived()); } /** \returns the matrix with the inverse permutation applied to the rows. */ template const Product operator*(const MatrixBase& matrix) const { return Product(derived(), matrix.derived()); } }; template const PermutationWrapper MatrixBase::asPermutation() const { return derived(); } namespace internal { template<> struct AssignmentKind { typedef EigenBase2EigenBase Kind; }; } // end namespace internal } // end namespace Eigen #endif // EIGEN_PERMUTATIONMATRIX_H