PermutationMatrix.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25 
26 #ifndef EIGEN_PERMUTATIONMATRIX_H
27 #define EIGEN_PERMUTATIONMATRIX_H
28 
29 namespace Eigen {
30 
31 template<int RowCol,typename IndicesType,typename MatrixType, typename StorageKind> class PermutedImpl;
32 
57 namespace internal {
58 
59 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
60 struct permut_matrix_product_retval;
61 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
62 struct permut_sparsematrix_product_retval;
64 
65 } // end namespace internal
66 
67 template<typename Derived>
68 class PermutationBase : public EigenBase<Derived>
69 {
70  typedef internal::traits<Derived> Traits;
71  typedef EigenBase<Derived> Base;
72  public:
73 
74  #ifndef EIGEN_PARSED_BY_DOXYGEN
75  typedef typename Traits::IndicesType IndicesType;
76  enum {
77  Flags = Traits::Flags,
78  CoeffReadCost = Traits::CoeffReadCost,
79  RowsAtCompileTime = Traits::RowsAtCompileTime,
80  ColsAtCompileTime = Traits::ColsAtCompileTime,
81  MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
82  MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
83  };
84  typedef typename Traits::Scalar Scalar;
85  typedef typename Traits::Index Index;
87  DenseMatrixType;
89  PlainPermutationType;
90  using Base::derived;
91  #endif
92 
94  template<typename OtherDerived>
96  {
97  indices() = other.indices();
98  return derived();
99  }
100 
102  template<typename OtherDerived>
104  {
105  setIdentity(tr.size());
106  for(Index k=size()-1; k>=0; --k)
108  return derived();
109  }
110 
111  #ifndef EIGEN_PARSED_BY_DOXYGEN
112 
115  Derived& operator=(const PermutationBase& other)
116  {
117  indices() = other.indices();
118  return derived();
119  }
120  #endif
121 
123  inline Index rows() const { return indices().size(); }
124 
126  inline Index cols() const { return indices().size(); }
127 
129  inline Index size() const { return indices().size(); }
130 
131  #ifndef EIGEN_PARSED_BY_DOXYGEN
132  template<typename DenseDerived>
133  void evalTo(MatrixBase<DenseDerived>& other) const
134  {
135  other.setZero();
136  for (int i=0; i<rows();++i)
137  other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
138  }
139  #endif
140 
145  DenseMatrixType toDenseMatrix() const
146  {
147  return derived();
148  }
149 
151  const IndicesType& indices() const { return derived().indices(); }
153  IndicesType& indices() { return derived().indices(); }
154 
157  inline void resize(Index size)
158  {
159  indices().resize(size);
160  }
161 
163  void setIdentity()
164  {
165  for(Index i = 0; i < size(); ++i)
166  indices().coeffRef(i) = i;
167  }
168 
172  {
173  resize(size);
174  setIdentity();
175  }
176 
187  {
188  eigen_assert(i>=0 && j>=0 && i<size() && j<size());
189  for(Index k = 0; k < size(); ++k)
190  {
191  if(indices().coeff(k) == i) indices().coeffRef(k) = j;
192  else if(indices().coeff(k) == j) indices().coeffRef(k) = i;
193  }
194  return derived();
195  }
196 
206  {
207  eigen_assert(i>=0 && j>=0 && i<size() && j<size());
208  std::swap(indices().coeffRef(i), indices().coeffRef(j));
209  return derived();
210  }
211 
217  { return derived(); }
223  { return derived(); }
224 
225  /**** multiplication helpers to hopefully get RVO ****/
226 
227 
228 #ifndef EIGEN_PARSED_BY_DOXYGEN
229  protected:
230  template<typename OtherDerived>
231  void assignTranspose(const PermutationBase<OtherDerived>& other)
232  {
233  for (int i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i;
234  }
235  template<typename Lhs,typename Rhs>
236  void assignProduct(const Lhs& lhs, const Rhs& rhs)
237  {
238  eigen_assert(lhs.cols() == rhs.rows());
239  for (int i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
240  }
241 #endif
242 
243  public:
244 
249  template<typename Other>
250  inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
251  { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
252 
257  template<typename Other>
258  inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other) const
259  { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
260 
265  template<typename Other> friend
266  inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other, const PermutationBase& perm)
267  { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
268 
269  protected:
270 
271 };
272 
287 namespace internal {
288 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
289 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
290  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
291 {
292  typedef IndexType Index;
293  typedef Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
294 };
295 }
296 
297 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
298 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
299 {
301  typedef internal::traits<PermutationMatrix> Traits;
302  public:
303 
304  #ifndef EIGEN_PARSED_BY_DOXYGEN
305  typedef typename Traits::IndicesType IndicesType;
306  #endif
307 
309  {}
310 
313  inline PermutationMatrix(int size) : m_indices(size)
314  {}
315 
317  template<typename OtherDerived>
319  : m_indices(other.indices()) {}
320 
321  #ifndef EIGEN_PARSED_BY_DOXYGEN
322 
324  inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
325  #endif
326 
334  template<typename Other>
335  explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
336  {}
337 
339  template<typename Other>
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>
357  {
358  return Base::operator=(tr.derived());
359  }
360 
361  #ifndef EIGEN_PARSED_BY_DOXYGEN
362 
366  {
367  m_indices = other.m_indices;
368  return *this;
369  }
370  #endif
371 
373  const IndicesType& indices() const { return m_indices; }
375  IndicesType& indices() { return m_indices; }
376 
377 
378  /**** multiplication helpers to hopefully get RVO ****/
379 
380 #ifndef EIGEN_PARSED_BY_DOXYGEN
381  template<typename Other>
383  : m_indices(other.nestedPermutation().size())
384  {
385  for (int i=0; i<m_indices.size();++i) m_indices.coeffRef(other.nestedPermutation().indices().coeff(i)) = i;
386  }
387  template<typename Lhs,typename Rhs>
388  PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
389  : m_indices(lhs.indices().size())
390  {
391  Base::assignProduct(lhs,rhs);
392  }
393 #endif
394 
395  protected:
396 
397  IndicesType m_indices;
398 };
399 
400 
401 namespace internal {
402 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
403 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
404  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
405 {
406  typedef IndexType Index;
408 };
409 }
410 
411 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
412 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess>
413  : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
414 {
415  typedef PermutationBase<Map> Base;
416  typedef internal::traits<Map> Traits;
417  public:
418 
419  #ifndef EIGEN_PARSED_BY_DOXYGEN
420  typedef typename Traits::IndicesType IndicesType;
421  typedef typename IndicesType::Scalar Index;
422  #endif
423 
424  inline Map(const Index* indices)
425  : m_indices(indices)
426  {}
427 
428  inline Map(const Index* indices, Index size)
429  : m_indices(indices,size)
430  {}
431 
433  template<typename Other>
434  Map& operator=(const PermutationBase<Other>& other)
435  { return Base::operator=(other.derived()); }
436 
438  template<typename Other>
439  Map& operator=(const TranspositionsBase<Other>& tr)
440  { return Base::operator=(tr.derived()); }
441 
442  #ifndef EIGEN_PARSED_BY_DOXYGEN
443 
446  Map& operator=(const Map& other)
447  {
448  m_indices = other.m_indices;
449  return *this;
450  }
451  #endif
452 
454  const IndicesType& indices() const { return m_indices; }
456  IndicesType& indices() { return m_indices; }
457 
458  protected:
459 
460  IndicesType m_indices;
461 };
462 
476 
477 template<typename _IndicesType> class TranspositionsWrapper;
478 namespace internal {
479 template<typename _IndicesType>
480 struct traits<PermutationWrapper<_IndicesType> >
481 {
482  typedef PermutationStorage StorageKind;
483  typedef typename _IndicesType::Scalar Scalar;
484  typedef typename _IndicesType::Scalar Index;
485  typedef _IndicesType IndicesType;
486  enum {
487  RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
488  ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
489  MaxRowsAtCompileTime = IndicesType::MaxRowsAtCompileTime,
490  MaxColsAtCompileTime = IndicesType::MaxColsAtCompileTime,
491  Flags = 0,
492  CoeffReadCost = _IndicesType::CoeffReadCost
493  };
494 };
495 }
496 
497 template<typename _IndicesType>
498 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
499 {
501  typedef internal::traits<PermutationWrapper> Traits;
502  public:
503 
504  #ifndef EIGEN_PARSED_BY_DOXYGEN
505  typedef typename Traits::IndicesType IndicesType;
506  #endif
507 
508  inline PermutationWrapper(const IndicesType& indices)
509  : m_indices(indices)
510  {}
511 
513  const typename internal::remove_all<typename IndicesType::Nested>::type&
514  indices() const { return m_indices; }
515 
516  protected:
517 
518  typename IndicesType::Nested m_indices;
519 };
520 
523 template<typename Derived, typename PermutationDerived>
524 inline const internal::permut_matrix_product_retval<PermutationDerived, Derived, OnTheRight>
526  const PermutationBase<PermutationDerived> &permutation)
527 {
528  return internal::permut_matrix_product_retval
529  <PermutationDerived, Derived, OnTheRight>
530  (permutation.derived(), matrix.derived());
531 }
532 
535 template<typename Derived, typename PermutationDerived>
536 inline const internal::permut_matrix_product_retval
537  <PermutationDerived, Derived, OnTheLeft>
539  const MatrixBase<Derived>& matrix)
540 {
541  return internal::permut_matrix_product_retval
542  <PermutationDerived, Derived, OnTheLeft>
543  (permutation.derived(), matrix.derived());
544 }
545 
546 namespace internal {
547 
548 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
549 struct traits<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
550 {
551  typedef typename MatrixType::PlainObject ReturnType;
552 };
553 
554 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
555 struct permut_matrix_product_retval
556  : public ReturnByValue<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
557 {
558  typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
559 
560  permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
561  : m_permutation(perm), m_matrix(matrix)
562  {}
563 
564  inline int rows() const { return m_matrix.rows(); }
565  inline int cols() const { return m_matrix.cols(); }
566 
567  template<typename Dest> inline void evalTo(Dest& dst) const
568  {
569  const int n = Side==OnTheLeft ? rows() : cols();
570 
571  if(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix))
572  {
573  // apply the permutation inplace
574  Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
575  mask.fill(false);
576  int r = 0;
577  while(r < m_permutation.size())
578  {
579  // search for the next seed
580  while(r<m_permutation.size() && mask[r]) r++;
581  if(r>=m_permutation.size())
582  break;
583  // we got one, let's follow it until we are back to the seed
584  int k0 = r++;
585  int kPrev = k0;
586  mask.coeffRef(k0) = true;
587  for(int k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
588  {
589  Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
590  .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
591  (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
592 
593  mask.coeffRef(k) = true;
594  kPrev = k;
595  }
596  }
597  }
598  else
599  {
600  for(int i = 0; i < n; ++i)
601  {
602  Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
603  (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
604 
605  =
606 
607  Block<const MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
608  (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
609  }
610  }
611  }
612 
613  protected:
614  const PermutationType& m_permutation;
615  typename MatrixType::Nested m_matrix;
616 };
617 
618 /* Template partial specialization for transposed/inverse permutations */
619 
620 template<typename Derived>
621 struct traits<Transpose<PermutationBase<Derived> > >
622  : traits<Derived>
623 {};
624 
625 } // end namespace internal
626 
627 template<typename Derived>
628 class Transpose<PermutationBase<Derived> >
629  : public EigenBase<Transpose<PermutationBase<Derived> > >
630 {
631  typedef Derived PermutationType;
632  typedef typename PermutationType::IndicesType IndicesType;
633  typedef typename PermutationType::PlainPermutationType PlainPermutationType;
634  public:
635 
636  #ifndef EIGEN_PARSED_BY_DOXYGEN
637  typedef internal::traits<PermutationType> Traits;
638  typedef typename Derived::DenseMatrixType DenseMatrixType;
639  enum {
640  Flags = Traits::Flags,
641  CoeffReadCost = Traits::CoeffReadCost,
642  RowsAtCompileTime = Traits::RowsAtCompileTime,
643  ColsAtCompileTime = Traits::ColsAtCompileTime,
644  MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
645  MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
646  };
647  typedef typename Traits::Scalar Scalar;
648  #endif
649 
650  Transpose(const PermutationType& p) : m_permutation(p) {}
651 
652  inline int rows() const { return m_permutation.rows(); }
653  inline int cols() const { return m_permutation.cols(); }
654 
655  #ifndef EIGEN_PARSED_BY_DOXYGEN
656  template<typename DenseDerived>
657  void evalTo(MatrixBase<DenseDerived>& other) const
658  {
659  other.setZero();
660  for (int i=0; i<rows();++i)
661  other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1);
662  }
663  #endif
664 
666  PlainPermutationType eval() const { return *this; }
667 
668  DenseMatrixType toDenseMatrix() const { return *this; }
669 
672  template<typename OtherDerived> friend
673  inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>
674  operator*(const MatrixBase<OtherDerived>& matrix, const Transpose& trPerm)
675  {
676  return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>(trPerm.m_permutation, matrix.derived());
677  }
678 
681  template<typename OtherDerived>
682  inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>
683  operator*(const MatrixBase<OtherDerived>& matrix) const
684  {
685  return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>(m_permutation, matrix.derived());
686  }
687 
688  const PermutationType& nestedPermutation() const { return m_permutation; }
689 
690  protected:
691  const PermutationType& m_permutation;
692 };
693 
694 template<typename Derived>
696 {
697  return derived();
698 }
699 
700 } // end namespace Eigen
701 
702 #endif // EIGEN_PERMUTATIONMATRIX_H