Eigen  3.2.1
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 template<int RowCol,typename IndicesType,typename MatrixType, typename StorageKind> class PermutedImpl;
17 
42 namespace internal {
43 
44 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
45 struct permut_matrix_product_retval;
46 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
47 struct permut_sparsematrix_product_retval;
48 enum PermPermProduct_t {PermPermProduct};
49 
50 } // end namespace internal
51 
52 template<typename Derived>
53 class PermutationBase : public EigenBase<Derived>
54 {
55  typedef internal::traits<Derived> Traits;
56  typedef EigenBase<Derived> Base;
57  public:
58 
59  #ifndef EIGEN_PARSED_BY_DOXYGEN
60  typedef typename Traits::IndicesType IndicesType;
61  enum {
62  Flags = Traits::Flags,
63  CoeffReadCost = Traits::CoeffReadCost,
64  RowsAtCompileTime = Traits::RowsAtCompileTime,
65  ColsAtCompileTime = Traits::ColsAtCompileTime,
66  MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
67  MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
68  };
69  typedef typename Traits::Scalar Scalar;
70  typedef typename Traits::Index Index;
72  DenseMatrixType;
74  PlainPermutationType;
75  using Base::derived;
76  #endif
77 
79  template<typename OtherDerived>
81  {
82  indices() = other.indices();
83  return derived();
84  }
85 
87  template<typename OtherDerived>
88  Derived& operator=(const TranspositionsBase<OtherDerived>& tr)
89  {
90  setIdentity(tr.size());
91  for(Index k=size()-1; k>=0; --k)
92  applyTranspositionOnTheRight(k,tr.coeff(k));
93  return derived();
94  }
95 
96  #ifndef EIGEN_PARSED_BY_DOXYGEN
97 
100  Derived& operator=(const PermutationBase& other)
101  {
102  indices() = other.indices();
103  return derived();
104  }
105  #endif
106 
108  inline Index rows() const { return Index(indices().size()); }
109 
111  inline Index cols() const { return Index(indices().size()); }
112 
114  inline Index size() const { return Index(indices().size()); }
115 
116  #ifndef EIGEN_PARSED_BY_DOXYGEN
117  template<typename DenseDerived>
118  void evalTo(MatrixBase<DenseDerived>& other) const
119  {
120  other.setZero();
121  for (int i=0; i<rows();++i)
122  other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
123  }
124  #endif
125 
130  DenseMatrixType toDenseMatrix() const
131  {
132  return derived();
133  }
134 
136  const IndicesType& indices() const { return derived().indices(); }
138  IndicesType& indices() { return derived().indices(); }
139 
142  inline void resize(Index newSize)
143  {
144  indices().resize(newSize);
145  }
146 
148  void setIdentity()
149  {
150  for(Index i = 0; i < size(); ++i)
151  indices().coeffRef(i) = i;
152  }
153 
156  void setIdentity(Index newSize)
157  {
158  resize(newSize);
159  setIdentity();
160  }
161 
171  Derived& applyTranspositionOnTheLeft(Index i, Index j)
172  {
173  eigen_assert(i>=0 && j>=0 && i<size() && j<size());
174  for(Index k = 0; k < size(); ++k)
175  {
176  if(indices().coeff(k) == i) indices().coeffRef(k) = j;
177  else if(indices().coeff(k) == j) indices().coeffRef(k) = i;
178  }
179  return derived();
180  }
181 
190  Derived& applyTranspositionOnTheRight(Index i, Index j)
191  {
192  eigen_assert(i>=0 && j>=0 && i<size() && j<size());
193  std::swap(indices().coeffRef(i), indices().coeffRef(j));
194  return derived();
195  }
196 
202  { return derived(); }
208  { return derived(); }
209 
210  /**** multiplication helpers to hopefully get RVO ****/
211 
212 
213 #ifndef EIGEN_PARSED_BY_DOXYGEN
214  protected:
215  template<typename OtherDerived>
216  void assignTranspose(const PermutationBase<OtherDerived>& other)
217  {
218  for (int i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i;
219  }
220  template<typename Lhs,typename Rhs>
221  void assignProduct(const Lhs& lhs, const Rhs& rhs)
222  {
223  eigen_assert(lhs.cols() == rhs.rows());
224  for (int i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
225  }
226 #endif
227 
228  public:
229 
234  template<typename Other>
235  inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
236  { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
237 
242  template<typename Other>
243  inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other) const
244  { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
245 
250  template<typename Other> friend
251  inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other, const PermutationBase& perm)
252  { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
253 
254  protected:
255 
256 };
257 
272 namespace internal {
273 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
274 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
275  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
276 {
277  typedef IndexType Index;
278  typedef Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
279 };
280 }
281 
282 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
283 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
284 {
286  typedef internal::traits<PermutationMatrix> Traits;
287  public:
288 
289  #ifndef EIGEN_PARSED_BY_DOXYGEN
290  typedef typename Traits::IndicesType IndicesType;
291  #endif
292 
293  inline PermutationMatrix()
294  {}
295 
298  inline PermutationMatrix(int size) : m_indices(size)
299  {}
300 
302  template<typename OtherDerived>
304  : m_indices(other.indices()) {}
305 
306  #ifndef EIGEN_PARSED_BY_DOXYGEN
307 
309  inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
310  #endif
311 
319  template<typename Other>
320  explicit inline PermutationMatrix(const MatrixBase<Other>& a_indices) : m_indices(a_indices)
321  {}
322 
324  template<typename Other>
325  explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
326  : m_indices(tr.size())
327  {
328  *this = tr;
329  }
330 
332  template<typename Other>
334  {
335  m_indices = other.indices();
336  return *this;
337  }
338 
340  template<typename Other>
341  PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
342  {
343  return Base::operator=(tr.derived());
344  }
345 
346  #ifndef EIGEN_PARSED_BY_DOXYGEN
347 
351  {
352  m_indices = other.m_indices;
353  return *this;
354  }
355  #endif
356 
358  const IndicesType& indices() const { return m_indices; }
360  IndicesType& indices() { return m_indices; }
361 
362 
363  /**** multiplication helpers to hopefully get RVO ****/
364 
365 #ifndef EIGEN_PARSED_BY_DOXYGEN
366  template<typename Other>
368  : m_indices(other.nestedPermutation().size())
369  {
370  for (int i=0; i<m_indices.size();++i) m_indices.coeffRef(other.nestedPermutation().indices().coeff(i)) = i;
371  }
372  template<typename Lhs,typename Rhs>
373  PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
374  : m_indices(lhs.indices().size())
375  {
376  Base::assignProduct(lhs,rhs);
377  }
378 #endif
379 
380  protected:
381 
382  IndicesType m_indices;
383 };
384 
385 
386 namespace internal {
387 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
388 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
389  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
390 {
391  typedef IndexType Index;
392  typedef Map<const Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType;
393 };
394 }
395 
396 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
397 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess>
398  : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
399 {
400  typedef PermutationBase<Map> Base;
401  typedef internal::traits<Map> Traits;
402  public:
403 
404  #ifndef EIGEN_PARSED_BY_DOXYGEN
405  typedef typename Traits::IndicesType IndicesType;
406  typedef typename IndicesType::Scalar Index;
407  #endif
408 
409  inline Map(const Index* indicesPtr)
410  : m_indices(indicesPtr)
411  {}
412 
413  inline Map(const Index* indicesPtr, Index size)
414  : m_indices(indicesPtr,size)
415  {}
416 
418  template<typename Other>
419  Map& operator=(const PermutationBase<Other>& other)
420  { return Base::operator=(other.derived()); }
421 
423  template<typename Other>
424  Map& operator=(const TranspositionsBase<Other>& tr)
425  { return Base::operator=(tr.derived()); }
426 
427  #ifndef EIGEN_PARSED_BY_DOXYGEN
428 
431  Map& operator=(const Map& other)
432  {
433  m_indices = other.m_indices;
434  return *this;
435  }
436  #endif
437 
439  const IndicesType& indices() const { return m_indices; }
441  IndicesType& indices() { return m_indices; }
442 
443  protected:
444 
445  IndicesType m_indices;
446 };
447 
460 struct PermutationStorage {};
461 
462 template<typename _IndicesType> class TranspositionsWrapper;
463 namespace internal {
464 template<typename _IndicesType>
465 struct traits<PermutationWrapper<_IndicesType> >
466 {
467  typedef PermutationStorage StorageKind;
468  typedef typename _IndicesType::Scalar Scalar;
469  typedef typename _IndicesType::Scalar Index;
470  typedef _IndicesType IndicesType;
471  enum {
472  RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
473  ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
474  MaxRowsAtCompileTime = IndicesType::MaxRowsAtCompileTime,
475  MaxColsAtCompileTime = IndicesType::MaxColsAtCompileTime,
476  Flags = 0,
477  CoeffReadCost = _IndicesType::CoeffReadCost
478  };
479 };
480 }
481 
482 template<typename _IndicesType>
483 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
484 {
486  typedef internal::traits<PermutationWrapper> Traits;
487  public:
488 
489  #ifndef EIGEN_PARSED_BY_DOXYGEN
490  typedef typename Traits::IndicesType IndicesType;
491  #endif
492 
493  inline PermutationWrapper(const IndicesType& a_indices)
494  : m_indices(a_indices)
495  {}
496 
498  const typename internal::remove_all<typename IndicesType::Nested>::type&
499  indices() const { return m_indices; }
500 
501  protected:
502 
503  typename IndicesType::Nested m_indices;
504 };
505 
508 template<typename Derived, typename PermutationDerived>
509 inline const internal::permut_matrix_product_retval<PermutationDerived, Derived, OnTheRight>
511  const PermutationBase<PermutationDerived> &permutation)
512 {
513  return internal::permut_matrix_product_retval
514  <PermutationDerived, Derived, OnTheRight>
515  (permutation.derived(), matrix.derived());
516 }
517 
520 template<typename Derived, typename PermutationDerived>
521 inline const internal::permut_matrix_product_retval
522  <PermutationDerived, Derived, OnTheLeft>
524  const MatrixBase<Derived>& matrix)
525 {
526  return internal::permut_matrix_product_retval
527  <PermutationDerived, Derived, OnTheLeft>
528  (permutation.derived(), matrix.derived());
529 }
530 
531 namespace internal {
532 
533 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
534 struct traits<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
535 {
536  typedef typename MatrixType::PlainObject ReturnType;
537 };
538 
539 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
540 struct permut_matrix_product_retval
541  : public ReturnByValue<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
542 {
543  typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
544  typedef typename MatrixType::Index Index;
545 
546  permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
547  : m_permutation(perm), m_matrix(matrix)
548  {}
549 
550  inline Index rows() const { return m_matrix.rows(); }
551  inline Index cols() const { return m_matrix.cols(); }
552 
553  template<typename Dest> inline void evalTo(Dest& dst) const
554  {
555  const Index n = Side==OnTheLeft ? rows() : cols();
556  // FIXME we need an is_same for expression that is not sensitive to constness. For instance
557  // is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true.
558  if(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix))
559  {
560  // apply the permutation inplace
561  Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
562  mask.fill(false);
563  Index r = 0;
564  while(r < m_permutation.size())
565  {
566  // search for the next seed
567  while(r<m_permutation.size() && mask[r]) r++;
568  if(r>=m_permutation.size())
569  break;
570  // we got one, let's follow it until we are back to the seed
571  Index k0 = r++;
572  Index kPrev = k0;
573  mask.coeffRef(k0) = true;
574  for(Index k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
575  {
576  Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
577  .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
578  (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
579 
580  mask.coeffRef(k) = true;
581  kPrev = k;
582  }
583  }
584  }
585  else
586  {
587  for(int i = 0; i < n; ++i)
588  {
589  Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
590  (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
591 
592  =
593 
594  Block<const MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
595  (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
596  }
597  }
598  }
599 
600  protected:
601  const PermutationType& m_permutation;
602  typename MatrixType::Nested m_matrix;
603 };
604 
605 /* Template partial specialization for transposed/inverse permutations */
606 
607 template<typename Derived>
608 struct traits<Transpose<PermutationBase<Derived> > >
609  : traits<Derived>
610 {};
611 
612 } // end namespace internal
613 
614 template<typename Derived>
615 class Transpose<PermutationBase<Derived> >
616  : public EigenBase<Transpose<PermutationBase<Derived> > >
617 {
618  typedef Derived PermutationType;
619  typedef typename PermutationType::IndicesType IndicesType;
620  typedef typename PermutationType::PlainPermutationType PlainPermutationType;
621  public:
622 
623  #ifndef EIGEN_PARSED_BY_DOXYGEN
624  typedef internal::traits<PermutationType> Traits;
625  typedef typename Derived::DenseMatrixType DenseMatrixType;
626  enum {
627  Flags = Traits::Flags,
628  CoeffReadCost = Traits::CoeffReadCost,
629  RowsAtCompileTime = Traits::RowsAtCompileTime,
630  ColsAtCompileTime = Traits::ColsAtCompileTime,
631  MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
632  MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
633  };
634  typedef typename Traits::Scalar Scalar;
635  #endif
636 
637  Transpose(const PermutationType& p) : m_permutation(p) {}
638 
639  inline int rows() const { return m_permutation.rows(); }
640  inline int cols() const { return m_permutation.cols(); }
641 
642  #ifndef EIGEN_PARSED_BY_DOXYGEN
643  template<typename DenseDerived>
644  void evalTo(MatrixBase<DenseDerived>& other) const
645  {
646  other.setZero();
647  for (int i=0; i<rows();++i)
648  other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1);
649  }
650  #endif
651 
653  PlainPermutationType eval() const { return *this; }
654 
655  DenseMatrixType toDenseMatrix() const { return *this; }
656 
659  template<typename OtherDerived> friend
660  inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>
661  operator*(const MatrixBase<OtherDerived>& matrix, const Transpose& trPerm)
662  {
663  return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>(trPerm.m_permutation, matrix.derived());
664  }
665 
668  template<typename OtherDerived>
669  inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>
670  operator*(const MatrixBase<OtherDerived>& matrix) const
671  {
672  return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>(m_permutation, matrix.derived());
673  }
674 
675  const PermutationType& nestedPermutation() const { return m_permutation; }
676 
677  protected:
678  const PermutationType& m_permutation;
679 };
680 
681 template<typename Derived>
682 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
683 {
684  return derived();
685 }
686 
687 } // end namespace Eigen
688 
689 #endif // EIGEN_PERMUTATIONMATRIX_H
void setIdentity()
Definition: PermutationMatrix.h:148
const IndicesType & indices() const
Definition: PermutationMatrix.h:358
Derived & applyTranspositionOnTheRight(Index i, Index j)
Definition: PermutationMatrix.h:190
Transpose< PermutationBase > inverse() const
Definition: PermutationMatrix.h:201
IndicesType & indices()
Definition: PermutationMatrix.h:138
Expression of the transpose of a matrix.
Definition: Transpose.h:57
Derived & setZero()
Definition: CwiseNullaryOp.h:499
Index rows() const
Definition: PermutationMatrix.h:108
const internal::permut_matrix_product_retval< PermutationDerived, Derived, OnTheRight > operator*(const MatrixBase< Derived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:510
friend PlainPermutationType operator*(const Transpose< PermutationBase< Other > > &other, const PermutationBase &perm)
Definition: PermutationMatrix.h:251
PermutationMatrix(const MatrixBase< Other > &a_indices)
Definition: PermutationMatrix.h:320
Derived & applyTranspositionOnTheLeft(Index i, Index j)
Definition: PermutationMatrix.h:171
PermutationMatrix & operator=(const TranspositionsBase< Other > &tr)
Definition: PermutationMatrix.h:341
PermutationMatrix(const PermutationBase< OtherDerived > &other)
Definition: PermutationMatrix.h:303
Base class for permutations.
Definition: PermutationMatrix.h:53
Definition: EigenBase.h:26
DenseMatrixType toDenseMatrix() const
Definition: PermutationMatrix.h:130
void resize(Index newSize)
Definition: PermutationMatrix.h:142
Permutation matrix.
Definition: PermutationMatrix.h:283
Derived & operator=(const TranspositionsBase< OtherDerived > &tr)
Definition: PermutationMatrix.h:88
Index size() const
Definition: PermutationMatrix.h:114
PermutationMatrix & operator=(const PermutationBase< Other > &other)
Definition: PermutationMatrix.h:333
Definition: Constants.h:279
Transpose< PermutationBase > transpose() const
Definition: PermutationMatrix.h:207
Derived & operator=(const PermutationBase< OtherDerived > &other)
Definition: PermutationMatrix.h:80
Derived & derived()
Definition: EigenBase.h:34
Index cols() const
Definition: PermutationMatrix.h:111
PermutationMatrix(int size)
Definition: PermutationMatrix.h:298
Class to view a vector of integers as a permutation matrix.
Definition: PermutationMatrix.h:483
Map(PointerArgType dataPtr, const StrideType &a_stride=StrideType())
Definition: Map.h:139
PlainPermutationType operator*(const PermutationBase< Other > &other) const
Definition: PermutationMatrix.h:235
void setIdentity(Index newSize)
Definition: PermutationMatrix.h:156
IndicesType & indices()
Definition: PermutationMatrix.h:360
const internal::remove_all< typename IndicesType::Nested >::type & indices() const
Definition: PermutationMatrix.h:499
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:127
Definition: Constants.h:277
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
PlainPermutationType operator*(const Transpose< PermutationBase< Other > > &other) const
Definition: PermutationMatrix.h:243
const IndicesType & indices() const
Definition: PermutationMatrix.h:136
PermutationMatrix(const TranspositionsBase< Other > &tr)
Definition: PermutationMatrix.h:325