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 PermutationType, typename MatrixType, int Side, bool Transposed>
20 struct traits<permut_sparsematrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
21 {
22   typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
23   typedef typename MatrixTypeNestedCleaned::Scalar Scalar;
24   typedef typename MatrixTypeNestedCleaned::Index Index;
25   enum {
26     SrcStorageOrder = MatrixTypeNestedCleaned::Flags&RowMajorBit ? RowMajor : ColMajor,
27     MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight
28   };
29 
30   typedef typename internal::conditional<MoveOuter,
31         SparseMatrix<Scalar,SrcStorageOrder,Index>,
32         SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,Index> >::type ReturnType;
33 };
34 
35 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
36 struct permut_sparsematrix_product_retval
37  : public ReturnByValue<permut_sparsematrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
38 {
39     typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
40     typedef typename MatrixTypeNestedCleaned::Scalar Scalar;
41     typedef typename MatrixTypeNestedCleaned::Index Index;
42 
43     enum {
44       SrcStorageOrder = MatrixTypeNestedCleaned::Flags&RowMajorBit ? RowMajor : ColMajor,
45       MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight
46     };
47 
48     permut_sparsematrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
49       : m_permutation(perm), m_matrix(matrix)
50     {}
51 
52     inline int rows() const { return m_matrix.rows(); }
53     inline int cols() const { return m_matrix.cols(); }
54 
55     template<typename Dest> inline void evalTo(Dest& dst) const
56     {
57       if(MoveOuter)
58       {
59         SparseMatrix<Scalar,SrcStorageOrder,Index> tmp(m_matrix.rows(), m_matrix.cols());
60         Matrix<Index,Dynamic,1> sizes(m_matrix.outerSize());
61         for(Index j=0; j<m_matrix.outerSize(); ++j)
62         {
63           Index jp = m_permutation.indices().coeff(j);
64           sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = m_matrix.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).nonZeros();
65         }
66         tmp.reserve(sizes);
67         for(Index j=0; j<m_matrix.outerSize(); ++j)
68         {
69           Index jp = m_permutation.indices().coeff(j);
70           Index jsrc = ((Side==OnTheRight) ^ Transposed) ? jp : j;
71           Index jdst = ((Side==OnTheLeft) ^ Transposed) ? jp : j;
72           for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,jsrc); it; ++it)
73             tmp.insertByOuterInner(jdst,it.index()) = it.value();
74         }
75         dst = tmp;
76       }
77       else
78       {
79         SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,Index> tmp(m_matrix.rows(), m_matrix.cols());
80         Matrix<Index,Dynamic,1> sizes(tmp.outerSize());
81         sizes.setZero();
82         PermutationMatrix<Dynamic,Dynamic,Index> perm;
83         if((Side==OnTheLeft) ^ Transposed)
84           perm = m_permutation;
85         else
86           perm = m_permutation.transpose();
87 
88         for(Index j=0; j<m_matrix.outerSize(); ++j)
89           for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,j); it; ++it)
90             sizes[perm.indices().coeff(it.index())]++;
91         tmp.reserve(sizes);
92         for(Index j=0; j<m_matrix.outerSize(); ++j)
93           for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,j); it; ++it)
94             tmp.insertByOuterInner(perm.indices().coeff(it.index()),j) = it.value();
95         dst = tmp;
96       }
97     }
98 
99   protected:
100     const PermutationType& m_permutation;
101     typename MatrixType::Nested m_matrix;
102 };
103 
104 }
105 
106 
107 
108 /** \returns the matrix with the permutation applied to the columns
109   */
110 template<typename SparseDerived, typename PermDerived>
111 inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, false>
112 operator*(const SparseMatrixBase<SparseDerived>& matrix, const PermutationBase<PermDerived>& perm)
113 {
114   return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, false>(perm, matrix.derived());
115 }
116 
117 /** \returns the matrix with the permutation applied to the rows
118   */
119 template<typename SparseDerived, typename PermDerived>
120 inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, false>
121 operator*( const PermutationBase<PermDerived>& perm, const SparseMatrixBase<SparseDerived>& matrix)
122 {
123   return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, false>(perm, matrix.derived());
124 }
125 
126 
127 
128 /** \returns the matrix with the inverse permutation applied to the columns.
129   */
130 template<typename SparseDerived, typename PermDerived>
131 inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, true>
132 operator*(const SparseMatrixBase<SparseDerived>& matrix, const Transpose<PermutationBase<PermDerived> >& tperm)
133 {
134   return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, true>(tperm.nestedPermutation(), matrix.derived());
135 }
136 
137 /** \returns the matrix with the inverse permutation applied to the rows.
138   */
139 template<typename SparseDerived, typename PermDerived>
140 inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, true>
141 operator*(const Transpose<PermutationBase<PermDerived> >& tperm, const SparseMatrixBase<SparseDerived>& matrix)
142 {
143   return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, true>(tperm.nestedPermutation(), matrix.derived());
144 }
145 
146 } // end namespace Eigen
147 
148 #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H
149