$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 Index(indices().size()); } 00109 00111 inline Index cols() const { return Index(indices().size()); } 00112 00114 inline Index size() const { return Index(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 newSize) 00143 { 00144 indices().resize(newSize); 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 newSize) 00157 { 00158 resize(newSize); 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 00258 Index determinant() const 00259 { 00260 Index res = 1; 00261 Index n = size(); 00262 Matrix<bool,RowsAtCompileTime,1,0,MaxRowsAtCompileTime> mask(n); 00263 mask.fill(false); 00264 Index r = 0; 00265 while(r < n) 00266 { 00267 // search for the next seed 00268 while(r<n && mask[r]) r++; 00269 if(r>=n) 00270 break; 00271 // we got one, let's follow it until we are back to the seed 00272 Index k0 = r++; 00273 mask.coeffRef(k0) = true; 00274 for(Index k=indices().coeff(k0); k!=k0; k=indices().coeff(k)) 00275 { 00276 mask.coeffRef(k) = true; 00277 res = -res; 00278 } 00279 } 00280 return res; 00281 } 00282 00283 protected: 00284 00285 }; 00286 00301 namespace internal { 00302 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType> 00303 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> > 00304 : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> > 00305 { 00306 typedef IndexType Index; 00307 typedef Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType; 00308 }; 00309 } 00310 00311 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType> 00312 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> > 00313 { 00314 typedef PermutationBase<PermutationMatrix> Base; 00315 typedef internal::traits<PermutationMatrix> Traits; 00316 public: 00317 00318 #ifndef EIGEN_PARSED_BY_DOXYGEN 00319 typedef typename Traits::IndicesType IndicesType; 00320 #endif 00321 00322 inline PermutationMatrix() 00323 {} 00324 00327 inline PermutationMatrix(int size) : m_indices(size) 00328 {} 00329 00331 template<typename OtherDerived> 00332 inline PermutationMatrix(const PermutationBase<OtherDerived>& other) 00333 : m_indices(other.indices()) {} 00334 00335 #ifndef EIGEN_PARSED_BY_DOXYGEN 00336 00338 inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {} 00339 #endif 00340 00348 template<typename Other> 00349 explicit inline PermutationMatrix(const MatrixBase<Other>& a_indices) : m_indices(a_indices) 00350 {} 00351 00353 template<typename Other> 00354 explicit PermutationMatrix(const TranspositionsBase<Other>& tr) 00355 : m_indices(tr.size()) 00356 { 00357 *this = tr; 00358 } 00359 00361 template<typename Other> 00362 PermutationMatrix& operator=(const PermutationBase<Other>& other) 00363 { 00364 m_indices = other.indices(); 00365 return *this; 00366 } 00367 00369 template<typename Other> 00370 PermutationMatrix& operator=(const TranspositionsBase<Other>& tr) 00371 { 00372 return Base::operator=(tr.derived()); 00373 } 00374 00375 #ifndef EIGEN_PARSED_BY_DOXYGEN 00376 00379 PermutationMatrix& operator=(const PermutationMatrix& other) 00380 { 00381 m_indices = other.m_indices; 00382 return *this; 00383 } 00384 #endif 00385 00387 const IndicesType& indices() const { return m_indices; } 00389 IndicesType& indices() { return m_indices; } 00390 00391 00392 /**** multiplication helpers to hopefully get RVO ****/ 00393 00394 #ifndef EIGEN_PARSED_BY_DOXYGEN 00395 template<typename Other> 00396 PermutationMatrix(const Transpose<PermutationBase<Other> >& other) 00397 : m_indices(other.nestedPermutation().size()) 00398 { 00399 for (int i=0; i<m_indices.size();++i) m_indices.coeffRef(other.nestedPermutation().indices().coeff(i)) = i; 00400 } 00401 template<typename Lhs,typename Rhs> 00402 PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs) 00403 : m_indices(lhs.indices().size()) 00404 { 00405 Base::assignProduct(lhs,rhs); 00406 } 00407 #endif 00408 00409 protected: 00410 00411 IndicesType m_indices; 00412 }; 00413 00414 00415 namespace internal { 00416 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess> 00417 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> > 00418 : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> > 00419 { 00420 typedef IndexType Index; 00421 typedef Map<const Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType; 00422 }; 00423 } 00424 00425 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess> 00426 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> 00427 : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> > 00428 { 00429 typedef PermutationBase<Map> Base; 00430 typedef internal::traits<Map> Traits; 00431 public: 00432 00433 #ifndef EIGEN_PARSED_BY_DOXYGEN 00434 typedef typename Traits::IndicesType IndicesType; 00435 typedef typename IndicesType::Scalar Index; 00436 #endif 00437 00438 inline Map(const Index* indicesPtr) 00439 : m_indices(indicesPtr) 00440 {} 00441 00442 inline Map(const Index* indicesPtr, Index size) 00443 : m_indices(indicesPtr,size) 00444 {} 00445 00447 template<typename Other> 00448 Map& operator=(const PermutationBase<Other>& other) 00449 { return Base::operator=(other.derived()); } 00450 00452 template<typename Other> 00453 Map& operator=(const TranspositionsBase<Other>& tr) 00454 { return Base::operator=(tr.derived()); } 00455 00456 #ifndef EIGEN_PARSED_BY_DOXYGEN 00457 00460 Map& operator=(const Map& other) 00461 { 00462 m_indices = other.m_indices; 00463 return *this; 00464 } 00465 #endif 00466 00468 const IndicesType& indices() const { return m_indices; } 00470 IndicesType& indices() { return m_indices; } 00471 00472 protected: 00473 00474 IndicesType m_indices; 00475 }; 00476 00489 struct PermutationStorage {}; 00490 00491 template<typename _IndicesType> class TranspositionsWrapper; 00492 namespace internal { 00493 template<typename _IndicesType> 00494 struct traits<PermutationWrapper<_IndicesType> > 00495 { 00496 typedef PermutationStorage StorageKind; 00497 typedef typename _IndicesType::Scalar Scalar; 00498 typedef typename _IndicesType::Scalar Index; 00499 typedef _IndicesType IndicesType; 00500 enum { 00501 RowsAtCompileTime = _IndicesType::SizeAtCompileTime, 00502 ColsAtCompileTime = _IndicesType::SizeAtCompileTime, 00503 MaxRowsAtCompileTime = IndicesType::MaxRowsAtCompileTime, 00504 MaxColsAtCompileTime = IndicesType::MaxColsAtCompileTime, 00505 Flags = 0, 00506 CoeffReadCost = _IndicesType::CoeffReadCost 00507 }; 00508 }; 00509 } 00510 00511 template<typename _IndicesType> 00512 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> > 00513 { 00514 typedef PermutationBase<PermutationWrapper> Base; 00515 typedef internal::traits<PermutationWrapper> Traits; 00516 public: 00517 00518 #ifndef EIGEN_PARSED_BY_DOXYGEN 00519 typedef typename Traits::IndicesType IndicesType; 00520 #endif 00521 00522 inline PermutationWrapper(const IndicesType& a_indices) 00523 : m_indices(a_indices) 00524 {} 00525 00527 const typename internal::remove_all<typename IndicesType::Nested>::type& 00528 indices() const { return m_indices; } 00529 00530 protected: 00531 00532 typename IndicesType::Nested m_indices; 00533 }; 00534 00537 template<typename Derived, typename PermutationDerived> 00538 inline const internal::permut_matrix_product_retval<PermutationDerived, Derived, OnTheRight> 00539 operator*(const MatrixBase<Derived>& matrix, 00540 const PermutationBase<PermutationDerived> &permutation) 00541 { 00542 return internal::permut_matrix_product_retval 00543 <PermutationDerived, Derived, OnTheRight> 00544 (permutation.derived(), matrix.derived()); 00545 } 00546 00549 template<typename Derived, typename PermutationDerived> 00550 inline const internal::permut_matrix_product_retval 00551 <PermutationDerived, Derived, OnTheLeft> 00552 operator*(const PermutationBase<PermutationDerived> &permutation, 00553 const MatrixBase<Derived>& matrix) 00554 { 00555 return internal::permut_matrix_product_retval 00556 <PermutationDerived, Derived, OnTheLeft> 00557 (permutation.derived(), matrix.derived()); 00558 } 00559 00560 namespace internal { 00561 00562 template<typename PermutationType, typename MatrixType, int Side, bool Transposed> 00563 struct traits<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> > 00564 { 00565 typedef typename MatrixType::PlainObject ReturnType; 00566 }; 00567 00568 template<typename PermutationType, typename MatrixType, int Side, bool Transposed> 00569 struct permut_matrix_product_retval 00570 : public ReturnByValue<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> > 00571 { 00572 typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned; 00573 typedef typename MatrixType::Index Index; 00574 00575 permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix) 00576 : m_permutation(perm), m_matrix(matrix) 00577 {} 00578 00579 inline Index rows() const { return m_matrix.rows(); } 00580 inline Index cols() const { return m_matrix.cols(); } 00581 00582 template<typename Dest> inline void evalTo(Dest& dst) const 00583 { 00584 const Index n = Side==OnTheLeft ? rows() : cols(); 00585 // FIXME we need an is_same for expression that is not sensitive to constness. For instance 00586 // is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true. 00587 if( is_same<MatrixTypeNestedCleaned,Dest>::value 00588 && blas_traits<MatrixTypeNestedCleaned>::HasUsableDirectAccess 00589 && blas_traits<Dest>::HasUsableDirectAccess 00590 && extract_data(dst) == extract_data(m_matrix)) 00591 { 00592 // apply the permutation inplace 00593 Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size()); 00594 mask.fill(false); 00595 Index r = 0; 00596 while(r < m_permutation.size()) 00597 { 00598 // search for the next seed 00599 while(r<m_permutation.size() && mask[r]) r++; 00600 if(r>=m_permutation.size()) 00601 break; 00602 // we got one, let's follow it until we are back to the seed 00603 Index k0 = r++; 00604 Index kPrev = k0; 00605 mask.coeffRef(k0) = true; 00606 for(Index k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k)) 00607 { 00608 Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k) 00609 .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime> 00610 (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev)); 00611 00612 mask.coeffRef(k) = true; 00613 kPrev = k; 00614 } 00615 } 00616 } 00617 else 00618 { 00619 for(int i = 0; i < n; ++i) 00620 { 00621 Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime> 00622 (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i) 00623 00624 = 00625 00626 Block<const MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime> 00627 (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i); 00628 } 00629 } 00630 } 00631 00632 protected: 00633 const PermutationType& m_permutation; 00634 typename MatrixType::Nested m_matrix; 00635 }; 00636 00637 /* Template partial specialization for transposed/inverse permutations */ 00638 00639 template<typename Derived> 00640 struct traits<Transpose<PermutationBase<Derived> > > 00641 : traits<Derived> 00642 {}; 00643 00644 } // end namespace internal 00645 00646 template<typename Derived> 00647 class Transpose<PermutationBase<Derived> > 00648 : public EigenBase<Transpose<PermutationBase<Derived> > > 00649 { 00650 typedef Derived PermutationType; 00651 typedef typename PermutationType::IndicesType IndicesType; 00652 typedef typename PermutationType::PlainPermutationType PlainPermutationType; 00653 public: 00654 00655 #ifndef EIGEN_PARSED_BY_DOXYGEN 00656 typedef internal::traits<PermutationType> Traits; 00657 typedef typename Derived::DenseMatrixType DenseMatrixType; 00658 enum { 00659 Flags = Traits::Flags, 00660 CoeffReadCost = Traits::CoeffReadCost, 00661 RowsAtCompileTime = Traits::RowsAtCompileTime, 00662 ColsAtCompileTime = Traits::ColsAtCompileTime, 00663 MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime, 00664 MaxColsAtCompileTime = Traits::MaxColsAtCompileTime 00665 }; 00666 typedef typename Traits::Scalar Scalar; 00667 #endif 00668 00669 Transpose(const PermutationType& p) : m_permutation(p) {} 00670 00671 inline int rows() const { return m_permutation.rows(); } 00672 inline int cols() const { return m_permutation.cols(); } 00673 00674 #ifndef EIGEN_PARSED_BY_DOXYGEN 00675 template<typename DenseDerived> 00676 void evalTo(MatrixBase<DenseDerived>& other) const 00677 { 00678 other.setZero(); 00679 for (int i=0; i<rows();++i) 00680 other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1); 00681 } 00682 #endif 00683 00685 PlainPermutationType eval() const { return *this; } 00686 00687 DenseMatrixType toDenseMatrix() const { return *this; } 00688 00691 template<typename OtherDerived> friend 00692 inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true> 00693 operator*(const MatrixBase<OtherDerived>& matrix, const Transpose& trPerm) 00694 { 00695 return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>(trPerm.m_permutation, matrix.derived()); 00696 } 00697 00700 template<typename OtherDerived> 00701 inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true> 00702 operator*(const MatrixBase<OtherDerived>& matrix) const 00703 { 00704 return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>(m_permutation, matrix.derived()); 00705 } 00706 00707 const PermutationType& nestedPermutation() const { return m_permutation; } 00708 00709 protected: 00710 const PermutationType& m_permutation; 00711 }; 00712 00713 template<typename Derived> 00714 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const 00715 { 00716 return derived(); 00717 } 00718 00719 } // end namespace Eigen 00720 00721 #endif // EIGEN_PERMUTATIONMATRIX_H