Eigen  3.2.91
SparsePermutation.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_SPARSE_PERMUTATION_H
11 #define EIGEN_SPARSE_PERMUTATION_H
12 
13 // This file implements sparse * permutation products
14 
15 namespace Eigen {
16 
17 namespace internal {
18 
19 template<typename MatrixType, int Side, bool Transposed>
20 struct permutation_matrix_product<MatrixType, Side, Transposed, SparseShape>
21 {
22  typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
23  typedef typename MatrixTypeNestedCleaned::Scalar Scalar;
24  typedef typename MatrixTypeNestedCleaned::StorageIndex StorageIndex;
25 
26  enum {
27  SrcStorageOrder = MatrixTypeNestedCleaned::Flags&RowMajorBit ? RowMajor : ColMajor,
28  MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight
29  };
30 
31  typedef typename internal::conditional<MoveOuter,
32  SparseMatrix<Scalar,SrcStorageOrder,StorageIndex>,
33  SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,StorageIndex> >::type ReturnType;
34 
35  template<typename Dest,typename PermutationType>
36  static inline void run(Dest& dst, const PermutationType& perm, const MatrixType& mat)
37  {
38  if(MoveOuter)
39  {
40  SparseMatrix<Scalar,SrcStorageOrder,StorageIndex> tmp(mat.rows(), mat.cols());
41  Matrix<StorageIndex,Dynamic,1> sizes(mat.outerSize());
42  for(Index j=0; j<mat.outerSize(); ++j)
43  {
44  Index jp = perm.indices().coeff(j);
45  sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = StorageIndex(mat.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).nonZeros());
46  }
47  tmp.reserve(sizes);
48  for(Index j=0; j<mat.outerSize(); ++j)
49  {
50  Index jp = perm.indices().coeff(j);
51  Index jsrc = ((Side==OnTheRight) ^ Transposed) ? jp : j;
52  Index jdst = ((Side==OnTheLeft) ^ Transposed) ? jp : j;
53  for(typename MatrixTypeNestedCleaned::InnerIterator it(mat,jsrc); it; ++it)
54  tmp.insertByOuterInner(jdst,it.index()) = it.value();
55  }
56  dst = tmp;
57  }
58  else
59  {
60  SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,StorageIndex> tmp(mat.rows(), mat.cols());
61  Matrix<StorageIndex,Dynamic,1> sizes(tmp.outerSize());
62  sizes.setZero();
63  PermutationMatrix<Dynamic,Dynamic,StorageIndex> perm_cpy;
64  if((Side==OnTheLeft) ^ Transposed)
65  perm_cpy = perm;
66  else
67  perm_cpy = perm.transpose();
68 
69  for(Index j=0; j<mat.outerSize(); ++j)
70  for(typename MatrixTypeNestedCleaned::InnerIterator it(mat,j); it; ++it)
71  sizes[perm_cpy.indices().coeff(it.index())]++;
72  tmp.reserve(sizes);
73  for(Index j=0; j<mat.outerSize(); ++j)
74  for(typename MatrixTypeNestedCleaned::InnerIterator it(mat,j); it; ++it)
75  tmp.insertByOuterInner(perm_cpy.indices().coeff(it.index()),j) = it.value();
76  dst = tmp;
77  }
78  }
79 };
80 
81 }
82 
83 namespace internal {
84 
85 template <int ProductTag> struct product_promote_storage_type<Sparse, PermutationStorage, ProductTag> { typedef Sparse ret; };
86 template <int ProductTag> struct product_promote_storage_type<PermutationStorage, Sparse, ProductTag> { typedef Sparse ret; };
87 
88 // TODO, the following two overloads are only needed to define the right temporary type through
89 // typename traits<permutation_sparse_matrix_product<Rhs,Lhs,OnTheRight,false> >::ReturnType
90 // whereas it should be correctly handled by traits<Product<> >::PlainObject
91 
92 template<typename Lhs, typename Rhs, int ProductTag>
93 struct product_evaluator<Product<Lhs, Rhs, AliasFreeProduct>, ProductTag, PermutationShape, SparseShape, typename traits<Lhs>::Scalar, typename traits<Rhs>::Scalar>
94  : public evaluator<typename permutation_matrix_product<Rhs,OnTheRight,false,SparseShape>::ReturnType>
95 {
96  typedef Product<Lhs, Rhs, AliasFreeProduct> XprType;
97  typedef typename permutation_matrix_product<Rhs,OnTheRight,false,SparseShape>::ReturnType PlainObject;
98  typedef evaluator<PlainObject> Base;
99 
100  explicit product_evaluator(const XprType& xpr)
101  : m_result(xpr.rows(), xpr.cols())
102  {
103  ::new (static_cast<Base*>(this)) Base(m_result);
104  generic_product_impl<Lhs, Rhs, PermutationShape, SparseShape, ProductTag>::evalTo(m_result, xpr.lhs(), xpr.rhs());
105  }
106 
107 protected:
108  PlainObject m_result;
109 };
110 
111 template<typename Lhs, typename Rhs, int ProductTag>
112 struct product_evaluator<Product<Lhs, Rhs, AliasFreeProduct>, ProductTag, SparseShape, PermutationShape, typename traits<Lhs>::Scalar, typename traits<Rhs>::Scalar>
113  : public evaluator<typename permutation_matrix_product<Lhs,OnTheRight,false,SparseShape>::ReturnType>
114 {
115  typedef Product<Lhs, Rhs, AliasFreeProduct> XprType;
116  typedef typename permutation_matrix_product<Lhs,OnTheRight,false,SparseShape>::ReturnType PlainObject;
117  typedef evaluator<PlainObject> Base;
118 
119  explicit product_evaluator(const XprType& xpr)
120  : m_result(xpr.rows(), xpr.cols())
121  {
122  ::new (static_cast<Base*>(this)) Base(m_result);
123  generic_product_impl<Lhs, Rhs, SparseShape, PermutationShape, ProductTag>::evalTo(m_result, xpr.lhs(), xpr.rhs());
124  }
125 
126 protected:
127  PlainObject m_result;
128 };
129 
130 } // end namespace internal
131 
134 template<typename SparseDerived, typename PermDerived>
135 inline const Product<SparseDerived, PermDerived>
136 operator*(const SparseMatrixBase<SparseDerived>& matrix, const PermutationBase<PermDerived>& perm)
137 { return Product<SparseDerived, PermDerived>(matrix.derived(), perm.derived()); }
138 
141 template<typename SparseDerived, typename PermDerived>
142 inline const Product<PermDerived, SparseDerived>
143 operator*( const PermutationBase<PermDerived>& perm, const SparseMatrixBase<SparseDerived>& matrix)
144 { return Product<PermDerived, SparseDerived>(perm.derived(), matrix.derived()); }
145 
146 
147 // TODO, the following specializations should not be needed as Transpose<Permutation*> should be a PermutationBase.
150 template<typename SparseDerived, typename PermDerived>
151 inline const Product<SparseDerived, Transpose<PermutationBase<PermDerived> > >
152 operator*(const SparseMatrixBase<SparseDerived>& matrix, const Transpose<PermutationBase<PermDerived> >& tperm)
153 {
154  return Product<SparseDerived, Transpose<PermutationBase<PermDerived> > >(matrix.derived(), tperm);
155 }
156 
159 template<typename SparseDerived, typename PermDerived>
160 inline const Product<Transpose<PermutationBase<PermDerived> >, SparseDerived>
161 operator*(const Transpose<PermutationBase<PermDerived> >& tperm, const SparseMatrixBase<SparseDerived>& matrix)
162 {
163  return Product<Transpose<PermutationBase<PermDerived> >, SparseDerived>(tperm, matrix.derived());
164 }
165 
166 } // end namespace Eigen
167 
168 #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H
Definition: Constants.h:314
Definition: LDLT.h:16
Definition: Constants.h:327
const unsigned int RowMajorBit
Definition: Constants.h:53
Definition: Eigen_Colamd.h:54
Definition: Constants.h:312
Definition: Constants.h:325