Eigen  3.2.10
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 
258  Index determinant() const
259  {
260  Index res = 1;
261  Index n = size();
263  mask.fill(false);
264  Index r = 0;
265  while(r < n)
266  {
267  // search for the next seed
268  while(r<n && mask[r]) r++;
269  if(r>=n)
270  break;
271  // we got one, let's follow it until we are back to the seed
272  Index k0 = r++;
273  mask.coeffRef(k0) = true;
274  for(Index k=indices().coeff(k0); k!=k0; k=indices().coeff(k))
275  {
276  mask.coeffRef(k) = true;
277  res = -res;
278  }
279  }
280  return res;
281  }
282 
283  protected:
284 
285 };
286 
301 namespace internal {
302 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
303 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
304  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
305 {
306  typedef IndexType Index;
308 };
309 }
310 
311 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
312 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
313 {
315  typedef internal::traits<PermutationMatrix> Traits;
316  public:
317 
318  #ifndef EIGEN_PARSED_BY_DOXYGEN
319  typedef typename Traits::IndicesType IndicesType;
320  #endif
321 
322  inline PermutationMatrix()
323  {}
324 
327  inline PermutationMatrix(int size) : m_indices(size)
328  {}
329 
331  template<typename OtherDerived>
333  : m_indices(other.indices()) {}
334 
335  #ifndef EIGEN_PARSED_BY_DOXYGEN
336 
338  inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
339  #endif
340 
348  template<typename Other>
349  explicit inline PermutationMatrix(const MatrixBase<Other>& a_indices) : m_indices(a_indices)
350  {}
351 
353  template<typename Other>
354  explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
355  : m_indices(tr.size())
356  {
357  *this = tr;
358  }
359 
361  template<typename Other>
363  {
364  m_indices = other.indices();
365  return *this;
366  }
367 
369  template<typename Other>
370  PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
371  {
372  return Base::operator=(tr.derived());
373  }
374 
375  #ifndef EIGEN_PARSED_BY_DOXYGEN
376 
379  PermutationMatrix& operator=(const PermutationMatrix& other)
380  {
381  m_indices = other.m_indices;
382  return *this;
383  }
384  #endif
385 
387  const IndicesType& indices() const { return m_indices; }
389  IndicesType& indices() { return m_indices; }
390 
391 
392  /**** multiplication helpers to hopefully get RVO ****/
393 
394 #ifndef EIGEN_PARSED_BY_DOXYGEN
395  template<typename Other>
397  : m_indices(other.nestedPermutation().size())
398  {
399  for (int i=0; i<m_indices.size();++i) m_indices.coeffRef(other.nestedPermutation().indices().coeff(i)) = i;
400  }
401  template<typename Lhs,typename Rhs>
402  PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
403  : m_indices(lhs.indices().size())
404  {
405  Base::assignProduct(lhs,rhs);
406  }
407 #endif
408 
409  protected:
410 
411  IndicesType m_indices;
412 };
413 
414 
415 namespace internal {
416 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
417 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
418  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
419 {
420  typedef IndexType Index;
422 };
423 }
424 
425 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
426 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess>
427  : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
428 {
429  typedef PermutationBase<Map> Base;
430  typedef internal::traits<Map> Traits;
431  public:
432 
433  #ifndef EIGEN_PARSED_BY_DOXYGEN
434  typedef typename Traits::IndicesType IndicesType;
435  typedef typename IndicesType::Scalar Index;
436  #endif
437 
438  inline Map(const Index* indicesPtr)
439  : m_indices(indicesPtr)
440  {}
441 
442  inline Map(const Index* indicesPtr, Index size)
443  : m_indices(indicesPtr,size)
444  {}
445 
447  template<typename Other>
448  Map& operator=(const PermutationBase<Other>& other)
449  { return Base::operator=(other.derived()); }
450 
452  template<typename Other>
453  Map& operator=(const TranspositionsBase<Other>& tr)
454  { return Base::operator=(tr.derived()); }
455 
456  #ifndef EIGEN_PARSED_BY_DOXYGEN
457 
460  Map& operator=(const Map& other)
461  {
462  m_indices = other.m_indices;
463  return *this;
464  }
465  #endif
466 
468  const IndicesType& indices() const { return m_indices; }
470  IndicesType& indices() { return m_indices; }
471 
472  protected:
473 
474  IndicesType m_indices;
475 };
476 
489 struct PermutationStorage {};
490 
491 template<typename _IndicesType> class TranspositionsWrapper;
492 namespace internal {
493 template<typename _IndicesType>
494 struct traits<PermutationWrapper<_IndicesType> >
495 {
496  typedef PermutationStorage StorageKind;
497  typedef typename _IndicesType::Scalar Scalar;
498  typedef typename _IndicesType::Scalar Index;
499  typedef _IndicesType IndicesType;
500  enum {
501  RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
502  ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
503  MaxRowsAtCompileTime = IndicesType::MaxRowsAtCompileTime,
504  MaxColsAtCompileTime = IndicesType::MaxColsAtCompileTime,
505  Flags = 0,
506  CoeffReadCost = _IndicesType::CoeffReadCost
507  };
508 };
509 }
510 
511 template<typename _IndicesType>
512 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
513 {
515  typedef internal::traits<PermutationWrapper> Traits;
516  public:
517 
518  #ifndef EIGEN_PARSED_BY_DOXYGEN
519  typedef typename Traits::IndicesType IndicesType;
520  #endif
521 
522  inline PermutationWrapper(const IndicesType& a_indices)
523  : m_indices(a_indices)
524  {}
525 
527  const typename internal::remove_all<typename IndicesType::Nested>::type&
528  indices() const { return m_indices; }
529 
530  protected:
531 
532  typename IndicesType::Nested m_indices;
533 };
534 
537 template<typename Derived, typename PermutationDerived>
538 inline const internal::permut_matrix_product_retval<PermutationDerived, Derived, OnTheRight>
539 operator*(const MatrixBase<Derived>& matrix,
540  const PermutationBase<PermutationDerived> &permutation)
541 {
542  return internal::permut_matrix_product_retval
543  <PermutationDerived, Derived, OnTheRight>
544  (permutation.derived(), matrix.derived());
545 }
546 
549 template<typename Derived, typename PermutationDerived>
550 inline const internal::permut_matrix_product_retval
551  <PermutationDerived, Derived, OnTheLeft>
552 operator*(const PermutationBase<PermutationDerived> &permutation,
553  const MatrixBase<Derived>& matrix)
554 {
555  return internal::permut_matrix_product_retval
556  <PermutationDerived, Derived, OnTheLeft>
557  (permutation.derived(), matrix.derived());
558 }
559 
560 namespace internal {
561 
562 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
563 struct traits<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
564 {
565  typedef typename MatrixType::PlainObject ReturnType;
566 };
567 
568 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
569 struct permut_matrix_product_retval
570  : public ReturnByValue<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
571 {
572  typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
573  typedef typename MatrixType::Index Index;
574 
575  permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
576  : m_permutation(perm), m_matrix(matrix)
577  {}
578 
579  inline Index rows() const { return m_matrix.rows(); }
580  inline Index cols() const { return m_matrix.cols(); }
581 
582  template<typename Dest> inline void evalTo(Dest& dst) const
583  {
584  const Index n = Side==OnTheLeft ? rows() : cols();
585  // FIXME we need an is_same for expression that is not sensitive to constness. For instance
586  // is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true.
587  const typename Dest::Scalar *dst_data = internal::extract_data(dst);
588  if( is_same<MatrixTypeNestedCleaned,Dest>::value
589  && blas_traits<MatrixTypeNestedCleaned>::HasUsableDirectAccess
590  && blas_traits<Dest>::HasUsableDirectAccess
591  && dst_data!=0 && dst_data == extract_data(m_matrix))
592  {
593  // apply the permutation inplace
595  mask.fill(false);
596  Index r = 0;
597  while(r < m_permutation.size())
598  {
599  // search for the next seed
600  while(r<m_permutation.size() && mask[r]) r++;
601  if(r>=m_permutation.size())
602  break;
603  // we got one, let's follow it until we are back to the seed
604  Index k0 = r++;
605  Index kPrev = k0;
606  mask.coeffRef(k0) = true;
607  for(Index k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
608  {
611  (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
612 
613  mask.coeffRef(k) = true;
614  kPrev = k;
615  }
616  }
617  }
618  else
619  {
620  for(int i = 0; i < n; ++i)
621  {
623  (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
624 
625  =
626 
628  (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
629  }
630  }
631  }
632 
633  protected:
634  const PermutationType& m_permutation;
635  typename MatrixType::Nested m_matrix;
636 };
637 
638 /* Template partial specialization for transposed/inverse permutations */
639 
640 template<typename Derived>
641 struct traits<Transpose<PermutationBase<Derived> > >
642  : traits<Derived>
643 {};
644 
645 } // end namespace internal
646 
647 template<typename Derived>
648 class Transpose<PermutationBase<Derived> >
649  : public EigenBase<Transpose<PermutationBase<Derived> > >
650 {
651  typedef Derived PermutationType;
652  typedef typename PermutationType::IndicesType IndicesType;
653  typedef typename PermutationType::PlainPermutationType PlainPermutationType;
654  public:
655 
656  #ifndef EIGEN_PARSED_BY_DOXYGEN
657  typedef internal::traits<PermutationType> Traits;
658  typedef typename Derived::DenseMatrixType DenseMatrixType;
659  enum {
660  Flags = Traits::Flags,
661  CoeffReadCost = Traits::CoeffReadCost,
662  RowsAtCompileTime = Traits::RowsAtCompileTime,
663  ColsAtCompileTime = Traits::ColsAtCompileTime,
664  MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
665  MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
666  };
667  typedef typename Traits::Scalar Scalar;
668  #endif
669 
670  Transpose(const PermutationType& p) : m_permutation(p) {}
671 
672  inline int rows() const { return m_permutation.rows(); }
673  inline int cols() const { return m_permutation.cols(); }
674 
675  #ifndef EIGEN_PARSED_BY_DOXYGEN
676  template<typename DenseDerived>
677  void evalTo(MatrixBase<DenseDerived>& other) const
678  {
679  other.setZero();
680  for (int i=0; i<rows();++i)
681  other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1);
682  }
683  #endif
684 
686  PlainPermutationType eval() const { return *this; }
687 
688  DenseMatrixType toDenseMatrix() const { return *this; }
689 
692  template<typename OtherDerived> friend
693  inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>
694  operator*(const MatrixBase<OtherDerived>& matrix, const Transpose& trPerm)
695  {
696  return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>(trPerm.m_permutation, matrix.derived());
697  }
698 
701  template<typename OtherDerived>
702  inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>
703  operator*(const MatrixBase<OtherDerived>& matrix) const
704  {
705  return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>(m_permutation, matrix.derived());
706  }
707 
708  const PermutationType& nestedPermutation() const { return m_permutation; }
709 
710  protected:
711  const PermutationType& m_permutation;
712 };
713 
714 template<typename Derived>
716 {
717  return derived();
718 }
719 
720 } // end namespace Eigen
721 
722 #endif // EIGEN_PERMUTATIONMATRIX_H
void setIdentity()
Definition: PermutationMatrix.h:148
const IndicesType & indices() const
Definition: PermutationMatrix.h:387
Derived & applyTranspositionOnTheRight(Index i, Index j)
Definition: PermutationMatrix.h:190
Transpose< PermutationBase > inverse() const
Definition: PermutationMatrix.h:201
IndicesType & indices()
Definition: PermutationMatrix.h:138
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:104
void fill(const Scalar &value)
Definition: CwiseNullaryOp.h:322
Expression of the transpose of a matrix.
Definition: Transpose.h:57
Definition: LDLT.h:16
Definition: Constants.h:279
Derived & setZero()
Definition: CwiseNullaryOp.h:499
Index rows() const
Definition: PermutationMatrix.h:108
friend PlainPermutationType operator*(const Transpose< PermutationBase< Other > > &other, const PermutationBase &perm)
Definition: PermutationMatrix.h:251
PermutationMatrix(const MatrixBase< Other > &a_indices)
Definition: PermutationMatrix.h:349
Derived & applyTranspositionOnTheLeft(Index i, Index j)
Definition: PermutationMatrix.h:171
PermutationMatrix & operator=(const TranspositionsBase< Other > &tr)
Definition: PermutationMatrix.h:370
PermutationMatrix(const PermutationBase< OtherDerived > &other)
Definition: PermutationMatrix.h:332
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:312
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:362
Definition: Constants.h:277
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:327
Class to view a vector of integers as a permutation matrix.
Definition: PermutationMatrix.h:512
Index determinant() const
Definition: PermutationMatrix.h:258
Definition: Eigen_Colamd.h:50
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:102
PlainPermutationType operator*(const PermutationBase< Other > &other) const
Definition: PermutationMatrix.h:235
void setIdentity(Index newSize)
Definition: PermutationMatrix.h:156
IndicesType & indices()
Definition: PermutationMatrix.h:389
const internal::remove_all< typename IndicesType::Nested >::type & indices() const
Definition: PermutationMatrix.h:528
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:127
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:354