SparsePermutation.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) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_SPARSE_PERMUTATION_H
26 #define EIGEN_SPARSE_PERMUTATION_H
27 
28 // This file implements sparse * permutation products
29 
30 namespace Eigen {
31 
32 namespace internal {
33 
34 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
35 struct traits<permut_sparsematrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
36 {
37  typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
38  typedef typename MatrixTypeNestedCleaned::Scalar Scalar;
39  typedef typename MatrixTypeNestedCleaned::Index Index;
40  enum {
41  SrcStorageOrder = MatrixTypeNestedCleaned::Flags&RowMajorBit ? RowMajor : ColMajor,
42  MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight
43  };
44 
45  typedef typename internal::conditional<MoveOuter,
46  SparseMatrix<Scalar,SrcStorageOrder,Index>,
47  SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,Index> >::type ReturnType;
48 };
49 
50 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
51 struct permut_sparsematrix_product_retval
52  : public ReturnByValue<permut_sparsematrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
53 {
54  typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
55  typedef typename MatrixTypeNestedCleaned::Scalar Scalar;
56  typedef typename MatrixTypeNestedCleaned::Index Index;
57 
58  enum {
59  SrcStorageOrder = MatrixTypeNestedCleaned::Flags&RowMajorBit ? RowMajor : ColMajor,
60  MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight
61  };
62 
63  permut_sparsematrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
64  : m_permutation(perm), m_matrix(matrix)
65  {}
66 
67  inline int rows() const { return m_matrix.rows(); }
68  inline int cols() const { return m_matrix.cols(); }
69 
70  template<typename Dest> inline void evalTo(Dest& dst) const
71  {
72  if(MoveOuter)
73  {
74  SparseMatrix<Scalar,SrcStorageOrder,Index> tmp(m_matrix.rows(), m_matrix.cols());
75  VectorXi sizes(m_matrix.outerSize());
76  for(Index j=0; j<m_matrix.outerSize(); ++j)
77  {
78  Index jp = m_permutation.indices().coeff(j);
79  sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = m_matrix.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).size();
80  }
81  tmp.reserve(sizes);
82  for(Index j=0; j<m_matrix.outerSize(); ++j)
83  {
84  Index jp = m_permutation.indices().coeff(j);
85  Index jsrc = ((Side==OnTheRight) ^ Transposed) ? jp : j;
86  Index jdst = ((Side==OnTheLeft) ^ Transposed) ? jp : j;
87  for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,jsrc); it; ++it)
88  tmp.insertByOuterInner(jdst,it.index()) = it.value();
89  }
90  dst = tmp;
91  }
92  else
93  {
94  SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,Index> tmp(m_matrix.rows(), m_matrix.cols());
95  VectorXi sizes(tmp.outerSize());
96  sizes.setZero();
97  PermutationMatrix<Dynamic,Dynamic,Index> perm;
98  if((Side==OnTheLeft) ^ Transposed)
99  perm = m_permutation;
100  else
101  perm = m_permutation.transpose();
102 
103  for(Index j=0; j<m_matrix.outerSize(); ++j)
104  for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,j); it; ++it)
105  sizes[perm.indices().coeff(it.index())]++;
106  tmp.reserve(sizes);
107  for(Index j=0; j<m_matrix.outerSize(); ++j)
108  for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,j); it; ++it)
109  tmp.insertByOuterInner(perm.indices().coeff(it.index()),j) = it.value();
110  dst = tmp;
111  }
112  }
113 
114  protected:
115  const PermutationType& m_permutation;
116  typename MatrixType::Nested m_matrix;
117 };
118 
119 }
120 
121 
122 
125 template<typename SparseDerived, typename PermDerived>
126 inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, false>
128 {
129  return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, false>(perm, matrix.derived());
130 }
131 
134 template<typename SparseDerived, typename PermDerived>
135 inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, false>
137 {
138  return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, false>(perm, matrix.derived());
139 }
140 
141 
142 
145 template<typename SparseDerived, typename PermDerived>
146 inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, true>
148 {
149  return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, true>(tperm.nestedPermutation(), matrix.derived());
150 }
151 
154 template<typename SparseDerived, typename PermDerived>
155 inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, true>
157 {
158  return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, true>(tperm.nestedPermutation(), matrix.derived());
159 }
160 
161 } // end namespace Eigen
162 
163 #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H