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
15namespace Eigen {
16
17namespace internal {
18
19template<typename ExpressionType, int Side, bool Transposed>
20struct permutation_matrix_product<ExpressionType, Side, Transposed, SparseShape>
21{
22 typedef typename nested_eval<ExpressionType, 1>::type MatrixType;
23 typedef typename remove_all<MatrixType>::type MatrixTypeCleaned;
24
25 typedef typename MatrixTypeCleaned::Scalar Scalar;
26 typedef typename MatrixTypeCleaned::StorageIndex StorageIndex;
27
28 enum {
29 SrcStorageOrder = MatrixTypeCleaned::Flags&RowMajorBit ? RowMajor : ColMajor,
30 MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight
31 };
32
33 typedef typename internal::conditional<MoveOuter,
34 SparseMatrix<Scalar,SrcStorageOrder,StorageIndex>,
35 SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,StorageIndex> >::type ReturnType;
36
37 template<typename Dest,typename PermutationType>
38 static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr)
39 {
40 MatrixType mat(xpr);
41 if(MoveOuter)
42 {
43 SparseMatrix<Scalar,SrcStorageOrder,StorageIndex> tmp(mat.rows(), mat.cols());
44 Matrix<StorageIndex,Dynamic,1> sizes(mat.outerSize());
45 for(Index j=0; j<mat.outerSize(); ++j)
46 {
47 Index jp = perm.indices().coeff(j);
48 sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = StorageIndex(mat.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).nonZeros());
49 }
50 tmp.reserve(sizes);
51 for(Index j=0; j<mat.outerSize(); ++j)
52 {
53 Index jp = perm.indices().coeff(j);
54 Index jsrc = ((Side==OnTheRight) ^ Transposed) ? jp : j;
55 Index jdst = ((Side==OnTheLeft) ^ Transposed) ? jp : j;
56 for(typename MatrixTypeCleaned::InnerIterator it(mat,jsrc); it; ++it)
57 tmp.insertByOuterInner(jdst,it.index()) = it.value();
58 }
59 dst = tmp;
60 }
61 else
62 {
63 SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,StorageIndex> tmp(mat.rows(), mat.cols());
64 Matrix<StorageIndex,Dynamic,1> sizes(tmp.outerSize());
65 sizes.setZero();
66 PermutationMatrix<Dynamic,Dynamic,StorageIndex> perm_cpy;
67 if((Side==OnTheLeft) ^ Transposed)
68 perm_cpy = perm;
69 else
70 perm_cpy = perm.transpose();
71
72 for(Index j=0; j<mat.outerSize(); ++j)
73 for(typename MatrixTypeCleaned::InnerIterator it(mat,j); it; ++it)
74 sizes[perm_cpy.indices().coeff(it.index())]++;
75 tmp.reserve(sizes);
76 for(Index j=0; j<mat.outerSize(); ++j)
77 for(typename MatrixTypeCleaned::InnerIterator it(mat,j); it; ++it)
78 tmp.insertByOuterInner(perm_cpy.indices().coeff(it.index()),j) = it.value();
79 dst = tmp;
80 }
81 }
82};
83
84}
85
86namespace internal {
87
88template <int ProductTag> struct product_promote_storage_type<Sparse, PermutationStorage, ProductTag> { typedef Sparse ret; };
89template <int ProductTag> struct product_promote_storage_type<PermutationStorage, Sparse, ProductTag> { typedef Sparse ret; };
90
91// TODO, the following two overloads are only needed to define the right temporary type through
92// typename traits<permutation_sparse_matrix_product<Rhs,Lhs,OnTheRight,false> >::ReturnType
93// whereas it should be correctly handled by traits<Product<> >::PlainObject
94
95template<typename Lhs, typename Rhs, int ProductTag>
96struct product_evaluator<Product<Lhs, Rhs, AliasFreeProduct>, ProductTag, PermutationShape, SparseShape>
97 : public evaluator<typename permutation_matrix_product<Rhs,OnTheLeft,false,SparseShape>::ReturnType>
98{
99 typedef Product<Lhs, Rhs, AliasFreeProduct> XprType;
100 typedef typename permutation_matrix_product<Rhs,OnTheLeft,false,SparseShape>::ReturnType PlainObject;
101 typedef evaluator<PlainObject> Base;
102
103 enum {
104 Flags = Base::Flags | EvalBeforeNestingBit
105 };
106
107 explicit product_evaluator(const XprType& xpr)
108 : m_result(xpr.rows(), xpr.cols())
109 {
110 ::new (static_cast<Base*>(this)) Base(m_result);
111 generic_product_impl<Lhs, Rhs, PermutationShape, SparseShape, ProductTag>::evalTo(m_result, xpr.lhs(), xpr.rhs());
112 }
113
114protected:
115 PlainObject m_result;
116};
117
118template<typename Lhs, typename Rhs, int ProductTag>
119struct product_evaluator<Product<Lhs, Rhs, AliasFreeProduct>, ProductTag, SparseShape, PermutationShape >
120 : public evaluator<typename permutation_matrix_product<Lhs,OnTheRight,false,SparseShape>::ReturnType>
121{
122 typedef Product<Lhs, Rhs, AliasFreeProduct> XprType;
123 typedef typename permutation_matrix_product<Lhs,OnTheRight,false,SparseShape>::ReturnType PlainObject;
124 typedef evaluator<PlainObject> Base;
125
126 enum {
127 Flags = Base::Flags | EvalBeforeNestingBit
128 };
129
130 explicit product_evaluator(const XprType& xpr)
131 : m_result(xpr.rows(), xpr.cols())
132 {
133 ::new (static_cast<Base*>(this)) Base(m_result);
134 generic_product_impl<Lhs, Rhs, SparseShape, PermutationShape, ProductTag>::evalTo(m_result, xpr.lhs(), xpr.rhs());
135 }
136
137protected:
138 PlainObject m_result;
139};
140
141} // end namespace internal
142
143/** \returns the matrix with the permutation applied to the columns
144 */
145template<typename SparseDerived, typename PermDerived>
146inline const Product<SparseDerived, PermDerived, AliasFreeProduct>
147operator*(const SparseMatrixBase<SparseDerived>& matrix, const PermutationBase<PermDerived>& perm)
148{ return Product<SparseDerived, PermDerived, AliasFreeProduct>(matrix.derived(), perm.derived()); }
149
150/** \returns the matrix with the permutation applied to the rows
151 */
152template<typename SparseDerived, typename PermDerived>
153inline const Product<PermDerived, SparseDerived, AliasFreeProduct>
154operator*( const PermutationBase<PermDerived>& perm, const SparseMatrixBase<SparseDerived>& matrix)
155{ return Product<PermDerived, SparseDerived, AliasFreeProduct>(perm.derived(), matrix.derived()); }
156
157
158/** \returns the matrix with the inverse permutation applied to the columns.
159 */
160template<typename SparseDerived, typename PermutationType>
161inline const Product<SparseDerived, Inverse<PermutationType>, AliasFreeProduct>
162operator*(const SparseMatrixBase<SparseDerived>& matrix, const InverseImpl<PermutationType, PermutationStorage>& tperm)
163{
164 return Product<SparseDerived, Inverse<PermutationType>, AliasFreeProduct>(matrix.derived(), tperm.derived());
165}
166
167/** \returns the matrix with the inverse permutation applied to the rows.
168 */
169template<typename SparseDerived, typename PermutationType>
170inline const Product<Inverse<PermutationType>, SparseDerived, AliasFreeProduct>
171operator*(const InverseImpl<PermutationType,PermutationStorage>& tperm, const SparseMatrixBase<SparseDerived>& matrix)
172{
173 return Product<Inverse<PermutationType>, SparseDerived, AliasFreeProduct>(tperm.derived(), matrix.derived());
174}
175
176} // end namespace Eigen
177
178#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H
179