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 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_PERMUTATIONMATRIX_H
00027 #define EIGEN_PERMUTATIONMATRIX_H
00028 
00029 template<int RowCol,typename IndicesType,typename MatrixType, typename StorageKind> class PermutedImpl;
00030 
00055 namespace internal {
00056 
00057 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
00058 struct permut_matrix_product_retval;
00059 enum PermPermProduct_t {PermPermProduct};
00060 
00061 } // end namespace internal
00062 
00063 template<typename Derived>
00064 class PermutationBase : public EigenBase<Derived>
00065 {
00066     typedef internal::traits<Derived> Traits;
00067     typedef EigenBase<Derived> Base;
00068   public:
00069 
00070     #ifndef EIGEN_PARSED_BY_DOXYGEN
00071     typedef typename Traits::IndicesType IndicesType;
00072     enum {
00073       Flags = Traits::Flags,
00074       CoeffReadCost = Traits::CoeffReadCost,
00075       RowsAtCompileTime = Traits::RowsAtCompileTime,
00076       ColsAtCompileTime = Traits::ColsAtCompileTime,
00077       MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
00078       MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
00079     };
00080     typedef typename Traits::Scalar Scalar;
00081     typedef typename Traits::Index Index;
00082     typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime,0,MaxRowsAtCompileTime,MaxColsAtCompileTime>
00083             DenseMatrixType;
00084     typedef PermutationMatrix<IndicesType::SizeAtCompileTime,IndicesType::MaxSizeAtCompileTime,Index>
00085             PlainPermutationType;
00086     using Base::derived;
00087     #endif
00088 
00090     template<typename OtherDerived>
00091     Derived& operator=(const PermutationBase<OtherDerived>& other)
00092     {
00093       indices() = other.indices();
00094       return derived();
00095     }
00096 
00098     template<typename OtherDerived>
00099     Derived& operator=(const TranspositionsBase<OtherDerived>& tr)
00100     {
00101       setIdentity(tr.size());
00102       for(Index k=size()-1; k>=0; --k)
00103         applyTranspositionOnTheRight(k,tr.coeff(k));
00104       return derived();
00105     }
00106 
00107     #ifndef EIGEN_PARSED_BY_DOXYGEN
00108 
00111     Derived& operator=(const PermutationBase& other)
00112     {
00113       indices() = other.indices();
00114       return derived();
00115     }
00116     #endif
00117 
00119     inline Index rows() const { return indices().size(); }
00120 
00122     inline Index cols() const { return indices().size(); }
00123 
00125     inline Index size() const { return indices().size(); }
00126 
00127     #ifndef EIGEN_PARSED_BY_DOXYGEN
00128     template<typename DenseDerived>
00129     void evalTo(MatrixBase<DenseDerived>& other) const
00130     {
00131       other.setZero();
00132       for (int i=0; i<rows();++i)
00133         other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
00134     }
00135     #endif
00136 
00141     DenseMatrixType toDenseMatrix() const
00142     {
00143       return derived();
00144     }
00145 
00147     const IndicesType& indices() const { return derived().indices(); }
00149     IndicesType& indices() { return derived().indices(); }
00150 
00153     inline void resize(Index size)
00154     {
00155       indices().resize(size);
00156     }
00157 
00159     void setIdentity()
00160     {
00161       for(Index i = 0; i < size(); ++i)
00162         indices().coeffRef(i) = i;
00163     }
00164 
00167     void setIdentity(Index size)
00168     {
00169       resize(size);
00170       setIdentity();
00171     }
00172 
00182     Derived& applyTranspositionOnTheLeft(Index i, Index j)
00183     {
00184       eigen_assert(i>=0 && j>=0 && i<size() && j<size());
00185       for(Index k = 0; k < size(); ++k)
00186       {
00187         if(indices().coeff(k) == i) indices().coeffRef(k) = j;
00188         else if(indices().coeff(k) == j) indices().coeffRef(k) = i;
00189       }
00190       return derived();
00191     }
00192 
00201     Derived& applyTranspositionOnTheRight(Index i, Index j)
00202     {
00203       eigen_assert(i>=0 && j>=0 && i<size() && j<size());
00204       std::swap(indices().coeffRef(i), indices().coeffRef(j));
00205       return derived();
00206     }
00207 
00212     inline Transpose<PermutationBase> inverse() const
00213     { return derived(); }
00218     inline Transpose<PermutationBase> transpose() const
00219     { return derived(); }
00220 
00221     /**** multiplication helpers to hopefully get RVO ****/
00222 
00223   
00224 #ifndef EIGEN_PARSED_BY_DOXYGEN
00225   protected:
00226     template<typename OtherDerived>
00227     void assignTranspose(const PermutationBase<OtherDerived>& other)
00228     {
00229       for (int i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i;
00230     }
00231     template<typename Lhs,typename Rhs>
00232     void assignProduct(const Lhs& lhs, const Rhs& rhs)
00233     {
00234       eigen_assert(lhs.cols() == rhs.rows());
00235       for (int i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
00236     }
00237 #endif
00238 
00239   public:
00240 
00245     template<typename Other>
00246     inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
00247     { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
00248 
00253     template<typename Other>
00254     inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other) const
00255     { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
00256 
00261     template<typename Other> friend
00262     inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other, const PermutationBase& perm)
00263     { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
00264 
00265   protected:
00266 
00267 };
00268 
00283 namespace internal {
00284 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
00285 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
00286  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
00287 {
00288   typedef IndexType Index;
00289   typedef Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
00290 };
00291 }
00292 
00293 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
00294 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
00295 {
00296     typedef PermutationBase<PermutationMatrix> Base;
00297     typedef internal::traits<PermutationMatrix> Traits;
00298   public:
00299 
00300     #ifndef EIGEN_PARSED_BY_DOXYGEN
00301     typedef typename Traits::IndicesType IndicesType;
00302     #endif
00303 
00304     inline PermutationMatrix()
00305     {}
00306 
00309     inline PermutationMatrix(int size) : m_indices(size)
00310     {}
00311 
00313     template<typename OtherDerived>
00314     inline PermutationMatrix(const PermutationBase<OtherDerived>& other)
00315       : m_indices(other.indices()) {}
00316 
00317     #ifndef EIGEN_PARSED_BY_DOXYGEN
00318 
00320     inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
00321     #endif
00322 
00330     template<typename Other>
00331     explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
00332     {}
00333 
00335     template<typename Other>
00336     explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
00337       : m_indices(tr.size())
00338     {
00339       *this = tr;
00340     }
00341 
00343     template<typename Other>
00344     PermutationMatrix& operator=(const PermutationBase<Other>& other)
00345     {
00346       m_indices = other.indices();
00347       return *this;
00348     }
00349 
00351     template<typename Other>
00352     PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
00353     {
00354       return Base::operator=(tr.derived());
00355     }
00356 
00357     #ifndef EIGEN_PARSED_BY_DOXYGEN
00358 
00361     PermutationMatrix& operator=(const PermutationMatrix& other)
00362     {
00363       m_indices = other.m_indices;
00364       return *this;
00365     }
00366     #endif
00367 
00369     const IndicesType& indices() const { return m_indices; }
00371     IndicesType& indices() { return m_indices; }
00372 
00373 
00374     /**** multiplication helpers to hopefully get RVO ****/
00375 
00376 #ifndef EIGEN_PARSED_BY_DOXYGEN
00377     template<typename Other>
00378     PermutationMatrix(const Transpose<PermutationBase<Other> >& other)
00379       : m_indices(other.nestedPermutation().size())
00380     {
00381       for (int i=0; i<m_indices.size();++i) m_indices.coeffRef(other.nestedPermutation().indices().coeff(i)) = i;
00382     }
00383     template<typename Lhs,typename Rhs>
00384     PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
00385       : m_indices(lhs.indices().size())
00386     {
00387       Base::assignProduct(lhs,rhs);
00388     }
00389 #endif
00390 
00391   protected:
00392 
00393     IndicesType m_indices;
00394 };
00395 
00396 
00397 namespace internal {
00398 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
00399 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
00400  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
00401 {
00402   typedef IndexType Index;
00403   typedef Map<const Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType;
00404 };
00405 }
00406 
00407 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
00408 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess>
00409   : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
00410 {
00411     typedef PermutationBase<Map> Base;
00412     typedef internal::traits<Map> Traits;
00413   public:
00414 
00415     #ifndef EIGEN_PARSED_BY_DOXYGEN
00416     typedef typename Traits::IndicesType IndicesType;
00417     typedef typename IndicesType::Scalar Index;
00418     #endif
00419 
00420     inline Map(const Index* indices)
00421       : m_indices(indices)
00422     {}
00423 
00424     inline Map(const Index* indices, Index size)
00425       : m_indices(indices,size)
00426     {}
00427 
00429     template<typename Other>
00430     Map& operator=(const PermutationBase<Other>& other)
00431     { return Base::operator=(other.derived()); }
00432 
00434     template<typename Other>
00435     Map& operator=(const TranspositionsBase<Other>& tr)
00436     { return Base::operator=(tr.derived()); }
00437 
00438     #ifndef EIGEN_PARSED_BY_DOXYGEN
00439 
00442     Map& operator=(const Map& other)
00443     {
00444       m_indices = other.m_indices;
00445       return *this;
00446     }
00447     #endif
00448 
00450     const IndicesType& indices() const { return m_indices; }
00452     IndicesType& indices() { return m_indices; }
00453 
00454   protected:
00455 
00456     IndicesType m_indices;
00457 };
00458 
00471 struct PermutationStorage {};
00472 
00473 template<typename _IndicesType> class TranspositionsWrapper;
00474 namespace internal {
00475 template<typename _IndicesType>
00476 struct traits<PermutationWrapper<_IndicesType> >
00477 {
00478   typedef PermutationStorage StorageKind;
00479   typedef typename _IndicesType::Scalar Scalar;
00480   typedef typename _IndicesType::Scalar Index;
00481   typedef _IndicesType IndicesType;
00482   enum {
00483     RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
00484     ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
00485     MaxRowsAtCompileTime = IndicesType::MaxRowsAtCompileTime,
00486     MaxColsAtCompileTime = IndicesType::MaxColsAtCompileTime,
00487     Flags = 0,
00488     CoeffReadCost = _IndicesType::CoeffReadCost
00489   };
00490 };
00491 }
00492 
00493 template<typename _IndicesType>
00494 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
00495 {
00496     typedef PermutationBase<PermutationWrapper> Base;
00497     typedef internal::traits<PermutationWrapper> Traits;
00498   public:
00499 
00500     #ifndef EIGEN_PARSED_BY_DOXYGEN
00501     typedef typename Traits::IndicesType IndicesType;
00502     #endif
00503 
00504     inline PermutationWrapper(const IndicesType& indices)
00505       : m_indices(indices)
00506     {}
00507 
00509     const typename internal::remove_all<typename IndicesType::Nested>::type&
00510     indices() const { return m_indices; }
00511 
00512   protected:
00513 
00514     const typename IndicesType::Nested m_indices;
00515 };
00516 
00519 template<typename Derived, typename PermutationDerived>
00520 inline const internal::permut_matrix_product_retval<PermutationDerived, Derived, OnTheRight>
00521 operator*(const MatrixBase<Derived>& matrix,
00522           const PermutationBase<PermutationDerived> &permutation)
00523 {
00524   return internal::permut_matrix_product_retval
00525            <PermutationDerived, Derived, OnTheRight>
00526            (permutation.derived(), matrix.derived());
00527 }
00528 
00531 template<typename Derived, typename PermutationDerived>
00532 inline const internal::permut_matrix_product_retval
00533                <PermutationDerived, Derived, OnTheLeft>
00534 operator*(const PermutationBase<PermutationDerived> &permutation,
00535           const MatrixBase<Derived>& matrix)
00536 {
00537   return internal::permut_matrix_product_retval
00538            <PermutationDerived, Derived, OnTheLeft>
00539            (permutation.derived(), matrix.derived());
00540 }
00541 
00542 namespace internal {
00543 
00544 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
00545 struct traits<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
00546 {
00547   typedef typename MatrixType::PlainObject ReturnType;
00548 };
00549 
00550 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
00551 struct permut_matrix_product_retval
00552  : public ReturnByValue<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
00553 {
00554     typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
00555 
00556     permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
00557       : m_permutation(perm), m_matrix(matrix)
00558     {}
00559 
00560     inline int rows() const { return m_matrix.rows(); }
00561     inline int cols() const { return m_matrix.cols(); }
00562 
00563     template<typename Dest> inline void evalTo(Dest& dst) const
00564     {
00565       const int n = Side==OnTheLeft ? rows() : cols();
00566 
00567       if(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix))
00568       {
00569         // apply the permutation inplace
00570         Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
00571         mask.fill(false);
00572         int r = 0;
00573         while(r < m_permutation.size())
00574         {
00575           // search for the next seed
00576           while(r<m_permutation.size() && mask[r]) r++;
00577           if(r>=m_permutation.size())
00578             break;
00579           // we got one, let's follow it until we are back to the seed
00580           int k0 = r++;
00581           int kPrev = k0;
00582           mask.coeffRef(k0) = true;
00583           for(int k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
00584           {
00585                   Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
00586             .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
00587                        (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
00588 
00589             mask.coeffRef(k) = true;
00590             kPrev = k;
00591           }
00592         }
00593       }
00594       else
00595       {
00596         for(int i = 0; i < n; ++i)
00597         {
00598           Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
00599                (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
00600 
00601           =
00602 
00603           Block<const MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
00604                (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
00605         }
00606       }
00607     }
00608 
00609   protected:
00610     const PermutationType& m_permutation;
00611     const typename MatrixType::Nested m_matrix;
00612 };
00613 
00614 /* Template partial specialization for transposed/inverse permutations */
00615 
00616 template<typename Derived>
00617 struct traits<Transpose<PermutationBase<Derived> > >
00618  : traits<Derived>
00619 {};
00620 
00621 } // end namespace internal
00622 
00623 template<typename Derived>
00624 class Transpose<PermutationBase<Derived> >
00625   : public EigenBase<Transpose<PermutationBase<Derived> > >
00626 {
00627     typedef Derived PermutationType;
00628     typedef typename PermutationType::IndicesType IndicesType;
00629     typedef typename PermutationType::PlainPermutationType PlainPermutationType;
00630   public:
00631 
00632     #ifndef EIGEN_PARSED_BY_DOXYGEN
00633     typedef internal::traits<PermutationType> Traits;
00634     typedef typename Derived::DenseMatrixType DenseMatrixType;
00635     enum {
00636       Flags = Traits::Flags,
00637       CoeffReadCost = Traits::CoeffReadCost,
00638       RowsAtCompileTime = Traits::RowsAtCompileTime,
00639       ColsAtCompileTime = Traits::ColsAtCompileTime,
00640       MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
00641       MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
00642     };
00643     typedef typename Traits::Scalar Scalar;
00644     #endif
00645 
00646     Transpose(const PermutationType& p) : m_permutation(p) {}
00647 
00648     inline int rows() const { return m_permutation.rows(); }
00649     inline int cols() const { return m_permutation.cols(); }
00650 
00651     #ifndef EIGEN_PARSED_BY_DOXYGEN
00652     template<typename DenseDerived>
00653     void evalTo(MatrixBase<DenseDerived>& other) const
00654     {
00655       other.setZero();
00656       for (int i=0; i<rows();++i)
00657         other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1);
00658     }
00659     #endif
00660 
00662     PlainPermutationType eval() const { return *this; }
00663 
00664     DenseMatrixType toDenseMatrix() const { return *this; }
00665 
00668     template<typename OtherDerived> friend
00669     inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>
00670     operator*(const MatrixBase<OtherDerived>& matrix, const Transpose& trPerm)
00671     {
00672       return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>(trPerm.m_permutation, matrix.derived());
00673     }
00674 
00677     template<typename OtherDerived>
00678     inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>
00679     operator*(const MatrixBase<OtherDerived>& matrix) const
00680     {
00681       return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>(m_permutation, matrix.derived());
00682     }
00683 
00684     const PermutationType& nestedPermutation() const { return m_permutation; }
00685 
00686   protected:
00687     const PermutationType& m_permutation;
00688 };
00689 
00690 template<typename Derived>
00691 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
00692 {
00693   return derived();
00694 }
00695 
00696 #endif // EIGEN_PERMUTATIONMATRIX_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:11