Eigen  3.3.7
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 namespace Eigen {
15 
16 namespace internal {
17 
18 enum PermPermProduct_t {PermPermProduct};
19 
20 } // end namespace internal
21 
45 template<typename Derived>
46 class PermutationBase : public EigenBase<Derived>
47 {
48  typedef internal::traits<Derived> Traits;
49  typedef EigenBase<Derived> Base;
50  public:
51 
52  #ifndef EIGEN_PARSED_BY_DOXYGEN
53  typedef typename Traits::IndicesType IndicesType;
54  enum {
55  Flags = Traits::Flags,
56  RowsAtCompileTime = Traits::RowsAtCompileTime,
57  ColsAtCompileTime = Traits::ColsAtCompileTime,
58  MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
59  MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
60  };
61  typedef typename Traits::StorageIndex StorageIndex;
63  DenseMatrixType;
65  PlainPermutationType;
66  typedef PlainPermutationType PlainObject;
67  using Base::derived;
68  typedef Inverse<Derived> InverseReturnType;
69  typedef void Scalar;
70  #endif
71 
73  template<typename OtherDerived>
75  {
76  indices() = other.indices();
77  return derived();
78  }
79 
81  template<typename OtherDerived>
82  Derived& operator=(const TranspositionsBase<OtherDerived>& tr)
83  {
84  setIdentity(tr.size());
85  for(Index k=size()-1; k>=0; --k)
86  applyTranspositionOnTheRight(k,tr.coeff(k));
87  return derived();
88  }
89 
90  #ifndef EIGEN_PARSED_BY_DOXYGEN
91 
94  Derived& operator=(const PermutationBase& other)
95  {
96  indices() = other.indices();
97  return derived();
98  }
99  #endif
100 
102  inline Index rows() const { return Index(indices().size()); }
103 
105  inline Index cols() const { return Index(indices().size()); }
106 
108  inline Index size() const { return Index(indices().size()); }
109 
110  #ifndef EIGEN_PARSED_BY_DOXYGEN
111  template<typename DenseDerived>
112  void evalTo(MatrixBase<DenseDerived>& other) const
113  {
114  other.setZero();
115  for (Index i=0; i<rows(); ++i)
116  other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
117  }
118  #endif
119 
124  DenseMatrixType toDenseMatrix() const
125  {
126  return derived();
127  }
128 
130  const IndicesType& indices() const { return derived().indices(); }
132  IndicesType& indices() { return derived().indices(); }
133 
136  inline void resize(Index newSize)
137  {
138  indices().resize(newSize);
139  }
140 
142  void setIdentity()
143  {
144  StorageIndex n = StorageIndex(size());
145  for(StorageIndex i = 0; i < n; ++i)
146  indices().coeffRef(i) = i;
147  }
148 
151  void setIdentity(Index newSize)
152  {
153  resize(newSize);
154  setIdentity();
155  }
156 
167  {
168  eigen_assert(i>=0 && j>=0 && i<size() && j<size());
169  for(Index k = 0; k < size(); ++k)
170  {
171  if(indices().coeff(k) == i) indices().coeffRef(k) = StorageIndex(j);
172  else if(indices().coeff(k) == j) indices().coeffRef(k) = StorageIndex(i);
173  }
174  return derived();
175  }
176 
186  {
187  eigen_assert(i>=0 && j>=0 && i<size() && j<size());
188  std::swap(indices().coeffRef(i), indices().coeffRef(j));
189  return derived();
190  }
191 
196  inline InverseReturnType inverse() const
197  { return InverseReturnType(derived()); }
202  inline InverseReturnType transpose() const
203  { return InverseReturnType(derived()); }
204 
205  /**** multiplication helpers to hopefully get RVO ****/
206 
207 
208 #ifndef EIGEN_PARSED_BY_DOXYGEN
209  protected:
210  template<typename OtherDerived>
211  void assignTranspose(const PermutationBase<OtherDerived>& other)
212  {
213  for (Index i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i;
214  }
215  template<typename Lhs,typename Rhs>
216  void assignProduct(const Lhs& lhs, const Rhs& rhs)
217  {
218  eigen_assert(lhs.cols() == rhs.rows());
219  for (Index i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
220  }
221 #endif
222 
223  public:
224 
229  template<typename Other>
230  inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
231  { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
232 
237  template<typename Other>
238  inline PlainPermutationType operator*(const InverseImpl<Other,PermutationStorage>& other) const
239  { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
240 
245  template<typename Other> friend
246  inline PlainPermutationType operator*(const InverseImpl<Other, PermutationStorage>& other, const PermutationBase& perm)
247  { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
248 
254  {
255  Index res = 1;
256  Index n = size();
258  mask.fill(false);
259  Index r = 0;
260  while(r < n)
261  {
262  // search for the next seed
263  while(r<n && mask[r]) r++;
264  if(r>=n)
265  break;
266  // we got one, let's follow it until we are back to the seed
267  Index k0 = r++;
268  mask.coeffRef(k0) = true;
269  for(Index k=indices().coeff(k0); k!=k0; k=indices().coeff(k))
270  {
271  mask.coeffRef(k) = true;
272  res = -res;
273  }
274  }
275  return res;
276  }
277 
278  protected:
279 
280 };
281 
282 namespace internal {
283 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex>
284 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> >
285  : traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
286 {
287  typedef PermutationStorage StorageKind;
288  typedef Matrix<_StorageIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
289  typedef _StorageIndex StorageIndex;
290  typedef void Scalar;
291 };
292 }
293 
307 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex>
308 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> >
309 {
311  typedef internal::traits<PermutationMatrix> Traits;
312  public:
313 
314  typedef const PermutationMatrix& Nested;
315 
316  #ifndef EIGEN_PARSED_BY_DOXYGEN
317  typedef typename Traits::IndicesType IndicesType;
318  typedef typename Traits::StorageIndex StorageIndex;
319  #endif
320 
321  inline PermutationMatrix()
322  {}
323 
326  explicit inline PermutationMatrix(Index size) : m_indices(size)
327  {
328  eigen_internal_assert(size <= NumTraits<StorageIndex>::highest());
329  }
330 
332  template<typename OtherDerived>
334  : m_indices(other.indices()) {}
335 
336  #ifndef EIGEN_PARSED_BY_DOXYGEN
337 
339  inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
340  #endif
341 
349  template<typename Other>
350  explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
351  {}
352 
354  template<typename Other>
355  explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
356  : m_indices(tr.size())
357  {
358  *this = tr;
359  }
360 
362  template<typename Other>
364  {
365  m_indices = other.indices();
366  return *this;
367  }
368 
370  template<typename Other>
371  PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
372  {
373  return Base::operator=(tr.derived());
374  }
375 
376  #ifndef EIGEN_PARSED_BY_DOXYGEN
377 
381  {
382  m_indices = other.m_indices;
383  return *this;
384  }
385  #endif
386 
388  const IndicesType& indices() const { return m_indices; }
390  IndicesType& indices() { return m_indices; }
391 
392 
393  /**** multiplication helpers to hopefully get RVO ****/
394 
395 #ifndef EIGEN_PARSED_BY_DOXYGEN
396  template<typename Other>
397  PermutationMatrix(const InverseImpl<Other,PermutationStorage>& other)
398  : m_indices(other.derived().nestedExpression().size())
399  {
400  eigen_internal_assert(m_indices.size() <= NumTraits<StorageIndex>::highest());
401  StorageIndex end = StorageIndex(m_indices.size());
402  for (StorageIndex i=0; i<end;++i)
403  m_indices.coeffRef(other.derived().nestedExpression().indices().coeff(i)) = i;
404  }
405  template<typename Lhs,typename Rhs>
406  PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
407  : m_indices(lhs.indices().size())
408  {
409  Base::assignProduct(lhs,rhs);
410  }
411 #endif
412 
413  protected:
414 
415  IndicesType m_indices;
416 };
417 
418 
419 namespace internal {
420 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess>
421 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> >
422  : traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
423 {
424  typedef PermutationStorage StorageKind;
425  typedef Map<const Matrix<_StorageIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType;
426  typedef _StorageIndex StorageIndex;
427  typedef void Scalar;
428 };
429 }
430 
431 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess>
432 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess>
433  : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> >
434 {
435  typedef PermutationBase<Map> Base;
436  typedef internal::traits<Map> Traits;
437  public:
438 
439  #ifndef EIGEN_PARSED_BY_DOXYGEN
440  typedef typename Traits::IndicesType IndicesType;
441  typedef typename IndicesType::Scalar StorageIndex;
442  #endif
443 
444  inline Map(const StorageIndex* indicesPtr)
445  : m_indices(indicesPtr)
446  {}
447 
448  inline Map(const StorageIndex* indicesPtr, Index size)
449  : m_indices(indicesPtr,size)
450  {}
451 
453  template<typename Other>
454  Map& operator=(const PermutationBase<Other>& other)
455  { return Base::operator=(other.derived()); }
456 
458  template<typename Other>
459  Map& operator=(const TranspositionsBase<Other>& tr)
460  { return Base::operator=(tr.derived()); }
461 
462  #ifndef EIGEN_PARSED_BY_DOXYGEN
463 
466  Map& operator=(const Map& other)
467  {
468  m_indices = other.m_indices;
469  return *this;
470  }
471  #endif
472 
474  const IndicesType& indices() const { return m_indices; }
476  IndicesType& indices() { return m_indices; }
477 
478  protected:
479 
480  IndicesType m_indices;
481 };
482 
483 template<typename _IndicesType> class TranspositionsWrapper;
484 namespace internal {
485 template<typename _IndicesType>
486 struct traits<PermutationWrapper<_IndicesType> >
487 {
488  typedef PermutationStorage StorageKind;
489  typedef void Scalar;
490  typedef typename _IndicesType::Scalar StorageIndex;
491  typedef _IndicesType IndicesType;
492  enum {
493  RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
494  ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
495  MaxRowsAtCompileTime = IndicesType::MaxSizeAtCompileTime,
496  MaxColsAtCompileTime = IndicesType::MaxSizeAtCompileTime,
497  Flags = 0
498  };
499 };
500 }
501 
513 template<typename _IndicesType>
514 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
515 {
517  typedef internal::traits<PermutationWrapper> Traits;
518  public:
519 
520  #ifndef EIGEN_PARSED_BY_DOXYGEN
521  typedef typename Traits::IndicesType IndicesType;
522  #endif
523 
524  inline PermutationWrapper(const IndicesType& indices)
525  : m_indices(indices)
526  {}
527 
529  const typename internal::remove_all<typename IndicesType::Nested>::type&
530  indices() const { return m_indices; }
531 
532  protected:
533 
534  typename IndicesType::Nested m_indices;
535 };
536 
537 
540 template<typename MatrixDerived, typename PermutationDerived>
541 EIGEN_DEVICE_FUNC
542 const Product<MatrixDerived, PermutationDerived, AliasFreeProduct>
544  const PermutationBase<PermutationDerived>& permutation)
545 {
547  (matrix.derived(), permutation.derived());
548 }
549 
552 template<typename PermutationDerived, typename MatrixDerived>
553 EIGEN_DEVICE_FUNC
554 const Product<PermutationDerived, MatrixDerived, AliasFreeProduct>
556  const MatrixBase<MatrixDerived>& matrix)
557 {
559  (permutation.derived(), matrix.derived());
560 }
561 
562 
563 template<typename PermutationType>
564 class InverseImpl<PermutationType, PermutationStorage>
565  : public EigenBase<Inverse<PermutationType> >
566 {
567  typedef typename PermutationType::PlainPermutationType PlainPermutationType;
568  typedef internal::traits<PermutationType> PermTraits;
569  protected:
570  InverseImpl() {}
571  public:
572  typedef Inverse<PermutationType> InverseType;
573  using EigenBase<Inverse<PermutationType> >::derived;
574 
575  #ifndef EIGEN_PARSED_BY_DOXYGEN
576  typedef typename PermutationType::DenseMatrixType DenseMatrixType;
577  enum {
578  RowsAtCompileTime = PermTraits::RowsAtCompileTime,
579  ColsAtCompileTime = PermTraits::ColsAtCompileTime,
580  MaxRowsAtCompileTime = PermTraits::MaxRowsAtCompileTime,
581  MaxColsAtCompileTime = PermTraits::MaxColsAtCompileTime
582  };
583  #endif
584 
585  #ifndef EIGEN_PARSED_BY_DOXYGEN
586  template<typename DenseDerived>
587  void evalTo(MatrixBase<DenseDerived>& other) const
588  {
589  other.setZero();
590  for (Index i=0; i<derived().rows();++i)
591  other.coeffRef(i, derived().nestedExpression().indices().coeff(i)) = typename DenseDerived::Scalar(1);
592  }
593  #endif
594 
596  PlainPermutationType eval() const { return derived(); }
597 
598  DenseMatrixType toDenseMatrix() const { return derived(); }
599 
602  template<typename OtherDerived> friend
603  const Product<OtherDerived, InverseType, AliasFreeProduct>
604  operator*(const MatrixBase<OtherDerived>& matrix, const InverseType& trPerm)
605  {
606  return Product<OtherDerived, InverseType, AliasFreeProduct>(matrix.derived(), trPerm.derived());
607  }
608 
611  template<typename OtherDerived>
612  const Product<InverseType, OtherDerived, AliasFreeProduct>
613  operator*(const MatrixBase<OtherDerived>& matrix) const
614  {
615  return Product<InverseType, OtherDerived, AliasFreeProduct>(derived(), matrix.derived());
616  }
617 };
618 
619 template<typename Derived>
620 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
621 {
622  return derived();
623 }
624 
625 namespace internal {
626 
627 template<> struct AssignmentKind<DenseShape,PermutationShape> { typedef EigenBase2EigenBase Kind; };
628 
629 } // end namespace internal
630 
631 } // end namespace Eigen
632 
633 #endif // EIGEN_PERMUTATIONMATRIX_H
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:71
PlainPermutationType operator*(const PermutationBase< Other > &other) const
Definition: PermutationMatrix.h:230
PermutationMatrix(const MatrixBase< Other > &indices)
Definition: PermutationMatrix.h:350
friend PlainPermutationType operator*(const InverseImpl< Other, PermutationStorage > &other, const PermutationBase &perm)
Definition: PermutationMatrix.h:246
InverseReturnType inverse() const
Definition: PermutationMatrix.h:196
PlainPermutationType operator*(const InverseImpl< Other, PermutationStorage > &other) const
Definition: PermutationMatrix.h:238
DenseMatrixType toDenseMatrix() const
Definition: PermutationMatrix.h:124
Derived & applyTranspositionOnTheRight(Index i, Index j)
Definition: PermutationMatrix.h:185
Namespace containing all symbols from the Eigen library.
Definition: Core:306
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
Derived & derived()
Definition: EigenBase.h:45
Derived & operator=(const PermutationBase< OtherDerived > &other)
Definition: PermutationMatrix.h:74
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:38
PermutationMatrix & operator=(const PermutationBase< Other > &other)
Definition: PermutationMatrix.h:363
Map(PointerArgType dataPtr, const StrideType &stride=StrideType())
Definition: Map.h:129
Base class for permutations.
Definition: PermutationMatrix.h:46
Definition: EigenBase.h:29
Expression of the inverse of another expression.
Definition: Inverse.h:43
Permutation matrix.
Definition: PermutationMatrix.h:308
Index size() const
Definition: PermutationMatrix.h:108
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:183
PermutationMatrix & operator=(const TranspositionsBase< Other > &tr)
Definition: PermutationMatrix.h:371
IndicesType & indices()
Definition: PermutationMatrix.h:132
Index cols() const
Definition: PermutationMatrix.h:105
const IndicesType & indices() const
Definition: PermutationMatrix.h:388
void fill(const Scalar &value)
Definition: CwiseNullaryOp.h:315
Derived & operator=(const TranspositionsBase< OtherDerived > &tr)
Definition: PermutationMatrix.h:82
Derived & derived()
Definition: EigenBase.h:45
const Product< MatrixDerived, PermutationDerived, AliasFreeProduct > operator*(const MatrixBase< MatrixDerived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:543
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Class to view a vector of integers as a permutation matrix.
Definition: PermutationMatrix.h:514
void setIdentity()
Definition: PermutationMatrix.h:142
const IndicesType & indices() const
Definition: PermutationMatrix.h:130
PermutationMatrix(const TranspositionsBase< Other > &tr)
Definition: PermutationMatrix.h:355
Definition: Eigen_Colamd.h:50
Derived & setZero()
Definition: CwiseNullaryOp.h:499
Derived & applyTranspositionOnTheLeft(Index i, Index j)
Definition: PermutationMatrix.h:166
PermutationMatrix(Index size)
Definition: PermutationMatrix.h:326
InverseReturnType transpose() const
Definition: PermutationMatrix.h:202
IndicesType & indices()
Definition: PermutationMatrix.h:390
const internal::remove_all< typename IndicesType::Nested >::type & indices() const
Definition: PermutationMatrix.h:530
void setIdentity(Index newSize)
Definition: PermutationMatrix.h:151
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
void resize(Index newSize)
Definition: PermutationMatrix.h:136
Scalar & coeffRef(Index row, Index col)
Definition: DenseCoeffsBase.h:340
Index determinant() const
Definition: PermutationMatrix.h:253
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
PermutationMatrix(const PermutationBase< OtherDerived > &other)
Definition: PermutationMatrix.h:333
Index rows() const
Definition: PermutationMatrix.h:102