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 PermutationType, typename MatrixType, int Side, bool Transposed> 20struct 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 35template<typename PermutationType, typename MatrixType, int Side, bool Transposed> 36struct 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 */ 110template<typename SparseDerived, typename PermDerived> 111inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, false> 112operator*(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 */ 119template<typename SparseDerived, typename PermDerived> 120inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, false> 121operator*( 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 */ 130template<typename SparseDerived, typename PermDerived> 131inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, true> 132operator*(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 */ 139template<typename SparseDerived, typename PermDerived> 140inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, true> 141operator*(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