![]() |
Eigen
3.3.3
|
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-2015 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 namespace internal { 00017 00018 enum PermPermProduct_t {PermPermProduct}; 00019 00020 } // end namespace internal 00021 00045 template<typename Derived> 00046 class PermutationBase : public EigenBase<Derived> 00047 { 00048 typedef internal::traits<Derived> Traits; 00049 typedef EigenBase<Derived> Base; 00050 public: 00051 00052 #ifndef EIGEN_PARSED_BY_DOXYGEN 00053 typedef typename Traits::IndicesType IndicesType; 00054 enum { 00055 Flags = Traits::Flags, 00056 RowsAtCompileTime = Traits::RowsAtCompileTime, 00057 ColsAtCompileTime = Traits::ColsAtCompileTime, 00058 MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime, 00059 MaxColsAtCompileTime = Traits::MaxColsAtCompileTime 00060 }; 00061 typedef typename Traits::StorageIndex StorageIndex; 00062 typedef Matrix<StorageIndex,RowsAtCompileTime,ColsAtCompileTime,0,MaxRowsAtCompileTime,MaxColsAtCompileTime> 00063 DenseMatrixType; 00064 typedef PermutationMatrix<IndicesType::SizeAtCompileTime,IndicesType::MaxSizeAtCompileTime,StorageIndex> 00065 PlainPermutationType; 00066 typedef PlainPermutationType PlainObject; 00067 using Base::derived; 00068 typedef Inverse<Derived> InverseReturnType; 00069 typedef void Scalar; 00070 #endif 00071 00073 template<typename OtherDerived> 00074 Derived& operator=(const PermutationBase<OtherDerived>& other) 00075 { 00076 indices() = other.indices(); 00077 return derived(); 00078 } 00079 00081 template<typename OtherDerived> 00082 Derived& operator=(const TranspositionsBase<OtherDerived>& tr) 00083 { 00084 setIdentity(tr.size()); 00085 for(Index k=size()-1; k>=0; --k) 00086 applyTranspositionOnTheRight(k,tr.coeff(k)); 00087 return derived(); 00088 } 00089 00090 #ifndef EIGEN_PARSED_BY_DOXYGEN 00091 00094 Derived& operator=(const PermutationBase& other) 00095 { 00096 indices() = other.indices(); 00097 return derived(); 00098 } 00099 #endif 00100 00102 inline Index rows() const { return Index(indices().size()); } 00103 00105 inline Index cols() const { return Index(indices().size()); } 00106 00108 inline Index size() const { return Index(indices().size()); } 00109 00110 #ifndef EIGEN_PARSED_BY_DOXYGEN 00111 template<typename DenseDerived> 00112 void evalTo(MatrixBase<DenseDerived>& other) const 00113 { 00114 other.setZero(); 00115 for (Index i=0; i<rows(); ++i) 00116 other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1); 00117 } 00118 #endif 00119 00124 DenseMatrixType toDenseMatrix() const 00125 { 00126 return derived(); 00127 } 00128 00130 const IndicesType& indices() const { return derived().indices(); } 00132 IndicesType& indices() { return derived().indices(); } 00133 00136 inline void resize(Index newSize) 00137 { 00138 indices().resize(newSize); 00139 } 00140 00142 void setIdentity() 00143 { 00144 StorageIndex n = StorageIndex(size()); 00145 for(StorageIndex i = 0; i < n; ++i) 00146 indices().coeffRef(i) = i; 00147 } 00148 00151 void setIdentity(Index newSize) 00152 { 00153 resize(newSize); 00154 setIdentity(); 00155 } 00156 00166 Derived& applyTranspositionOnTheLeft(Index i, Index j) 00167 { 00168 eigen_assert(i>=0 && j>=0 && i<size() && j<size()); 00169 for(Index k = 0; k < size(); ++k) 00170 { 00171 if(indices().coeff(k) == i) indices().coeffRef(k) = StorageIndex(j); 00172 else if(indices().coeff(k) == j) indices().coeffRef(k) = StorageIndex(i); 00173 } 00174 return derived(); 00175 } 00176 00185 Derived& applyTranspositionOnTheRight(Index i, Index j) 00186 { 00187 eigen_assert(i>=0 && j>=0 && i<size() && j<size()); 00188 std::swap(indices().coeffRef(i), indices().coeffRef(j)); 00189 return derived(); 00190 } 00191 00196 inline InverseReturnType inverse() const 00197 { return InverseReturnType(derived()); } 00202 inline InverseReturnType transpose() const 00203 { return InverseReturnType(derived()); } 00204 00205 /**** multiplication helpers to hopefully get RVO ****/ 00206 00207 00208 #ifndef EIGEN_PARSED_BY_DOXYGEN 00209 protected: 00210 template<typename OtherDerived> 00211 void assignTranspose(const PermutationBase<OtherDerived>& other) 00212 { 00213 for (Index i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i; 00214 } 00215 template<typename Lhs,typename Rhs> 00216 void assignProduct(const Lhs& lhs, const Rhs& rhs) 00217 { 00218 eigen_assert(lhs.cols() == rhs.rows()); 00219 for (Index i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i)); 00220 } 00221 #endif 00222 00223 public: 00224 00229 template<typename Other> 00230 inline PlainPermutationType operator*(const PermutationBase<Other>& other) const 00231 { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); } 00232 00237 template<typename Other> 00238 inline PlainPermutationType operator*(const InverseImpl<Other,PermutationStorage>& other) const 00239 { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); } 00240 00245 template<typename Other> friend 00246 inline PlainPermutationType operator*(const InverseImpl<Other, PermutationStorage>& other, const PermutationBase& perm) 00247 { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); } 00248 00253 Index determinant() const 00254 { 00255 Index res = 1; 00256 Index n = size(); 00257 Matrix<bool,RowsAtCompileTime,1,0,MaxRowsAtCompileTime> mask(n); 00258 mask.fill(false); 00259 Index r = 0; 00260 while(r < n) 00261 { 00262 // search for the next seed 00263 while(r<n && mask[r]) r++; 00264 if(r>=n) 00265 break; 00266 // we got one, let's follow it until we are back to the seed 00267 Index k0 = r++; 00268 mask.coeffRef(k0) = true; 00269 for(Index k=indices().coeff(k0); k!=k0; k=indices().coeff(k)) 00270 { 00271 mask.coeffRef(k) = true; 00272 res = -res; 00273 } 00274 } 00275 return res; 00276 } 00277 00278 protected: 00279 00280 }; 00281 00282 namespace internal { 00283 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex> 00284 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> > 00285 : traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> > 00286 { 00287 typedef PermutationStorage StorageKind; 00288 typedef Matrix<_StorageIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType; 00289 typedef _StorageIndex StorageIndex; 00290 typedef void Scalar; 00291 }; 00292 } 00293 00307 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex> 00308 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> > 00309 { 00310 typedef PermutationBase<PermutationMatrix> Base; 00311 typedef internal::traits<PermutationMatrix> Traits; 00312 public: 00313 00314 typedef const PermutationMatrix& Nested; 00315 00316 #ifndef EIGEN_PARSED_BY_DOXYGEN 00317 typedef typename Traits::IndicesType IndicesType; 00318 typedef typename Traits::StorageIndex StorageIndex; 00319 #endif 00320 00321 inline PermutationMatrix() 00322 {} 00323 00326 explicit inline PermutationMatrix(Index size) : m_indices(size) 00327 { 00328 eigen_internal_assert(size <= NumTraits<StorageIndex>::highest()); 00329 } 00330 00332 template<typename OtherDerived> 00333 inline PermutationMatrix(const PermutationBase<OtherDerived>& other) 00334 : m_indices(other.indices()) {} 00335 00336 #ifndef EIGEN_PARSED_BY_DOXYGEN 00337 00339 inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {} 00340 #endif 00341 00349 template<typename Other> 00350 explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices) 00351 {} 00352 00354 template<typename Other> 00355 explicit PermutationMatrix(const TranspositionsBase<Other>& tr) 00356 : m_indices(tr.size()) 00357 { 00358 *this = tr; 00359 } 00360 00362 template<typename Other> 00363 PermutationMatrix& operator=(const PermutationBase<Other>& other) 00364 { 00365 m_indices = other.indices(); 00366 return *this; 00367 } 00368 00370 template<typename Other> 00371 PermutationMatrix& operator=(const TranspositionsBase<Other>& tr) 00372 { 00373 return Base::operator=(tr.derived()); 00374 } 00375 00376 #ifndef EIGEN_PARSED_BY_DOXYGEN 00377 00380 PermutationMatrix& operator=(const PermutationMatrix& other) 00381 { 00382 m_indices = other.m_indices; 00383 return *this; 00384 } 00385 #endif 00386 00388 const IndicesType& indices() const { return m_indices; } 00390 IndicesType& indices() { return m_indices; } 00391 00392 00393 /**** multiplication helpers to hopefully get RVO ****/ 00394 00395 #ifndef EIGEN_PARSED_BY_DOXYGEN 00396 template<typename Other> 00397 PermutationMatrix(const InverseImpl<Other,PermutationStorage>& other) 00398 : m_indices(other.derived().nestedExpression().size()) 00399 { 00400 eigen_internal_assert(m_indices.size() <= NumTraits<StorageIndex>::highest()); 00401 StorageIndex end = StorageIndex(m_indices.size()); 00402 for (StorageIndex i=0; i<end;++i) 00403 m_indices.coeffRef(other.derived().nestedExpression().indices().coeff(i)) = i; 00404 } 00405 template<typename Lhs,typename Rhs> 00406 PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs) 00407 : m_indices(lhs.indices().size()) 00408 { 00409 Base::assignProduct(lhs,rhs); 00410 } 00411 #endif 00412 00413 protected: 00414 00415 IndicesType m_indices; 00416 }; 00417 00418 00419 namespace internal { 00420 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess> 00421 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> > 00422 : traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> > 00423 { 00424 typedef PermutationStorage StorageKind; 00425 typedef Map<const Matrix<_StorageIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType; 00426 typedef _StorageIndex StorageIndex; 00427 typedef void Scalar; 00428 }; 00429 } 00430 00431 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess> 00432 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> 00433 : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> > 00434 { 00435 typedef PermutationBase<Map> Base; 00436 typedef internal::traits<Map> Traits; 00437 public: 00438 00439 #ifndef EIGEN_PARSED_BY_DOXYGEN 00440 typedef typename Traits::IndicesType IndicesType; 00441 typedef typename IndicesType::Scalar StorageIndex; 00442 #endif 00443 00444 inline Map(const StorageIndex* indicesPtr) 00445 : m_indices(indicesPtr) 00446 {} 00447 00448 inline Map(const StorageIndex* indicesPtr, Index size) 00449 : m_indices(indicesPtr,size) 00450 {} 00451 00453 template<typename Other> 00454 Map& operator=(const PermutationBase<Other>& other) 00455 { return Base::operator=(other.derived()); } 00456 00458 template<typename Other> 00459 Map& operator=(const TranspositionsBase<Other>& tr) 00460 { return Base::operator=(tr.derived()); } 00461 00462 #ifndef EIGEN_PARSED_BY_DOXYGEN 00463 00466 Map& operator=(const Map& other) 00467 { 00468 m_indices = other.m_indices; 00469 return *this; 00470 } 00471 #endif 00472 00474 const IndicesType& indices() const { return m_indices; } 00476 IndicesType& indices() { return m_indices; } 00477 00478 protected: 00479 00480 IndicesType m_indices; 00481 }; 00482 00483 template<typename _IndicesType> class TranspositionsWrapper; 00484 namespace internal { 00485 template<typename _IndicesType> 00486 struct traits<PermutationWrapper<_IndicesType> > 00487 { 00488 typedef PermutationStorage StorageKind; 00489 typedef void Scalar; 00490 typedef typename _IndicesType::Scalar StorageIndex; 00491 typedef _IndicesType IndicesType; 00492 enum { 00493 RowsAtCompileTime = _IndicesType::SizeAtCompileTime, 00494 ColsAtCompileTime = _IndicesType::SizeAtCompileTime, 00495 MaxRowsAtCompileTime = IndicesType::MaxSizeAtCompileTime, 00496 MaxColsAtCompileTime = IndicesType::MaxSizeAtCompileTime, 00497 Flags = 0 00498 }; 00499 }; 00500 } 00501 00513 template<typename _IndicesType> 00514 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> > 00515 { 00516 typedef PermutationBase<PermutationWrapper> Base; 00517 typedef internal::traits<PermutationWrapper> Traits; 00518 public: 00519 00520 #ifndef EIGEN_PARSED_BY_DOXYGEN 00521 typedef typename Traits::IndicesType IndicesType; 00522 #endif 00523 00524 inline PermutationWrapper(const IndicesType& indices) 00525 : m_indices(indices) 00526 {} 00527 00529 const typename internal::remove_all<typename IndicesType::Nested>::type& 00530 indices() const { return m_indices; } 00531 00532 protected: 00533 00534 typename IndicesType::Nested m_indices; 00535 }; 00536 00537 00540 template<typename MatrixDerived, typename PermutationDerived> 00541 EIGEN_DEVICE_FUNC 00542 const Product<MatrixDerived, PermutationDerived, AliasFreeProduct> 00543 operator*(const MatrixBase<MatrixDerived> &matrix, 00544 const PermutationBase<PermutationDerived>& permutation) 00545 { 00546 return Product<MatrixDerived, PermutationDerived, AliasFreeProduct> 00547 (matrix.derived(), permutation.derived()); 00548 } 00549 00552 template<typename PermutationDerived, typename MatrixDerived> 00553 EIGEN_DEVICE_FUNC 00554 const Product<PermutationDerived, MatrixDerived, AliasFreeProduct> 00555 operator*(const PermutationBase<PermutationDerived> &permutation, 00556 const MatrixBase<MatrixDerived>& matrix) 00557 { 00558 return Product<PermutationDerived, MatrixDerived, AliasFreeProduct> 00559 (permutation.derived(), matrix.derived()); 00560 } 00561 00562 00563 template<typename PermutationType> 00564 class InverseImpl<PermutationType, PermutationStorage> 00565 : public EigenBase<Inverse<PermutationType> > 00566 { 00567 typedef typename PermutationType::PlainPermutationType PlainPermutationType; 00568 typedef internal::traits<PermutationType> PermTraits; 00569 protected: 00570 InverseImpl() {} 00571 public: 00572 typedef Inverse<PermutationType> InverseType; 00573 using EigenBase<Inverse<PermutationType> >::derived; 00574 00575 #ifndef EIGEN_PARSED_BY_DOXYGEN 00576 typedef typename PermutationType::DenseMatrixType DenseMatrixType; 00577 enum { 00578 RowsAtCompileTime = PermTraits::RowsAtCompileTime, 00579 ColsAtCompileTime = PermTraits::ColsAtCompileTime, 00580 MaxRowsAtCompileTime = PermTraits::MaxRowsAtCompileTime, 00581 MaxColsAtCompileTime = PermTraits::MaxColsAtCompileTime 00582 }; 00583 #endif 00584 00585 #ifndef EIGEN_PARSED_BY_DOXYGEN 00586 template<typename DenseDerived> 00587 void evalTo(MatrixBase<DenseDerived>& other) const 00588 { 00589 other.setZero(); 00590 for (Index i=0; i<derived().rows();++i) 00591 other.coeffRef(i, derived().nestedExpression().indices().coeff(i)) = typename DenseDerived::Scalar(1); 00592 } 00593 #endif 00594 00596 PlainPermutationType eval() const { return derived(); } 00597 00598 DenseMatrixType toDenseMatrix() const { return derived(); } 00599 00602 template<typename OtherDerived> friend 00603 const Product<OtherDerived, InverseType, AliasFreeProduct> 00604 operator*(const MatrixBase<OtherDerived>& matrix, const InverseType& trPerm) 00605 { 00606 return Product<OtherDerived, InverseType, AliasFreeProduct>(matrix.derived(), trPerm.derived()); 00607 } 00608 00611 template<typename OtherDerived> 00612 const Product<InverseType, OtherDerived, AliasFreeProduct> 00613 operator*(const MatrixBase<OtherDerived>& matrix) const 00614 { 00615 return Product<InverseType, OtherDerived, AliasFreeProduct>(derived(), matrix.derived()); 00616 } 00617 }; 00618 00619 template<typename Derived> 00620 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const 00621 { 00622 return derived(); 00623 } 00624 00625 namespace internal { 00626 00627 template<> struct AssignmentKind<DenseShape,PermutationShape> { typedef EigenBase2EigenBase Kind; }; 00628 00629 } // end namespace internal 00630 00631 } // end namespace Eigen 00632 00633 #endif // EIGEN_PERMUTATIONMATRIX_H