Eigen  3.2.92
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-2015 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 
40 namespace internal {
41 
42 enum PermPermProduct_t {PermPermProduct};
43 
44 } // end namespace internal
45 
46 template<typename Derived>
47 class PermutationBase : public EigenBase<Derived>
48 {
49  typedef internal::traits<Derived> Traits;
50  typedef EigenBase<Derived> Base;
51  public:
52 
53  #ifndef EIGEN_PARSED_BY_DOXYGEN
54  typedef typename Traits::IndicesType IndicesType;
55  enum {
56  Flags = Traits::Flags,
57  RowsAtCompileTime = Traits::RowsAtCompileTime,
58  ColsAtCompileTime = Traits::ColsAtCompileTime,
59  MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
60  MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
61  };
62  typedef typename Traits::StorageIndex StorageIndex;
64  DenseMatrixType;
66  PlainPermutationType;
67  typedef PlainPermutationType PlainObject;
68  using Base::derived;
69  typedef Inverse<Derived> InverseReturnType;
70  typedef void Scalar;
71  #endif
72 
74  template<typename OtherDerived>
76  {
77  indices() = other.indices();
78  return derived();
79  }
80 
82  template<typename OtherDerived>
83  Derived& operator=(const TranspositionsBase<OtherDerived>& tr)
84  {
85  setIdentity(tr.size());
86  for(Index k=size()-1; k>=0; --k)
87  applyTranspositionOnTheRight(k,tr.coeff(k));
88  return derived();
89  }
90 
91  #ifndef EIGEN_PARSED_BY_DOXYGEN
92 
95  Derived& operator=(const PermutationBase& other)
96  {
97  indices() = other.indices();
98  return derived();
99  }
100  #endif
101 
103  inline Index rows() const { return Index(indices().size()); }
104 
106  inline Index cols() const { return Index(indices().size()); }
107 
109  inline Index size() const { return Index(indices().size()); }
110 
111  #ifndef EIGEN_PARSED_BY_DOXYGEN
112  template<typename DenseDerived>
113  void evalTo(MatrixBase<DenseDerived>& other) const
114  {
115  other.setZero();
116  for (Index i=0; i<rows(); ++i)
117  other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
118  }
119  #endif
120 
125  DenseMatrixType toDenseMatrix() const
126  {
127  return derived();
128  }
129 
131  const IndicesType& indices() const { return derived().indices(); }
133  IndicesType& indices() { return derived().indices(); }
134 
137  inline void resize(Index newSize)
138  {
139  indices().resize(newSize);
140  }
141 
143  void setIdentity()
144  {
145  StorageIndex n = StorageIndex(size());
146  for(StorageIndex i = 0; i < n; ++i)
147  indices().coeffRef(i) = i;
148  }
149 
152  void setIdentity(Index newSize)
153  {
154  resize(newSize);
155  setIdentity();
156  }
157 
168  {
169  eigen_assert(i>=0 && j>=0 && i<size() && j<size());
170  for(Index k = 0; k < size(); ++k)
171  {
172  if(indices().coeff(k) == i) indices().coeffRef(k) = StorageIndex(j);
173  else if(indices().coeff(k) == j) indices().coeffRef(k) = StorageIndex(i);
174  }
175  return derived();
176  }
177 
187  {
188  eigen_assert(i>=0 && j>=0 && i<size() && j<size());
189  std::swap(indices().coeffRef(i), indices().coeffRef(j));
190  return derived();
191  }
192 
197  inline InverseReturnType inverse() const
198  { return InverseReturnType(derived()); }
203  inline InverseReturnType transpose() const
204  { return InverseReturnType(derived()); }
205 
206  /**** multiplication helpers to hopefully get RVO ****/
207 
208 
209 #ifndef EIGEN_PARSED_BY_DOXYGEN
210  protected:
211  template<typename OtherDerived>
212  void assignTranspose(const PermutationBase<OtherDerived>& other)
213  {
214  for (Index i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i;
215  }
216  template<typename Lhs,typename Rhs>
217  void assignProduct(const Lhs& lhs, const Rhs& rhs)
218  {
219  eigen_assert(lhs.cols() == rhs.rows());
220  for (Index i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
221  }
222 #endif
223 
224  public:
225 
230  template<typename Other>
231  inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
232  { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
233 
238  template<typename Other>
239  inline PlainPermutationType operator*(const InverseImpl<Other,PermutationStorage>& other) const
240  { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
241 
246  template<typename Other> friend
247  inline PlainPermutationType operator*(const InverseImpl<Other, PermutationStorage>& other, const PermutationBase& perm)
248  { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
249 
255  {
256  Index res = 1;
257  Index n = size();
259  mask.fill(false);
260  Index r = 0;
261  while(r < n)
262  {
263  // search for the next seed
264  while(r<n && mask[r]) r++;
265  if(r>=n)
266  break;
267  // we got one, let's follow it until we are back to the seed
268  Index k0 = r++;
269  mask.coeffRef(k0) = true;
270  for(Index k=indices().coeff(k0); k!=k0; k=indices().coeff(k))
271  {
272  mask.coeffRef(k) = true;
273  res = -res;
274  }
275  }
276  return res;
277  }
278 
279  protected:
280 
281 };
282 
297 namespace internal {
298 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex>
299 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> >
300  : traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
301 {
302  typedef PermutationStorage StorageKind;
304  typedef _StorageIndex StorageIndex;
305  typedef void Scalar;
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 
382  PermutationMatrix& operator=(const PermutationMatrix& other)
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>
399  PermutationMatrix(const InverseImpl<Other,PermutationStorage>& other)
400  : m_indices(other.derived().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.derived().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;
428  typedef _StorageIndex StorageIndex;
429  typedef void Scalar;
430 };
431 }
432 
433 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess>
434 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess>
435  : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> >
436 {
437  typedef PermutationBase<Map> Base;
438  typedef internal::traits<Map> Traits;
439  public:
440 
441  #ifndef EIGEN_PARSED_BY_DOXYGEN
442  typedef typename Traits::IndicesType IndicesType;
443  typedef typename IndicesType::Scalar StorageIndex;
444  #endif
445 
446  inline Map(const StorageIndex* indicesPtr)
447  : m_indices(indicesPtr)
448  {}
449 
450  inline Map(const StorageIndex* indicesPtr, Index size)
451  : m_indices(indicesPtr,size)
452  {}
453 
455  template<typename Other>
456  Map& operator=(const PermutationBase<Other>& other)
457  { return Base::operator=(other.derived()); }
458 
460  template<typename Other>
461  Map& operator=(const TranspositionsBase<Other>& tr)
462  { return Base::operator=(tr.derived()); }
463 
464  #ifndef EIGEN_PARSED_BY_DOXYGEN
465 
468  Map& operator=(const Map& other)
469  {
470  m_indices = other.m_indices;
471  return *this;
472  }
473  #endif
474 
476  const IndicesType& indices() const { return m_indices; }
478  IndicesType& indices() { return m_indices; }
479 
480  protected:
481 
482  IndicesType m_indices;
483 };
484 
497 template<typename _IndicesType> class TranspositionsWrapper;
498 namespace internal {
499 template<typename _IndicesType>
500 struct traits<PermutationWrapper<_IndicesType> >
501 {
502  typedef PermutationStorage StorageKind;
503  typedef void Scalar;
504  typedef typename _IndicesType::Scalar StorageIndex;
505  typedef _IndicesType IndicesType;
506  enum {
507  RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
508  ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
509  MaxRowsAtCompileTime = IndicesType::MaxSizeAtCompileTime,
510  MaxColsAtCompileTime = IndicesType::MaxSizeAtCompileTime,
511  Flags = 0
512  };
513 };
514 }
515 
516 template<typename _IndicesType>
517 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
518 {
520  typedef internal::traits<PermutationWrapper> Traits;
521  public:
522 
523  #ifndef EIGEN_PARSED_BY_DOXYGEN
524  typedef typename Traits::IndicesType IndicesType;
525  #endif
526 
527  inline PermutationWrapper(const IndicesType& indices)
528  : m_indices(indices)
529  {}
530 
532  const typename internal::remove_all<typename IndicesType::Nested>::type&
533  indices() const { return m_indices; }
534 
535  protected:
536 
537  typename IndicesType::Nested m_indices;
538 };
539 
540 
543 template<typename MatrixDerived, typename PermutationDerived>
544 EIGEN_DEVICE_FUNC
546 operator*(const MatrixBase<MatrixDerived> &matrix,
547  const PermutationBase<PermutationDerived>& permutation)
548 {
550  (matrix.derived(), permutation.derived());
551 }
552 
555 template<typename PermutationDerived, typename MatrixDerived>
556 EIGEN_DEVICE_FUNC
558 operator*(const PermutationBase<PermutationDerived> &permutation,
559  const MatrixBase<MatrixDerived>& matrix)
560 {
562  (permutation.derived(), matrix.derived());
563 }
564 
565 
566 template<typename PermutationType>
567 class InverseImpl<PermutationType, PermutationStorage>
568  : public EigenBase<Inverse<PermutationType> >
569 {
570  typedef typename PermutationType::PlainPermutationType PlainPermutationType;
571  typedef internal::traits<PermutationType> PermTraits;
572  protected:
573  InverseImpl() {}
574  public:
575  typedef Inverse<PermutationType> InverseType;
576  using EigenBase<Inverse<PermutationType> >::derived;
577 
578  #ifndef EIGEN_PARSED_BY_DOXYGEN
579  typedef typename PermutationType::DenseMatrixType DenseMatrixType;
580  enum {
581  RowsAtCompileTime = PermTraits::RowsAtCompileTime,
582  ColsAtCompileTime = PermTraits::ColsAtCompileTime,
583  MaxRowsAtCompileTime = PermTraits::MaxRowsAtCompileTime,
584  MaxColsAtCompileTime = PermTraits::MaxColsAtCompileTime
585  };
586  #endif
587 
588  #ifndef EIGEN_PARSED_BY_DOXYGEN
589  template<typename DenseDerived>
590  void evalTo(MatrixBase<DenseDerived>& other) const
591  {
592  other.setZero();
593  for (Index i=0; i<derived().rows();++i)
594  other.coeffRef(i, derived().nestedExpression().indices().coeff(i)) = typename DenseDerived::Scalar(1);
595  }
596  #endif
597 
599  PlainPermutationType eval() const { return derived(); }
600 
601  DenseMatrixType toDenseMatrix() const { return derived(); }
602 
605  template<typename OtherDerived> friend
607  operator*(const MatrixBase<OtherDerived>& matrix, const InverseType& trPerm)
608  {
609  return Product<OtherDerived, InverseType, AliasFreeProduct>(matrix.derived(), trPerm.derived());
610  }
611 
614  template<typename OtherDerived>
616  operator*(const MatrixBase<OtherDerived>& matrix) const
617  {
618  return Product<InverseType, OtherDerived, AliasFreeProduct>(derived(), matrix.derived());
619  }
620 };
621 
622 template<typename Derived>
624 {
625  return derived();
626 }
627 
628 namespace internal {
629 
630 template<> struct AssignmentKind<DenseShape,PermutationShape> { typedef EigenBase2EigenBase Kind; };
631 
632 } // end namespace internal
633 
634 } // end namespace Eigen
635 
636 #endif // EIGEN_PERMUTATIONMATRIX_H
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:107
PermutationMatrix(const MatrixBase< Other > &indices)
Definition: PermutationMatrix.h:352
friend PlainPermutationType operator*(const InverseImpl< Other, PermutationStorage > &other, const PermutationBase &perm)
Definition: PermutationMatrix.h:247
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:89
const internal::remove_all< typename IndicesType::Nested >::type & indices() const
Definition: PermutationMatrix.h:533
Definition: Constants.h:499
Derived & applyTranspositionOnTheRight(Index i, Index j)
Definition: PermutationMatrix.h:186
Index size() const
Definition: PermutationMatrix.h:109
Definition: LDLT.h:16
const IndicesType & indices() const
Definition: PermutationMatrix.h:131
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
Derived & derived()
Definition: EigenBase.h:44
Derived & operator=(const PermutationBase< OtherDerived > &other)
Definition: PermutationMatrix.h:75
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:37
PermutationMatrix & operator=(const PermutationBase< Other > &other)
Definition: PermutationMatrix.h:365
Index rows() const
Definition: PermutationMatrix.h:103
Base class for permutations.
Definition: PermutationMatrix.h:47
Definition: EigenBase.h:28
Expression of the inverse of another expression.
Definition: Inverse.h:43
Permutation matrix.
Definition: PermutationMatrix.h:310
PermutationMatrix & operator=(const TranspositionsBase< Other > &tr)
Definition: PermutationMatrix.h:373
IndicesType & indices()
Definition: PermutationMatrix.h:133
PlainPermutationType operator*(const PermutationBase< Other > &other) const
Definition: PermutationMatrix.h:231
void fill(const Scalar &value)
Definition: CwiseNullaryOp.h:326
Derived & operator=(const TranspositionsBase< OtherDerived > &tr)
Definition: PermutationMatrix.h:83
DenseMatrixType toDenseMatrix() const
Definition: PermutationMatrix.h:125
Index determinant() const
Definition: PermutationMatrix.h:254
Class to view a vector of integers as a permutation matrix.
Definition: PermutationMatrix.h:517
InverseReturnType transpose() const
Definition: PermutationMatrix.h:203
void setIdentity()
Definition: PermutationMatrix.h:143
PermutationMatrix(const TranspositionsBase< Other > &tr)
Definition: PermutationMatrix.h:357
InverseReturnType inverse() const
Definition: PermutationMatrix.h:197
Definition: Eigen_Colamd.h:54
Derived & setZero()
Definition: CwiseNullaryOp.h:504
Derived & applyTranspositionOnTheLeft(Index i, Index j)
Definition: PermutationMatrix.h:167
PermutationMatrix(Index size)
Definition: PermutationMatrix.h:328
const IndicesType & indices() const
Definition: PermutationMatrix.h:390
IndicesType & indices()
Definition: PermutationMatrix.h:392
void setIdentity(Index newSize)
Definition: PermutationMatrix.h:152
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
void resize(Index newSize)
Definition: PermutationMatrix.h:137
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Index cols() const
Definition: PermutationMatrix.h:106
PermutationMatrix(const PermutationBase< OtherDerived > &other)
Definition: PermutationMatrix.h:335
PlainPermutationType operator*(const InverseImpl< Other, PermutationStorage > &other) const
Definition: PermutationMatrix.h:239