Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
PermutationMatrix.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
5 // Copyright (C) 2009-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_PERMUTATIONMATRIX_H
12 #define EIGEN_PERMUTATIONMATRIX_H
13 
14 #include "./InternalHeaderCheck.h"
15 
16 namespace Eigen {
17 
18 namespace internal {
19 
20 enum PermPermProduct_t {PermPermProduct};
21 
22 } // end namespace internal
23 
47 template<typename Derived>
48 class PermutationBase : public EigenBase<Derived>
49 {
50  typedef internal::traits<Derived> Traits;
51  typedef EigenBase<Derived> Base;
52  public:
53 
54  #ifndef EIGEN_PARSED_BY_DOXYGEN
55  typedef typename Traits::IndicesType IndicesType;
56  enum {
57  Flags = Traits::Flags,
58  RowsAtCompileTime = Traits::RowsAtCompileTime,
59  ColsAtCompileTime = Traits::ColsAtCompileTime,
60  MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
61  MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
62  };
63  typedef typename Traits::StorageIndex StorageIndex;
65  DenseMatrixType;
67  PlainPermutationType;
68  typedef PlainPermutationType PlainObject;
69  using Base::derived;
70  typedef Inverse<Derived> InverseReturnType;
71  typedef void Scalar;
72  #endif
73 
75  template<typename OtherDerived>
77  {
78  indices() = other.indices();
79  return derived();
80  }
81 
83  template<typename OtherDerived>
84  Derived& operator=(const TranspositionsBase<OtherDerived>& tr)
85  {
86  setIdentity(tr.size());
87  for(Index k=size()-1; k>=0; --k)
88  applyTranspositionOnTheRight(k,tr.coeff(k));
89  return derived();
90  }
91 
93  inline EIGEN_DEVICE_FUNC Index rows() const { return Index(indices().size()); }
94 
96  inline EIGEN_DEVICE_FUNC Index cols() const { return Index(indices().size()); }
97 
99  inline EIGEN_DEVICE_FUNC Index size() const { return Index(indices().size()); }
100 
101  #ifndef EIGEN_PARSED_BY_DOXYGEN
102  template<typename DenseDerived>
103  void evalTo(MatrixBase<DenseDerived>& other) const
104  {
105  other.setZero();
106  for (Index i=0; i<rows(); ++i)
107  other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
108  }
109  #endif
110 
115  DenseMatrixType toDenseMatrix() const
116  {
117  return derived();
118  }
119 
121  const IndicesType& indices() const { return derived().indices(); }
123  IndicesType& indices() { return derived().indices(); }
124 
127  inline void resize(Index newSize)
128  {
129  indices().resize(newSize);
130  }
131 
133  void setIdentity()
134  {
135  StorageIndex n = StorageIndex(size());
136  for(StorageIndex i = 0; i < n; ++i)
137  indices().coeffRef(i) = i;
138  }
139 
142  void setIdentity(Index newSize)
143  {
144  resize(newSize);
145  setIdentity();
146  }
147 
158  {
159  eigen_assert(i>=0 && j>=0 && i<size() && j<size());
160  for(Index k = 0; k < size(); ++k)
161  {
162  if(indices().coeff(k) == i) indices().coeffRef(k) = StorageIndex(j);
163  else if(indices().coeff(k) == j) indices().coeffRef(k) = StorageIndex(i);
164  }
165  return derived();
166  }
167 
177  {
178  eigen_assert(i>=0 && j>=0 && i<size() && j<size());
179  std::swap(indices().coeffRef(i), indices().coeffRef(j));
180  return derived();
181  }
182 
187  inline InverseReturnType inverse() const
188  { return InverseReturnType(derived()); }
193  inline InverseReturnType transpose() const
194  { return InverseReturnType(derived()); }
195 
196  /**** multiplication helpers to hopefully get RVO ****/
197 
198 
199 #ifndef EIGEN_PARSED_BY_DOXYGEN
200  protected:
201  template<typename OtherDerived>
202  void assignTranspose(const PermutationBase<OtherDerived>& other)
203  {
204  for (Index i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i;
205  }
206  template<typename Lhs,typename Rhs>
207  void assignProduct(const Lhs& lhs, const Rhs& rhs)
208  {
209  eigen_assert(lhs.cols() == rhs.rows());
210  for (Index i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
211  }
212 #endif
213 
214  public:
215 
220  template<typename Other>
221  inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
222  { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
223 
228  template<typename Other>
229  inline PlainPermutationType operator*(const InverseImpl<Other,PermutationStorage>& other) const
230  { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
231 
236  template<typename Other> friend
237  inline PlainPermutationType operator*(const InverseImpl<Other, PermutationStorage>& other, const PermutationBase& perm)
238  { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
239 
245  {
246  Index res = 1;
247  Index n = size();
249  mask.fill(false);
250  Index r = 0;
251  while(r < n)
252  {
253  // search for the next seed
254  while(r<n && mask[r]) r++;
255  if(r>=n)
256  break;
257  // we got one, let's follow it until we are back to the seed
258  Index k0 = r++;
259  mask.coeffRef(k0) = true;
260  for(Index k=indices().coeff(k0); k!=k0; k=indices().coeff(k))
261  {
262  mask.coeffRef(k) = true;
263  res = -res;
264  }
265  }
266  return res;
267  }
268 
269  protected:
270 
271 };
272 
273 namespace internal {
274 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename StorageIndex_>
275 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_> >
276  : traits<Matrix<StorageIndex_,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
277 {
278  typedef PermutationStorage StorageKind;
279  typedef Matrix<StorageIndex_, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
280  typedef StorageIndex_ StorageIndex;
281  typedef void Scalar;
282 };
283 }
284 
298 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename StorageIndex_>
299 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_> >
300 {
302  typedef internal::traits<PermutationMatrix> Traits;
303  public:
304 
305  typedef const PermutationMatrix& Nested;
306 
307  #ifndef EIGEN_PARSED_BY_DOXYGEN
308  typedef typename Traits::IndicesType IndicesType;
309  typedef typename Traits::StorageIndex StorageIndex;
310  #endif
311 
312  inline PermutationMatrix()
313  {}
314 
317  explicit inline PermutationMatrix(Index size) : m_indices(size)
318  {
319  eigen_internal_assert(size <= NumTraits<StorageIndex>::highest());
320  }
321 
323  template<typename OtherDerived>
325  : m_indices(other.indices()) {}
326 
334  template<typename Other>
335  explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
336  {}
337 
339  template<typename Other>
340  explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
341  : m_indices(tr.size())
342  {
343  *this = tr;
344  }
345 
347  template<typename Other>
349  {
350  m_indices = other.indices();
351  return *this;
352  }
353 
355  template<typename Other>
356  PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
357  {
358  return Base::operator=(tr.derived());
359  }
360 
362  const IndicesType& indices() const { return m_indices; }
364  IndicesType& indices() { return m_indices; }
365 
366 
367  /**** multiplication helpers to hopefully get RVO ****/
368 
369 #ifndef EIGEN_PARSED_BY_DOXYGEN
370  template<typename Other>
371  PermutationMatrix(const InverseImpl<Other,PermutationStorage>& other)
372  : m_indices(other.derived().nestedExpression().size())
373  {
374  eigen_internal_assert(m_indices.size() <= NumTraits<StorageIndex>::highest());
375  StorageIndex end = StorageIndex(m_indices.size());
376  for (StorageIndex i=0; i<end;++i)
377  m_indices.coeffRef(other.derived().nestedExpression().indices().coeff(i)) = i;
378  }
379  template<typename Lhs,typename Rhs>
380  PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
381  : m_indices(lhs.indices().size())
382  {
383  Base::assignProduct(lhs,rhs);
384  }
385 #endif
386 
387  protected:
388 
389  IndicesType m_indices;
390 };
391 
392 
393 namespace internal {
394 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename StorageIndex_, int PacketAccess_>
395 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_>,PacketAccess_> >
396  : traits<Matrix<StorageIndex_,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
397 {
398  typedef PermutationStorage StorageKind;
399  typedef Map<const Matrix<StorageIndex_, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, PacketAccess_> IndicesType;
400  typedef StorageIndex_ StorageIndex;
401  typedef void Scalar;
402 };
403 }
404 
405 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename StorageIndex_, int PacketAccess_>
406 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_>,PacketAccess_>
407  : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_>,PacketAccess_> >
408 {
409  typedef PermutationBase<Map> Base;
410  typedef internal::traits<Map> Traits;
411  public:
412 
413  #ifndef EIGEN_PARSED_BY_DOXYGEN
414  typedef typename Traits::IndicesType IndicesType;
415  typedef typename IndicesType::Scalar StorageIndex;
416  #endif
417 
418  inline Map(const StorageIndex* indicesPtr)
419  : m_indices(indicesPtr)
420  {}
421 
422  inline Map(const StorageIndex* indicesPtr, Index size)
423  : m_indices(indicesPtr,size)
424  {}
425 
427  template<typename Other>
428  Map& operator=(const PermutationBase<Other>& other)
429  { return Base::operator=(other.derived()); }
430 
432  template<typename Other>
433  Map& operator=(const TranspositionsBase<Other>& tr)
434  { return Base::operator=(tr.derived()); }
435 
436  #ifndef EIGEN_PARSED_BY_DOXYGEN
440  Map& operator=(const Map& other)
441  {
442  m_indices = other.m_indices;
443  return *this;
444  }
445  #endif
446 
448  const IndicesType& indices() const { return m_indices; }
450  IndicesType& indices() { return m_indices; }
451 
452  protected:
453 
454  IndicesType m_indices;
455 };
456 
457 template<typename IndicesType_> class TranspositionsWrapper;
458 namespace internal {
459 template<typename IndicesType_>
460 struct traits<PermutationWrapper<IndicesType_> >
461 {
462  typedef PermutationStorage StorageKind;
463  typedef void Scalar;
464  typedef typename IndicesType_::Scalar StorageIndex;
465  typedef IndicesType_ IndicesType;
466  enum {
467  RowsAtCompileTime = IndicesType_::SizeAtCompileTime,
468  ColsAtCompileTime = IndicesType_::SizeAtCompileTime,
469  MaxRowsAtCompileTime = IndicesType::MaxSizeAtCompileTime,
470  MaxColsAtCompileTime = IndicesType::MaxSizeAtCompileTime,
471  Flags = 0
472  };
473 };
474 }
475 
487 template<typename IndicesType_>
488 class PermutationWrapper : public PermutationBase<PermutationWrapper<IndicesType_> >
489 {
491  typedef internal::traits<PermutationWrapper> Traits;
492  public:
493 
494  #ifndef EIGEN_PARSED_BY_DOXYGEN
495  typedef typename Traits::IndicesType IndicesType;
496  #endif
497 
498  inline PermutationWrapper(const IndicesType& indices)
499  : m_indices(indices)
500  {}
501 
503  const internal::remove_all_t<typename IndicesType::Nested>&
504  indices() const { return m_indices; }
505 
506  protected:
507 
508  typename IndicesType::Nested m_indices;
509 };
510 
511 
514 template<typename MatrixDerived, typename PermutationDerived>
515 EIGEN_DEVICE_FUNC
516 const Product<MatrixDerived, PermutationDerived, AliasFreeProduct>
518  const PermutationBase<PermutationDerived>& permutation)
519 {
521  (matrix.derived(), permutation.derived());
522 }
523 
526 template<typename PermutationDerived, typename MatrixDerived>
527 EIGEN_DEVICE_FUNC
528 const Product<PermutationDerived, MatrixDerived, AliasFreeProduct>
530  const MatrixBase<MatrixDerived>& matrix)
531 {
533  (permutation.derived(), matrix.derived());
534 }
535 
536 
537 template<typename PermutationType>
538 class InverseImpl<PermutationType, PermutationStorage>
539  : public EigenBase<Inverse<PermutationType> >
540 {
541  typedef typename PermutationType::PlainPermutationType PlainPermutationType;
542  typedef internal::traits<PermutationType> PermTraits;
543  protected:
544  InverseImpl() {}
545  public:
546  typedef Inverse<PermutationType> InverseType;
547  using EigenBase<Inverse<PermutationType> >::derived;
548 
549  #ifndef EIGEN_PARSED_BY_DOXYGEN
550  typedef typename PermutationType::DenseMatrixType DenseMatrixType;
551  enum {
552  RowsAtCompileTime = PermTraits::RowsAtCompileTime,
553  ColsAtCompileTime = PermTraits::ColsAtCompileTime,
554  MaxRowsAtCompileTime = PermTraits::MaxRowsAtCompileTime,
555  MaxColsAtCompileTime = PermTraits::MaxColsAtCompileTime
556  };
557  #endif
558 
559  #ifndef EIGEN_PARSED_BY_DOXYGEN
560  template<typename DenseDerived>
561  void evalTo(MatrixBase<DenseDerived>& other) const
562  {
563  other.setZero();
564  for (Index i=0; i<derived().rows();++i)
565  other.coeffRef(i, derived().nestedExpression().indices().coeff(i)) = typename DenseDerived::Scalar(1);
566  }
567  #endif
568 
570  PlainPermutationType eval() const { return derived(); }
571 
572  DenseMatrixType toDenseMatrix() const { return derived(); }
573 
576  template<typename OtherDerived> friend
577  const Product<OtherDerived, InverseType, AliasFreeProduct>
578  operator*(const MatrixBase<OtherDerived>& matrix, const InverseType& trPerm)
579  {
580  return Product<OtherDerived, InverseType, AliasFreeProduct>(matrix.derived(), trPerm.derived());
581  }
582 
585  template<typename OtherDerived>
586  const Product<InverseType, OtherDerived, AliasFreeProduct>
587  operator*(const MatrixBase<OtherDerived>& matrix) const
588  {
589  return Product<InverseType, OtherDerived, AliasFreeProduct>(derived(), matrix.derived());
590  }
591 };
592 
593 template<typename Derived>
594 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
595 {
596  return derived();
597 }
598 
599 namespace internal {
600 
601 template<> struct AssignmentKind<DenseShape,PermutationShape> { typedef EigenBase2EigenBase Kind; };
602 
603 } // end namespace internal
604 
605 } // end namespace Eigen
606 
607 #endif // EIGEN_PERMUTATIONMATRIX_H
void fill(const Scalar &value)
Definition: CwiseNullaryOp.h:337
Derived & setZero()
Definition: CwiseNullaryOp.h:548
Derived & derived()
Definition: EigenBase.h:48
Scalar & coeffRef(Index row, Index col)
Definition: DenseCoeffsBase.h:344
Expression of the inverse of another expression.
Definition: Inverse.h:46
Map(PointerArgType dataPtr, const StrideType &stride=StrideType())
Definition: Map.h:131
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:182
Base class for permutations.
Definition: PermutationMatrix.h:49
InverseReturnType transpose() const
Definition: PermutationMatrix.h:193
void resize(Index newSize)
Definition: PermutationMatrix.h:127
IndicesType & indices()
Definition: PermutationMatrix.h:123
Index determinant() const
Definition: PermutationMatrix.h:244
Index size() const
Definition: PermutationMatrix.h:99
Index cols() const
Definition: PermutationMatrix.h:96
friend PlainPermutationType operator*(const InverseImpl< Other, PermutationStorage > &other, const PermutationBase &perm)
Definition: PermutationMatrix.h:237
Derived & applyTranspositionOnTheLeft(Index i, Index j)
Definition: PermutationMatrix.h:157
Derived & applyTranspositionOnTheRight(Index i, Index j)
Definition: PermutationMatrix.h:176
void setIdentity()
Definition: PermutationMatrix.h:133
void setIdentity(Index newSize)
Definition: PermutationMatrix.h:142
PlainPermutationType operator*(const InverseImpl< Other, PermutationStorage > &other) const
Definition: PermutationMatrix.h:229
Derived & operator=(const PermutationBase< OtherDerived > &other)
Definition: PermutationMatrix.h:76
Derived & operator=(const TranspositionsBase< OtherDerived > &tr)
Definition: PermutationMatrix.h:84
Index rows() const
Definition: PermutationMatrix.h:93
InverseReturnType inverse() const
Definition: PermutationMatrix.h:187
DenseMatrixType toDenseMatrix() const
Definition: PermutationMatrix.h:115
const IndicesType & indices() const
Definition: PermutationMatrix.h:121
PlainPermutationType operator*(const PermutationBase< Other > &other) const
Definition: PermutationMatrix.h:221
Permutation matrix.
Definition: PermutationMatrix.h:300
PermutationMatrix & operator=(const TranspositionsBase< Other > &tr)
Definition: PermutationMatrix.h:356
PermutationMatrix(const PermutationBase< OtherDerived > &other)
Definition: PermutationMatrix.h:324
PermutationMatrix(const TranspositionsBase< Other > &tr)
Definition: PermutationMatrix.h:340
const IndicesType & indices() const
Definition: PermutationMatrix.h:362
IndicesType & indices()
Definition: PermutationMatrix.h:364
PermutationMatrix(const MatrixBase< Other > &indices)
Definition: PermutationMatrix.h:335
PermutationMatrix(Index size)
Definition: PermutationMatrix.h:317
PermutationMatrix & operator=(const PermutationBase< Other > &other)
Definition: PermutationMatrix.h:348
Class to view a vector of integers as a permutation matrix.
Definition: PermutationMatrix.h:489
const internal::remove_all_t< typename IndicesType::Nested > & indices() const
Definition: PermutationMatrix.h:504
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:187
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:77
static const lastp1_t end
Definition: IndexedViewHelper.h:183
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Product< MatrixDerived, PermutationDerived, AliasFreeProduct > operator*(const MatrixBase< MatrixDerived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:517
Definition: EigenBase.h:32
Derived & derived()
Definition: EigenBase.h:48
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231