PermutationMatrix.h
Go to the documentation of this file.
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 
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 
00079     template<typename OtherDerived>
00080     Derived& operator=(const PermutationBase<OtherDerived>& other)
00081     {
00082       indices() = other.indices();
00083       return derived();
00084     }
00085 
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 
00100     Derived& operator=(const PermutationBase& other)
00101     {
00102       indices() = other.indices();
00103       return derived();
00104     }
00105     #endif
00106 
00108     inline Index rows() const { return indices().size(); }
00109 
00111     inline Index cols() const { return indices().size(); }
00112 
00114     inline Index size() const { return 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 
00130     DenseMatrixType toDenseMatrix() const
00131     {
00132       return derived();
00133     }
00134 
00136     const IndicesType& indices() const { return derived().indices(); }
00138     IndicesType& indices() { return derived().indices(); }
00139 
00142     inline void resize(Index size)
00143     {
00144       indices().resize(size);
00145     }
00146 
00148     void setIdentity()
00149     {
00150       for(Index i = 0; i < size(); ++i)
00151         indices().coeffRef(i) = i;
00152     }
00153 
00156     void setIdentity(Index size)
00157     {
00158       resize(size);
00159       setIdentity();
00160     }
00161 
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 
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 
00201     inline Transpose<PermutationBase> inverse() const
00202     { return derived(); }
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 
00234     template<typename Other>
00235     inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
00236     { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
00237 
00242     template<typename Other>
00243     inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other) const
00244     { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
00245 
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   protected:
00255 
00256 };
00257 
00272 namespace internal {
00273 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
00274 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
00275  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
00276 {
00277   typedef IndexType Index;
00278   typedef Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
00279 };
00280 }
00281 
00282 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
00283 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
00284 {
00285     typedef PermutationBase<PermutationMatrix> Base;
00286     typedef internal::traits<PermutationMatrix> Traits;
00287   public:
00288 
00289     #ifndef EIGEN_PARSED_BY_DOXYGEN
00290     typedef typename Traits::IndicesType IndicesType;
00291     #endif
00292 
00293     inline PermutationMatrix()
00294     {}
00295 
00298     inline PermutationMatrix(int size) : m_indices(size)
00299     {}
00300 
00302     template<typename OtherDerived>
00303     inline PermutationMatrix(const PermutationBase<OtherDerived>& other)
00304       : m_indices(other.indices()) {}
00305 
00306     #ifndef EIGEN_PARSED_BY_DOXYGEN
00307 
00309     inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
00310     #endif
00311 
00319     template<typename Other>
00320     explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
00321     {}
00322 
00324     template<typename Other>
00325     explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
00326       : m_indices(tr.size())
00327     {
00328       *this = tr;
00329     }
00330 
00332     template<typename Other>
00333     PermutationMatrix& operator=(const PermutationBase<Other>& other)
00334     {
00335       m_indices = other.indices();
00336       return *this;
00337     }
00338 
00340     template<typename Other>
00341     PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
00342     {
00343       return Base::operator=(tr.derived());
00344     }
00345 
00346     #ifndef EIGEN_PARSED_BY_DOXYGEN
00347 
00350     PermutationMatrix& operator=(const PermutationMatrix& other)
00351     {
00352       m_indices = other.m_indices;
00353       return *this;
00354     }
00355     #endif
00356 
00358     const IndicesType& indices() const { return m_indices; }
00360     IndicesType& indices() { return m_indices; }
00361 
00362 
00363     /**** multiplication helpers to hopefully get RVO ****/
00364 
00365 #ifndef EIGEN_PARSED_BY_DOXYGEN
00366     template<typename Other>
00367     PermutationMatrix(const Transpose<PermutationBase<Other> >& other)
00368       : m_indices(other.nestedPermutation().size())
00369     {
00370       for (int i=0; i<m_indices.size();++i) m_indices.coeffRef(other.nestedPermutation().indices().coeff(i)) = i;
00371     }
00372     template<typename Lhs,typename Rhs>
00373     PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
00374       : m_indices(lhs.indices().size())
00375     {
00376       Base::assignProduct(lhs,rhs);
00377     }
00378 #endif
00379 
00380   protected:
00381 
00382     IndicesType m_indices;
00383 };
00384 
00385 
00386 namespace internal {
00387 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
00388 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
00389  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
00390 {
00391   typedef IndexType Index;
00392   typedef Map<const Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType;
00393 };
00394 }
00395 
00396 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
00397 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess>
00398   : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
00399 {
00400     typedef PermutationBase<Map> Base;
00401     typedef internal::traits<Map> Traits;
00402   public:
00403 
00404     #ifndef EIGEN_PARSED_BY_DOXYGEN
00405     typedef typename Traits::IndicesType IndicesType;
00406     typedef typename IndicesType::Scalar Index;
00407     #endif
00408 
00409     inline Map(const Index* indices)
00410       : m_indices(indices)
00411     {}
00412 
00413     inline Map(const Index* indices, Index size)
00414       : m_indices(indices,size)
00415     {}
00416 
00418     template<typename Other>
00419     Map& operator=(const PermutationBase<Other>& other)
00420     { return Base::operator=(other.derived()); }
00421 
00423     template<typename Other>
00424     Map& operator=(const TranspositionsBase<Other>& tr)
00425     { return Base::operator=(tr.derived()); }
00426 
00427     #ifndef EIGEN_PARSED_BY_DOXYGEN
00428 
00431     Map& operator=(const Map& other)
00432     {
00433       m_indices = other.m_indices;
00434       return *this;
00435     }
00436     #endif
00437 
00439     const IndicesType& indices() const { return m_indices; }
00441     IndicesType& indices() { return m_indices; }
00442 
00443   protected:
00444 
00445     IndicesType m_indices;
00446 };
00447 
00460 struct PermutationStorage {};
00461 
00462 template<typename _IndicesType> class TranspositionsWrapper;
00463 namespace internal {
00464 template<typename _IndicesType>
00465 struct traits<PermutationWrapper<_IndicesType> >
00466 {
00467   typedef PermutationStorage StorageKind;
00468   typedef typename _IndicesType::Scalar Scalar;
00469   typedef typename _IndicesType::Scalar Index;
00470   typedef _IndicesType IndicesType;
00471   enum {
00472     RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
00473     ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
00474     MaxRowsAtCompileTime = IndicesType::MaxRowsAtCompileTime,
00475     MaxColsAtCompileTime = IndicesType::MaxColsAtCompileTime,
00476     Flags = 0,
00477     CoeffReadCost = _IndicesType::CoeffReadCost
00478   };
00479 };
00480 }
00481 
00482 template<typename _IndicesType>
00483 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
00484 {
00485     typedef PermutationBase<PermutationWrapper> Base;
00486     typedef internal::traits<PermutationWrapper> Traits;
00487   public:
00488 
00489     #ifndef EIGEN_PARSED_BY_DOXYGEN
00490     typedef typename Traits::IndicesType IndicesType;
00491     #endif
00492 
00493     inline PermutationWrapper(const IndicesType& indices)
00494       : m_indices(indices)
00495     {}
00496 
00498     const typename internal::remove_all<typename IndicesType::Nested>::type&
00499     indices() const { return m_indices; }
00500 
00501   protected:
00502 
00503     typename IndicesType::Nested m_indices;
00504 };
00505 
00508 template<typename Derived, typename PermutationDerived>
00509 inline const internal::permut_matrix_product_retval<PermutationDerived, Derived, OnTheRight>
00510 operator*(const MatrixBase<Derived>& matrix,
00511           const PermutationBase<PermutationDerived> &permutation)
00512 {
00513   return internal::permut_matrix_product_retval
00514            <PermutationDerived, Derived, OnTheRight>
00515            (permutation.derived(), matrix.derived());
00516 }
00517 
00520 template<typename Derived, typename PermutationDerived>
00521 inline const internal::permut_matrix_product_retval
00522                <PermutationDerived, Derived, OnTheLeft>
00523 operator*(const PermutationBase<PermutationDerived> &permutation,
00524           const MatrixBase<Derived>& matrix)
00525 {
00526   return internal::permut_matrix_product_retval
00527            <PermutationDerived, Derived, OnTheLeft>
00528            (permutation.derived(), matrix.derived());
00529 }
00530 
00531 namespace internal {
00532 
00533 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
00534 struct traits<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
00535 {
00536   typedef typename MatrixType::PlainObject ReturnType;
00537 };
00538 
00539 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
00540 struct permut_matrix_product_retval
00541  : public ReturnByValue<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
00542 {
00543     typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
00544 
00545     permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
00546       : m_permutation(perm), m_matrix(matrix)
00547     {}
00548 
00549     inline int rows() const { return m_matrix.rows(); }
00550     inline int cols() const { return m_matrix.cols(); }
00551 
00552     template<typename Dest> inline void evalTo(Dest& dst) const
00553     {
00554       const int n = Side==OnTheLeft ? rows() : cols();
00555 
00556       if(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix))
00557       {
00558         // apply the permutation inplace
00559         Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
00560         mask.fill(false);
00561         int r = 0;
00562         while(r < m_permutation.size())
00563         {
00564           // search for the next seed
00565           while(r<m_permutation.size() && mask[r]) r++;
00566           if(r>=m_permutation.size())
00567             break;
00568           // we got one, let's follow it until we are back to the seed
00569           int k0 = r++;
00570           int kPrev = k0;
00571           mask.coeffRef(k0) = true;
00572           for(int k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
00573           {
00574                   Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
00575             .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
00576                        (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
00577 
00578             mask.coeffRef(k) = true;
00579             kPrev = k;
00580           }
00581         }
00582       }
00583       else
00584       {
00585         for(int i = 0; i < n; ++i)
00586         {
00587           Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
00588                (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
00589 
00590           =
00591 
00592           Block<const MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
00593                (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
00594         }
00595       }
00596     }
00597 
00598   protected:
00599     const PermutationType& m_permutation;
00600     typename MatrixType::Nested m_matrix;
00601 };
00602 
00603 /* Template partial specialization for transposed/inverse permutations */
00604 
00605 template<typename Derived>
00606 struct traits<Transpose<PermutationBase<Derived> > >
00607  : traits<Derived>
00608 {};
00609 
00610 } // end namespace internal
00611 
00612 template<typename Derived>
00613 class Transpose<PermutationBase<Derived> >
00614   : public EigenBase<Transpose<PermutationBase<Derived> > >
00615 {
00616     typedef Derived PermutationType;
00617     typedef typename PermutationType::IndicesType IndicesType;
00618     typedef typename PermutationType::PlainPermutationType PlainPermutationType;
00619   public:
00620 
00621     #ifndef EIGEN_PARSED_BY_DOXYGEN
00622     typedef internal::traits<PermutationType> Traits;
00623     typedef typename Derived::DenseMatrixType DenseMatrixType;
00624     enum {
00625       Flags = Traits::Flags,
00626       CoeffReadCost = Traits::CoeffReadCost,
00627       RowsAtCompileTime = Traits::RowsAtCompileTime,
00628       ColsAtCompileTime = Traits::ColsAtCompileTime,
00629       MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
00630       MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
00631     };
00632     typedef typename Traits::Scalar Scalar;
00633     #endif
00634 
00635     Transpose(const PermutationType& p) : m_permutation(p) {}
00636 
00637     inline int rows() const { return m_permutation.rows(); }
00638     inline int cols() const { return m_permutation.cols(); }
00639 
00640     #ifndef EIGEN_PARSED_BY_DOXYGEN
00641     template<typename DenseDerived>
00642     void evalTo(MatrixBase<DenseDerived>& other) const
00643     {
00644       other.setZero();
00645       for (int i=0; i<rows();++i)
00646         other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1);
00647     }
00648     #endif
00649 
00651     PlainPermutationType eval() const { return *this; }
00652 
00653     DenseMatrixType toDenseMatrix() const { return *this; }
00654 
00657     template<typename OtherDerived> friend
00658     inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>
00659     operator*(const MatrixBase<OtherDerived>& matrix, const Transpose& trPerm)
00660     {
00661       return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>(trPerm.m_permutation, matrix.derived());
00662     }
00663 
00666     template<typename OtherDerived>
00667     inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>
00668     operator*(const MatrixBase<OtherDerived>& matrix) const
00669     {
00670       return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>(m_permutation, matrix.derived());
00671     }
00672 
00673     const PermutationType& nestedPermutation() const { return m_permutation; }
00674 
00675   protected:
00676     const PermutationType& m_permutation;
00677 };
00678 
00679 template<typename Derived>
00680 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
00681 {
00682   return derived();
00683 }
00684 
00685 } // end namespace Eigen
00686 
00687 #endif // EIGEN_PERMUTATIONMATRIX_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:25:44