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 ExpressionType, int Side, bool Transposed> 20 struct 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 86 namespace internal { 87 88 template <int ProductTag> struct product_promote_storage_type<Sparse, PermutationStorage, ProductTag> { typedef Sparse ret; }; 89 template <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 95 template<typename Lhs, typename Rhs, int ProductTag> 96 struct 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 114 protected: 115 PlainObject m_result; 116 }; 117 118 template<typename Lhs, typename Rhs, int ProductTag> 119 struct 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 137 protected: 138 PlainObject m_result; 139 }; 140 141 } // end namespace internal 142 143 /** \returns the matrix with the permutation applied to the columns 144 */ 145 template<typename SparseDerived, typename PermDerived> 146 inline const Product<SparseDerived, PermDerived, AliasFreeProduct> 147 operator*(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 */ 152 template<typename SparseDerived, typename PermDerived> 153 inline const Product<PermDerived, SparseDerived, AliasFreeProduct> 154 operator*( 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 */ 160 template<typename SparseDerived, typename PermutationType> 161 inline const Product<SparseDerived, Inverse<PermutationType>, AliasFreeProduct> 162 operator*(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 */ 169 template<typename SparseDerived, typename PermutationType> 170 inline const Product<Inverse<PermutationType>, SparseDerived, AliasFreeProduct> 171 operator*(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