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 indices().size(); }
109 
111  inline Index cols() const { return indices().size(); }
112 
114  inline Index size() const { return 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 size)
143  {
144  indices().resize(size);
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 size)
157  {
158  resize(size);
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>& indices) : m_indices(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* indices)
410  : m_indices(indices)
411  {}
412 
413  inline Map(const Index* indices, Index size)
414  : m_indices(indices,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& indices)
494  : m_indices(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 
545  permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
546  : m_permutation(perm), m_matrix(matrix)
547  {}
548 
549  inline int rows() const { return m_matrix.rows(); }
550  inline int cols() const { return m_matrix.cols(); }
551 
552  template<typename Dest> inline void evalTo(Dest& dst) const
553  {
554  const int n = Side==OnTheLeft ? rows() : cols();
555 
556  if(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix))
557  {
558  // apply the permutation inplace
559  Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
560  mask.fill(false);
561  int r = 0;
562  while(r < m_permutation.size())
563  {
564  // search for the next seed
565  while(r<m_permutation.size() && mask[r]) r++;
566  if(r>=m_permutation.size())
567  break;
568  // we got one, let's follow it until we are back to the seed
569  int k0 = r++;
570  int kPrev = k0;
571  mask.coeffRef(k0) = true;
572  for(int k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
573  {
574  Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
575  .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
576  (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
577 
578  mask.coeffRef(k) = true;
579  kPrev = k;
580  }
581  }
582  }
583  else
584  {
585  for(int i = 0; i < n; ++i)
586  {
587  Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
588  (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
589 
590  =
591 
592  Block<const MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
593  (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
594  }
595  }
596  }
597 
598  protected:
599  const PermutationType& m_permutation;
600  typename MatrixType::Nested m_matrix;
601 };
602 
603 /* Template partial specialization for transposed/inverse permutations */
604 
605 template<typename Derived>
606 struct traits<Transpose<PermutationBase<Derived> > >
607  : traits<Derived>
608 {};
609 
610 } // end namespace internal
611 
612 template<typename Derived>
613 class Transpose<PermutationBase<Derived> >
614  : public EigenBase<Transpose<PermutationBase<Derived> > >
615 {
616  typedef Derived PermutationType;
617  typedef typename PermutationType::IndicesType IndicesType;
618  typedef typename PermutationType::PlainPermutationType PlainPermutationType;
619  public:
620 
621  #ifndef EIGEN_PARSED_BY_DOXYGEN
622  typedef internal::traits<PermutationType> Traits;
623  typedef typename Derived::DenseMatrixType DenseMatrixType;
624  enum {
625  Flags = Traits::Flags,
626  CoeffReadCost = Traits::CoeffReadCost,
627  RowsAtCompileTime = Traits::RowsAtCompileTime,
628  ColsAtCompileTime = Traits::ColsAtCompileTime,
629  MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
630  MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
631  };
632  typedef typename Traits::Scalar Scalar;
633  #endif
634 
635  Transpose(const PermutationType& p) : m_permutation(p) {}
636 
637  inline int rows() const { return m_permutation.rows(); }
638  inline int cols() const { return m_permutation.cols(); }
639 
640  #ifndef EIGEN_PARSED_BY_DOXYGEN
641  template<typename DenseDerived>
642  void evalTo(MatrixBase<DenseDerived>& other) const
643  {
644  other.setZero();
645  for (int i=0; i<rows();++i)
646  other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1);
647  }
648  #endif
649 
651  PlainPermutationType eval() const { return *this; }
652 
653  DenseMatrixType toDenseMatrix() const { return *this; }
654 
657  template<typename OtherDerived> friend
658  inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>
659  operator*(const MatrixBase<OtherDerived>& matrix, const Transpose& trPerm)
660  {
661  return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>(trPerm.m_permutation, matrix.derived());
662  }
663 
666  template<typename OtherDerived>
667  inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>
668  operator*(const MatrixBase<OtherDerived>& matrix) const
669  {
670  return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>(m_permutation, matrix.derived());
671  }
672 
673  const PermutationType& nestedPermutation() const { return m_permutation; }
674 
675  protected:
676  const PermutationType& m_permutation;
677 };
678 
679 template<typename Derived>
680 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
681 {
682  return derived();
683 }
684 
685 } // end namespace Eigen
686 
687 #endif // EIGEN_PERMUTATIONMATRIX_H