00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef EIGEN_PERMUTATIONMATRIX_H
00027 #define EIGEN_PERMUTATIONMATRIX_H
00028
00029 namespace Eigen {
00030
00031 template<int RowCol,typename IndicesType,typename MatrixType, typename StorageKind> class PermutedImpl;
00032
00057 namespace internal {
00058
00059 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
00060 struct permut_matrix_product_retval;
00061 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
00062 struct permut_sparsematrix_product_retval;
00063 enum PermPermProduct_t {PermPermProduct};
00064
00065 }
00066
00067 template<typename Derived>
00068 class PermutationBase : public EigenBase<Derived>
00069 {
00070 typedef internal::traits<Derived> Traits;
00071 typedef EigenBase<Derived> Base;
00072 public:
00073
00074 #ifndef EIGEN_PARSED_BY_DOXYGEN
00075 typedef typename Traits::IndicesType IndicesType;
00076 enum {
00077 Flags = Traits::Flags,
00078 CoeffReadCost = Traits::CoeffReadCost,
00079 RowsAtCompileTime = Traits::RowsAtCompileTime,
00080 ColsAtCompileTime = Traits::ColsAtCompileTime,
00081 MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
00082 MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
00083 };
00084 typedef typename Traits::Scalar Scalar;
00085 typedef typename Traits::Index Index;
00086 typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime,0,MaxRowsAtCompileTime,MaxColsAtCompileTime>
00087 DenseMatrixType;
00088 typedef PermutationMatrix<IndicesType::SizeAtCompileTime,IndicesType::MaxSizeAtCompileTime,Index>
00089 PlainPermutationType;
00090 using Base::derived;
00091 #endif
00092
00094 template<typename OtherDerived>
00095 Derived& operator=(const PermutationBase<OtherDerived>& other)
00096 {
00097 indices() = other.indices();
00098 return derived();
00099 }
00100
00102 template<typename OtherDerived>
00103 Derived& operator=(const TranspositionsBase<OtherDerived>& tr)
00104 {
00105 setIdentity(tr.size());
00106 for(Index k=size()-1; k>=0; --k)
00107 applyTranspositionOnTheRight(k,tr.coeff(k));
00108 return derived();
00109 }
00110
00111 #ifndef EIGEN_PARSED_BY_DOXYGEN
00112
00115 Derived& operator=(const PermutationBase& other)
00116 {
00117 indices() = other.indices();
00118 return derived();
00119 }
00120 #endif
00121
00123 inline Index rows() const { return indices().size(); }
00124
00126 inline Index cols() const { return indices().size(); }
00127
00129 inline Index size() const { return indices().size(); }
00130
00131 #ifndef EIGEN_PARSED_BY_DOXYGEN
00132 template<typename DenseDerived>
00133 void evalTo(MatrixBase<DenseDerived>& other) const
00134 {
00135 other.setZero();
00136 for (int i=0; i<rows();++i)
00137 other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
00138 }
00139 #endif
00140
00145 DenseMatrixType toDenseMatrix() const
00146 {
00147 return derived();
00148 }
00149
00151 const IndicesType& indices() const { return derived().indices(); }
00153 IndicesType& indices() { return derived().indices(); }
00154
00157 inline void resize(Index size)
00158 {
00159 indices().resize(size);
00160 }
00161
00163 void setIdentity()
00164 {
00165 for(Index i = 0; i < size(); ++i)
00166 indices().coeffRef(i) = i;
00167 }
00168
00171 void setIdentity(Index size)
00172 {
00173 resize(size);
00174 setIdentity();
00175 }
00176
00186 Derived& applyTranspositionOnTheLeft(Index i, Index j)
00187 {
00188 eigen_assert(i>=0 && j>=0 && i<size() && j<size());
00189 for(Index k = 0; k < size(); ++k)
00190 {
00191 if(indices().coeff(k) == i) indices().coeffRef(k) = j;
00192 else if(indices().coeff(k) == j) indices().coeffRef(k) = i;
00193 }
00194 return derived();
00195 }
00196
00205 Derived& applyTranspositionOnTheRight(Index i, Index j)
00206 {
00207 eigen_assert(i>=0 && j>=0 && i<size() && j<size());
00208 std::swap(indices().coeffRef(i), indices().coeffRef(j));
00209 return derived();
00210 }
00211
00216 inline Transpose<PermutationBase> inverse() const
00217 { return derived(); }
00222 inline Transpose<PermutationBase> transpose() const
00223 { return derived(); }
00224
00225
00226
00227
00228 #ifndef EIGEN_PARSED_BY_DOXYGEN
00229 protected:
00230 template<typename OtherDerived>
00231 void assignTranspose(const PermutationBase<OtherDerived>& other)
00232 {
00233 for (int i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i;
00234 }
00235 template<typename Lhs,typename Rhs>
00236 void assignProduct(const Lhs& lhs, const Rhs& rhs)
00237 {
00238 eigen_assert(lhs.cols() == rhs.rows());
00239 for (int i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
00240 }
00241 #endif
00242
00243 public:
00244
00249 template<typename Other>
00250 inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
00251 { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
00252
00257 template<typename Other>
00258 inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other) const
00259 { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
00260
00265 template<typename Other> friend
00266 inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other, const PermutationBase& perm)
00267 { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
00268
00269 protected:
00270
00271 };
00272
00287 namespace internal {
00288 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
00289 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
00290 : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
00291 {
00292 typedef IndexType Index;
00293 typedef Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
00294 };
00295 }
00296
00297 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
00298 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
00299 {
00300 typedef PermutationBase<PermutationMatrix> Base;
00301 typedef internal::traits<PermutationMatrix> Traits;
00302 public:
00303
00304 #ifndef EIGEN_PARSED_BY_DOXYGEN
00305 typedef typename Traits::IndicesType IndicesType;
00306 #endif
00307
00308 inline PermutationMatrix()
00309 {}
00310
00313 inline PermutationMatrix(int size) : m_indices(size)
00314 {}
00315
00317 template<typename OtherDerived>
00318 inline PermutationMatrix(const PermutationBase<OtherDerived>& other)
00319 : m_indices(other.indices()) {}
00320
00321 #ifndef EIGEN_PARSED_BY_DOXYGEN
00322
00324 inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
00325 #endif
00326
00334 template<typename Other>
00335 explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
00336 {}
00337
00339 template<typename Other>
00340 explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
00341 : m_indices(tr.size())
00342 {
00343 *this = tr;
00344 }
00345
00347 template<typename Other>
00348 PermutationMatrix& operator=(const PermutationBase<Other>& other)
00349 {
00350 m_indices = other.indices();
00351 return *this;
00352 }
00353
00355 template<typename Other>
00356 PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
00357 {
00358 return Base::operator=(tr.derived());
00359 }
00360
00361 #ifndef EIGEN_PARSED_BY_DOXYGEN
00362
00365 PermutationMatrix& operator=(const PermutationMatrix& other)
00366 {
00367 m_indices = other.m_indices;
00368 return *this;
00369 }
00370 #endif
00371
00373 const IndicesType& indices() const { return m_indices; }
00375 IndicesType& indices() { return m_indices; }
00376
00377
00378
00379
00380 #ifndef EIGEN_PARSED_BY_DOXYGEN
00381 template<typename Other>
00382 PermutationMatrix(const Transpose<PermutationBase<Other> >& other)
00383 : m_indices(other.nestedPermutation().size())
00384 {
00385 for (int i=0; i<m_indices.size();++i) m_indices.coeffRef(other.nestedPermutation().indices().coeff(i)) = i;
00386 }
00387 template<typename Lhs,typename Rhs>
00388 PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
00389 : m_indices(lhs.indices().size())
00390 {
00391 Base::assignProduct(lhs,rhs);
00392 }
00393 #endif
00394
00395 protected:
00396
00397 IndicesType m_indices;
00398 };
00399
00400
00401 namespace internal {
00402 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
00403 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
00404 : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
00405 {
00406 typedef IndexType Index;
00407 typedef Map<const Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType;
00408 };
00409 }
00410
00411 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
00412 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess>
00413 : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
00414 {
00415 typedef PermutationBase<Map> Base;
00416 typedef internal::traits<Map> Traits;
00417 public:
00418
00419 #ifndef EIGEN_PARSED_BY_DOXYGEN
00420 typedef typename Traits::IndicesType IndicesType;
00421 typedef typename IndicesType::Scalar Index;
00422 #endif
00423
00424 inline Map(const Index* indices)
00425 : m_indices(indices)
00426 {}
00427
00428 inline Map(const Index* indices, Index size)
00429 : m_indices(indices,size)
00430 {}
00431
00433 template<typename Other>
00434 Map& operator=(const PermutationBase<Other>& other)
00435 { return Base::operator=(other.derived()); }
00436
00438 template<typename Other>
00439 Map& operator=(const TranspositionsBase<Other>& tr)
00440 { return Base::operator=(tr.derived()); }
00441
00442 #ifndef EIGEN_PARSED_BY_DOXYGEN
00443
00446 Map& operator=(const Map& other)
00447 {
00448 m_indices = other.m_indices;
00449 return *this;
00450 }
00451 #endif
00452
00454 const IndicesType& indices() const { return m_indices; }
00456 IndicesType& indices() { return m_indices; }
00457
00458 protected:
00459
00460 IndicesType m_indices;
00461 };
00462
00475 struct PermutationStorage {};
00476
00477 template<typename _IndicesType> class TranspositionsWrapper;
00478 namespace internal {
00479 template<typename _IndicesType>
00480 struct traits<PermutationWrapper<_IndicesType> >
00481 {
00482 typedef PermutationStorage StorageKind;
00483 typedef typename _IndicesType::Scalar Scalar;
00484 typedef typename _IndicesType::Scalar Index;
00485 typedef _IndicesType IndicesType;
00486 enum {
00487 RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
00488 ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
00489 MaxRowsAtCompileTime = IndicesType::MaxRowsAtCompileTime,
00490 MaxColsAtCompileTime = IndicesType::MaxColsAtCompileTime,
00491 Flags = 0,
00492 CoeffReadCost = _IndicesType::CoeffReadCost
00493 };
00494 };
00495 }
00496
00497 template<typename _IndicesType>
00498 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
00499 {
00500 typedef PermutationBase<PermutationWrapper> Base;
00501 typedef internal::traits<PermutationWrapper> Traits;
00502 public:
00503
00504 #ifndef EIGEN_PARSED_BY_DOXYGEN
00505 typedef typename Traits::IndicesType IndicesType;
00506 #endif
00507
00508 inline PermutationWrapper(const IndicesType& indices)
00509 : m_indices(indices)
00510 {}
00511
00513 const typename internal::remove_all<typename IndicesType::Nested>::type&
00514 indices() const { return m_indices; }
00515
00516 protected:
00517
00518 typename IndicesType::Nested m_indices;
00519 };
00520
00523 template<typename Derived, typename PermutationDerived>
00524 inline const internal::permut_matrix_product_retval<PermutationDerived, Derived, OnTheRight>
00525 operator*(const MatrixBase<Derived>& matrix,
00526 const PermutationBase<PermutationDerived> &permutation)
00527 {
00528 return internal::permut_matrix_product_retval
00529 <PermutationDerived, Derived, OnTheRight>
00530 (permutation.derived(), matrix.derived());
00531 }
00532
00535 template<typename Derived, typename PermutationDerived>
00536 inline const internal::permut_matrix_product_retval
00537 <PermutationDerived, Derived, OnTheLeft>
00538 operator*(const PermutationBase<PermutationDerived> &permutation,
00539 const MatrixBase<Derived>& matrix)
00540 {
00541 return internal::permut_matrix_product_retval
00542 <PermutationDerived, Derived, OnTheLeft>
00543 (permutation.derived(), matrix.derived());
00544 }
00545
00546 namespace internal {
00547
00548 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
00549 struct traits<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
00550 {
00551 typedef typename MatrixType::PlainObject ReturnType;
00552 };
00553
00554 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
00555 struct permut_matrix_product_retval
00556 : public ReturnByValue<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
00557 {
00558 typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
00559
00560 permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
00561 : m_permutation(perm), m_matrix(matrix)
00562 {}
00563
00564 inline int rows() const { return m_matrix.rows(); }
00565 inline int cols() const { return m_matrix.cols(); }
00566
00567 template<typename Dest> inline void evalTo(Dest& dst) const
00568 {
00569 const int n = Side==OnTheLeft ? rows() : cols();
00570
00571 if(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix))
00572 {
00573
00574 Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
00575 mask.fill(false);
00576 int r = 0;
00577 while(r < m_permutation.size())
00578 {
00579
00580 while(r<m_permutation.size() && mask[r]) r++;
00581 if(r>=m_permutation.size())
00582 break;
00583
00584 int k0 = r++;
00585 int kPrev = k0;
00586 mask.coeffRef(k0) = true;
00587 for(int k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
00588 {
00589 Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
00590 .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
00591 (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
00592
00593 mask.coeffRef(k) = true;
00594 kPrev = k;
00595 }
00596 }
00597 }
00598 else
00599 {
00600 for(int i = 0; i < n; ++i)
00601 {
00602 Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
00603 (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
00604
00605 =
00606
00607 Block<const MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
00608 (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
00609 }
00610 }
00611 }
00612
00613 protected:
00614 const PermutationType& m_permutation;
00615 typename MatrixType::Nested m_matrix;
00616 };
00617
00618
00619
00620 template<typename Derived>
00621 struct traits<Transpose<PermutationBase<Derived> > >
00622 : traits<Derived>
00623 {};
00624
00625 }
00626
00627 template<typename Derived>
00628 class Transpose<PermutationBase<Derived> >
00629 : public EigenBase<Transpose<PermutationBase<Derived> > >
00630 {
00631 typedef Derived PermutationType;
00632 typedef typename PermutationType::IndicesType IndicesType;
00633 typedef typename PermutationType::PlainPermutationType PlainPermutationType;
00634 public:
00635
00636 #ifndef EIGEN_PARSED_BY_DOXYGEN
00637 typedef internal::traits<PermutationType> Traits;
00638 typedef typename Derived::DenseMatrixType DenseMatrixType;
00639 enum {
00640 Flags = Traits::Flags,
00641 CoeffReadCost = Traits::CoeffReadCost,
00642 RowsAtCompileTime = Traits::RowsAtCompileTime,
00643 ColsAtCompileTime = Traits::ColsAtCompileTime,
00644 MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
00645 MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
00646 };
00647 typedef typename Traits::Scalar Scalar;
00648 #endif
00649
00650 Transpose(const PermutationType& p) : m_permutation(p) {}
00651
00652 inline int rows() const { return m_permutation.rows(); }
00653 inline int cols() const { return m_permutation.cols(); }
00654
00655 #ifndef EIGEN_PARSED_BY_DOXYGEN
00656 template<typename DenseDerived>
00657 void evalTo(MatrixBase<DenseDerived>& other) const
00658 {
00659 other.setZero();
00660 for (int i=0; i<rows();++i)
00661 other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1);
00662 }
00663 #endif
00664
00666 PlainPermutationType eval() const { return *this; }
00667
00668 DenseMatrixType toDenseMatrix() const { return *this; }
00669
00672 template<typename OtherDerived> friend
00673 inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>
00674 operator*(const MatrixBase<OtherDerived>& matrix, const Transpose& trPerm)
00675 {
00676 return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>(trPerm.m_permutation, matrix.derived());
00677 }
00678
00681 template<typename OtherDerived>
00682 inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>
00683 operator*(const MatrixBase<OtherDerived>& matrix) const
00684 {
00685 return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>(m_permutation, matrix.derived());
00686 }
00687
00688 const PermutationType& nestedPermutation() const { return m_permutation; }
00689
00690 protected:
00691 const PermutationType& m_permutation;
00692 };
00693
00694 template<typename Derived>
00695 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
00696 {
00697 return derived();
00698 }
00699
00700 }
00701
00702 #endif // EIGEN_PERMUTATIONMATRIX_H