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 (at) 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 VectorXi 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).size(); 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 VectorXi 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
