Eigen  3.2.91
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-2011 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 // TODO: this does not seems to be needed at all:
17 // template<int RowCol,typename IndicesType,typename MatrixType, typename StorageKind> class PermutedImpl;
18 
43 namespace internal {
44 
45 enum PermPermProduct_t {PermPermProduct};
46 
47 } // end namespace internal
48 
49 template<typename Derived>
50 class PermutationBase : public EigenBase<Derived>
51 {
52  typedef internal::traits<Derived> Traits;
53  typedef EigenBase<Derived> Base;
54  public:
55 
56  #ifndef EIGEN_PARSED_BY_DOXYGEN
57  typedef typename Traits::IndicesType IndicesType;
58  enum {
59  Flags = Traits::Flags,
60  RowsAtCompileTime = Traits::RowsAtCompileTime,
61  ColsAtCompileTime = Traits::ColsAtCompileTime,
62  MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
63  MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
64  };
65  typedef typename Traits::StorageIndex StorageIndex;
67  DenseMatrixType;
69  PlainPermutationType;
70  using Base::derived;
71  typedef Transpose<PermutationBase> TransposeReturnType;
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 
92  #ifndef EIGEN_PARSED_BY_DOXYGEN
93 
96  Derived& operator=(const PermutationBase& other)
97  {
98  indices() = other.indices();
99  return derived();
100  }
101  #endif
102 
104  inline Index rows() const { return Index(indices().size()); }
105 
107  inline Index cols() const { return Index(indices().size()); }
108 
110  inline Index size() const { return Index(indices().size()); }
111 
112  #ifndef EIGEN_PARSED_BY_DOXYGEN
113  template<typename DenseDerived>
114  void evalTo(MatrixBase<DenseDerived>& other) const
115  {
116  other.setZero();
117  for (Index i=0; i<rows(); ++i)
118  other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
119  }
120  #endif
121 
126  DenseMatrixType toDenseMatrix() const
127  {
128  return derived();
129  }
130 
132  const IndicesType& indices() const { return derived().indices(); }
134  IndicesType& indices() { return derived().indices(); }
135 
138  inline void resize(Index newSize)
139  {
140  indices().resize(newSize);
141  }
142 
144  void setIdentity()
145  {
146  StorageIndex n = StorageIndex(size());
147  for(StorageIndex i = 0; i < n; ++i)
148  indices().coeffRef(i) = i;
149  }
150 
153  void setIdentity(Index newSize)
154  {
155  resize(newSize);
156  setIdentity();
157  }
158 
169  {
170  eigen_assert(i>=0 && j>=0 && i<size() && j<size());
171  for(Index k = 0; k < size(); ++k)
172  {
173  if(indices().coeff(k) == i) indices().coeffRef(k) = StorageIndex(j);
174  else if(indices().coeff(k) == j) indices().coeffRef(k) = StorageIndex(i);
175  }
176  return derived();
177  }
178 
188  {
189  eigen_assert(i>=0 && j>=0 && i<size() && j<size());
190  std::swap(indices().coeffRef(i), indices().coeffRef(j));
191  return derived();
192  }
193 
198  inline TransposeReturnType inverse() const
199  { return TransposeReturnType(derived()); }
204  inline TransposeReturnType transpose() const
205  { return TransposeReturnType(derived()); }
206 
207  /**** multiplication helpers to hopefully get RVO ****/
208 
209 
210 #ifndef EIGEN_PARSED_BY_DOXYGEN
211  protected:
212  template<typename OtherDerived>
213  void assignTranspose(const PermutationBase<OtherDerived>& other)
214  {
215  for (Index i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i;
216  }
217  template<typename Lhs,typename Rhs>
218  void assignProduct(const Lhs& lhs, const Rhs& rhs)
219  {
220  eigen_assert(lhs.cols() == rhs.rows());
221  for (Index i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
222  }
223 #endif
224 
225  public:
226 
231  template<typename Other>
232  inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
233  { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
234 
239  template<typename Other>
240  inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other) const
241  { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
242 
247  template<typename Other> friend
248  inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other, const PermutationBase& perm)
249  { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
250 
256  {
257  Index res = 1;
258  Index n = size();
260  mask.fill(false);
261  Index r = 0;
262  while(r < n)
263  {
264  // search for the next seed
265  while(r<n && mask[r]) r++;
266  if(r>=n)
267  break;
268  // we got one, let's follow it until we are back to the seed
269  Index k0 = r++;
270  mask.coeffRef(k0) = true;
271  for(Index k=indices().coeff(k0); k!=k0; k=indices().coeff(k))
272  {
273  mask.coeffRef(k) = true;
274  res = -res;
275  }
276  }
277  return res;
278  }
279 
280  protected:
281 
282 };
283 
298 namespace internal {
299 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex>
300 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> >
301  : traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
302 {
303  typedef PermutationStorage StorageKind;
304  typedef Matrix<_StorageIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
305  typedef _StorageIndex StorageIndex;
306 };
307 }
308 
309 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex>
310 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> >
311 {
313  typedef internal::traits<PermutationMatrix> Traits;
314  public:
315 
316  typedef const PermutationMatrix& Nested;
317 
318  #ifndef EIGEN_PARSED_BY_DOXYGEN
319  typedef typename Traits::IndicesType IndicesType;
320  typedef typename Traits::StorageIndex StorageIndex;
321  #endif
322 
323  inline PermutationMatrix()
324  {}
325 
328  explicit inline PermutationMatrix(Index size) : m_indices(size)
329  {
330  eigen_internal_assert(size <= NumTraits<StorageIndex>::highest());
331  }
332 
334  template<typename OtherDerived>
336  : m_indices(other.indices()) {}
337 
338  #ifndef EIGEN_PARSED_BY_DOXYGEN
339 
341  inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
342  #endif
343 
351  template<typename Other>
352  explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
353  {}
354 
356  template<typename Other>
357  explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
358  : m_indices(tr.size())
359  {
360  *this = tr;
361  }
362 
364  template<typename Other>
366  {
367  m_indices = other.indices();
368  return *this;
369  }
370 
372  template<typename Other>
373  PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
374  {
375  return Base::operator=(tr.derived());
376  }
377 
378  #ifndef EIGEN_PARSED_BY_DOXYGEN
379 
383  {
384  m_indices = other.m_indices;
385  return *this;
386  }
387  #endif
388 
390  const IndicesType& indices() const { return m_indices; }
392  IndicesType& indices() { return m_indices; }
393 
394 
395  /**** multiplication helpers to hopefully get RVO ****/
396 
397 #ifndef EIGEN_PARSED_BY_DOXYGEN
398  template<typename Other>
400  : m_indices(other.nestedExpression().size())
401  {
402  eigen_internal_assert(m_indices.size() <= NumTraits<StorageIndex>::highest());
403  StorageIndex end = StorageIndex(m_indices.size());
404  for (StorageIndex i=0; i<end;++i)
405  m_indices.coeffRef(other.nestedExpression().indices().coeff(i)) = i;
406  }
407  template<typename Lhs,typename Rhs>
408  PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
409  : m_indices(lhs.indices().size())
410  {
411  Base::assignProduct(lhs,rhs);
412  }
413 #endif
414 
415  protected:
416 
417  IndicesType m_indices;
418 };
419 
420 
421 namespace internal {
422 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess>
423 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> >
424  : traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
425 {
426  typedef PermutationStorage StorageKind;
427  typedef Map<const Matrix<_StorageIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType;
428  typedef _StorageIndex StorageIndex;
429 };
430 }
431 
432 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess>
433 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess>
434  : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> >
435 {
436  typedef PermutationBase<Map> Base;
437  typedef internal::traits<Map> Traits;
438  public:
439 
440  #ifndef EIGEN_PARSED_BY_DOXYGEN
441  typedef typename Traits::IndicesType IndicesType;
442  typedef typename IndicesType::Scalar StorageIndex;
443  #endif
444 
445  inline Map(const StorageIndex* indicesPtr)
446  : m_indices(indicesPtr)
447  {}
448 
449  inline Map(const StorageIndex* indicesPtr, Index size)
450  : m_indices(indicesPtr,size)
451  {}
452 
454  template<typename Other>
455  Map& operator=(const PermutationBase<Other>& other)
456  { return Base::operator=(other.derived()); }
457 
459  template<typename Other>
460  Map& operator=(const TranspositionsBase<Other>& tr)
461  { return Base::operator=(tr.derived()); }
462 
463  #ifndef EIGEN_PARSED_BY_DOXYGEN
464 
467  Map& operator=(const Map& other)
468  {
469  m_indices = other.m_indices;
470  return *this;
471  }
472  #endif
473 
475  const IndicesType& indices() const { return m_indices; }
477  IndicesType& indices() { return m_indices; }
478 
479  protected:
480 
481  IndicesType m_indices;
482 };
483 
496 template<typename _IndicesType> class TranspositionsWrapper;
497 namespace internal {
498 template<typename _IndicesType>
499 struct traits<PermutationWrapper<_IndicesType> >
500 {
501  typedef PermutationStorage StorageKind;
502  typedef typename _IndicesType::Scalar Scalar;
503  typedef typename _IndicesType::Scalar StorageIndex;
504  typedef _IndicesType IndicesType;
505  enum {
506  RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
507  ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
508  MaxRowsAtCompileTime = IndicesType::MaxSizeAtCompileTime,
509  MaxColsAtCompileTime = IndicesType::MaxSizeAtCompileTime,
510  Flags = 0
511  };
512 };
513 }
514 
515 template<typename _IndicesType>
516 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
517 {
519  typedef internal::traits<PermutationWrapper> Traits;
520  public:
521 
522  #ifndef EIGEN_PARSED_BY_DOXYGEN
523  typedef typename Traits::IndicesType IndicesType;
524  #endif
525 
526  inline PermutationWrapper(const IndicesType& indices)
527  : m_indices(indices)
528  {}
529 
531  const typename internal::remove_all<typename IndicesType::Nested>::type&
532  indices() const { return m_indices; }
533 
534  protected:
535 
536  typename IndicesType::Nested m_indices;
537 };
538 
539 
542 template<typename MatrixDerived, typename PermutationDerived>
543 EIGEN_DEVICE_FUNC
544 const Product<MatrixDerived, PermutationDerived, AliasFreeProduct>
545 operator*(const MatrixBase<MatrixDerived> &matrix,
546  const PermutationBase<PermutationDerived>& permutation)
547 {
548  return Product<MatrixDerived, PermutationDerived, AliasFreeProduct>
549  (matrix.derived(), permutation.derived());
550 }
551 
554 template<typename PermutationDerived, typename MatrixDerived>
555 EIGEN_DEVICE_FUNC
556 const Product<PermutationDerived, MatrixDerived, AliasFreeProduct>
557 operator*(const PermutationBase<PermutationDerived> &permutation,
558  const MatrixBase<MatrixDerived>& matrix)
559 {
560  return Product<PermutationDerived, MatrixDerived, AliasFreeProduct>
561  (permutation.derived(), matrix.derived());
562 }
563 
564 namespace internal {
565 
566 /* Template partial specialization for transposed/inverse permutations */
567 
568 template<typename Derived>
569 struct traits<Transpose<PermutationBase<Derived> > >
570  : traits<Derived>
571 {};
572 
573 } // end namespace internal
574 
575 // TODO: the specificties should be handled by the evaluator,
576 // at the very least we should only specialize TransposeImpl
577 template<typename Derived>
578 class Transpose<PermutationBase<Derived> >
579  : public EigenBase<Transpose<PermutationBase<Derived> > >
580 {
581  typedef Derived PermutationType;
582  typedef typename PermutationType::IndicesType IndicesType;
583  typedef typename PermutationType::PlainPermutationType PlainPermutationType;
584  public:
585 
586  #ifndef EIGEN_PARSED_BY_DOXYGEN
587  typedef internal::traits<PermutationType> Traits;
588  typedef typename Derived::DenseMatrixType DenseMatrixType;
589  enum {
590  Flags = Traits::Flags,
591  RowsAtCompileTime = Traits::RowsAtCompileTime,
592  ColsAtCompileTime = Traits::ColsAtCompileTime,
593  MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
594  MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
595  };
596  typedef typename Traits::Scalar Scalar;
597  typedef typename Traits::StorageIndex StorageIndex;
598  #endif
599 
600  Transpose(const PermutationType& p) : m_permutation(p) {}
601 
602  inline Index rows() const { return m_permutation.rows(); }
603  inline Index cols() const { return m_permutation.cols(); }
604 
605  #ifndef EIGEN_PARSED_BY_DOXYGEN
606  template<typename DenseDerived>
607  void evalTo(MatrixBase<DenseDerived>& other) const
608  {
609  other.setZero();
610  for (Index i=0; i<rows();++i)
611  other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1);
612  }
613  #endif
614 
616  PlainPermutationType eval() const { return *this; }
617 
618  DenseMatrixType toDenseMatrix() const { return *this; }
619 
622  template<typename OtherDerived> friend
623  const Product<OtherDerived, Transpose, AliasFreeProduct>
624  operator*(const MatrixBase<OtherDerived>& matrix, const Transpose& trPerm)
625  {
626  return Product<OtherDerived, Transpose, AliasFreeProduct>(matrix.derived(), trPerm.derived());
627  }
628 
631  template<typename OtherDerived>
632  const Product<Transpose, OtherDerived, AliasFreeProduct>
633  operator*(const MatrixBase<OtherDerived>& matrix) const
634  {
635  return Product<Transpose, OtherDerived, AliasFreeProduct>(*this, matrix.derived());
636  }
637 
638  const PermutationType& nestedExpression() const { return m_permutation; }
639 
640  protected:
641  const PermutationType& m_permutation;
642 };
643 
644 template<typename Derived>
645 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
646 {
647  return derived();
648 }
649 
650 namespace internal {
651 
652 template<> struct AssignmentKind<DenseShape,PermutationShape> { typedef EigenBase2EigenBase Kind; };
653 
654 } // end namespace internal
655 
656 } // end namespace Eigen
657 
658 #endif // EIGEN_PERMUTATIONMATRIX_H
PermutationMatrix(const MatrixBase< Other > &indices)
Definition: PermutationMatrix.h:352
TransposeReturnType transpose() const
Definition: PermutationMatrix.h:204
const internal::remove_all< typename IndicesType::Nested >::type & indices() const
Definition: PermutationMatrix.h:532
Expression of the transpose of a matrix.
Definition: Transpose.h:53
Derived & applyTranspositionOnTheRight(Index i, Index j)
Definition: PermutationMatrix.h:187
Index size() const
Definition: PermutationMatrix.h:110
Definition: LDLT.h:16
const IndicesType & indices() const
Definition: PermutationMatrix.h:132
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
friend PlainPermutationType operator*(const Transpose< PermutationBase< Other > > &other, const PermutationBase &perm)
Definition: PermutationMatrix.h:248
Derived & derived()
Definition: EigenBase.h:44
Derived & operator=(const PermutationBase< OtherDerived > &other)
Definition: PermutationMatrix.h:76
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:37
PermutationMatrix & operator=(const PermutationBase< Other > &other)
Definition: PermutationMatrix.h:365
Map(PointerArgType dataPtr, const StrideType &stride=StrideType())
Definition: Map.h:123
Index rows() const
Definition: PermutationMatrix.h:104
Base class for permutations.
Definition: PermutationMatrix.h:50
Definition: EigenBase.h:28
Permutation matrix.
Definition: PermutationMatrix.h:310
PermutationMatrix & operator=(const TranspositionsBase< Other > &tr)
Definition: PermutationMatrix.h:373
IndicesType & indices()
Definition: PermutationMatrix.h:134
PlainPermutationType operator*(const PermutationBase< Other > &other) const
Definition: PermutationMatrix.h:232
void fill(const Scalar &value)
Definition: CwiseNullaryOp.h:326
Derived & operator=(const TranspositionsBase< OtherDerived > &tr)
Definition: PermutationMatrix.h:84
DenseMatrixType toDenseMatrix() const
Definition: PermutationMatrix.h:126
Index determinant() const
Definition: PermutationMatrix.h:255
Class to view a vector of integers as a permutation matrix.
Definition: PermutationMatrix.h:516
void setIdentity()
Definition: PermutationMatrix.h:144
PermutationMatrix(const TranspositionsBase< Other > &tr)
Definition: PermutationMatrix.h:357
Definition: Eigen_Colamd.h:54
Derived & setZero()
Definition: CwiseNullaryOp.h:504
Derived & applyTranspositionOnTheLeft(Index i, Index j)
Definition: PermutationMatrix.h:168
PermutationMatrix(Index size)
Definition: PermutationMatrix.h:328
const IndicesType & indices() const
Definition: PermutationMatrix.h:390
IndicesType & indices()
Definition: PermutationMatrix.h:392
TransposeReturnType inverse() const
Definition: PermutationMatrix.h:198
void setIdentity(Index newSize)
Definition: PermutationMatrix.h:153
PlainPermutationType operator*(const Transpose< PermutationBase< Other > > &other) const
Definition: PermutationMatrix.h:240
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
void resize(Index newSize)
Definition: PermutationMatrix.h:138
const internal::remove_all< typename MatrixType::Nested >::type & nestedExpression() const
Definition: Transpose.h:73
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Index cols() const
Definition: PermutationMatrix.h:107
PermutationMatrix(const PermutationBase< OtherDerived > &other)
Definition: PermutationMatrix.h:335