![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_SPARSE_PERMUTATION_H 00011 #define EIGEN_SPARSE_PERMUTATION_H 00012 00013 // This file implements sparse * permutation products 00014 00015 namespace Eigen { 00016 00017 namespace internal { 00018 00019 template<typename ExpressionType, int Side, bool Transposed> 00020 struct permutation_matrix_product<ExpressionType, Side, Transposed, SparseShape> 00021 { 00022 typedef typename nested_eval<ExpressionType, 1>::type MatrixType; 00023 typedef typename remove_all<MatrixType>::type MatrixTypeCleaned; 00024 00025 typedef typename MatrixTypeCleaned::Scalar Scalar; 00026 typedef typename MatrixTypeCleaned::StorageIndex StorageIndex; 00027 00028 enum { 00029 SrcStorageOrder = MatrixTypeCleaned::Flags&RowMajorBit ? RowMajor : ColMajor, 00030 MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight 00031 }; 00032 00033 typedef typename internal::conditional<MoveOuter, 00034 SparseMatrix<Scalar,SrcStorageOrder,StorageIndex>, 00035 SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,StorageIndex> >::type ReturnType; 00036 00037 template<typename Dest,typename PermutationType> 00038 static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) 00039 { 00040 MatrixType mat(xpr); 00041 if(MoveOuter) 00042 { 00043 SparseMatrix<Scalar,SrcStorageOrder,StorageIndex> tmp(mat.rows(), mat.cols()); 00044 Matrix<StorageIndex,Dynamic,1> sizes(mat.outerSize()); 00045 for(Index j=0; j<mat.outerSize(); ++j) 00046 { 00047 Index jp = perm.indices().coeff(j); 00048 sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = StorageIndex(mat.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).nonZeros()); 00049 } 00050 tmp.reserve(sizes); 00051 for(Index j=0; j<mat.outerSize(); ++j) 00052 { 00053 Index jp = perm.indices().coeff(j); 00054 Index jsrc = ((Side==OnTheRight) ^ Transposed) ? jp : j; 00055 Index jdst = ((Side==OnTheLeft) ^ Transposed) ? jp : j; 00056 for(typename MatrixTypeCleaned::InnerIterator it(mat,jsrc); it; ++it) 00057 tmp.insertByOuterInner(jdst,it.index()) = it.value(); 00058 } 00059 dst = tmp; 00060 } 00061 else 00062 { 00063 SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,StorageIndex> tmp(mat.rows(), mat.cols()); 00064 Matrix<StorageIndex,Dynamic,1> sizes(tmp.outerSize()); 00065 sizes.setZero(); 00066 PermutationMatrix<Dynamic,Dynamic,StorageIndex> perm_cpy; 00067 if((Side==OnTheLeft) ^ Transposed) 00068 perm_cpy = perm; 00069 else 00070 perm_cpy = perm.transpose(); 00071 00072 for(Index j=0; j<mat.outerSize(); ++j) 00073 for(typename MatrixTypeCleaned::InnerIterator it(mat,j); it; ++it) 00074 sizes[perm_cpy.indices().coeff(it.index())]++; 00075 tmp.reserve(sizes); 00076 for(Index j=0; j<mat.outerSize(); ++j) 00077 for(typename MatrixTypeCleaned::InnerIterator it(mat,j); it; ++it) 00078 tmp.insertByOuterInner(perm_cpy.indices().coeff(it.index()),j) = it.value(); 00079 dst = tmp; 00080 } 00081 } 00082 }; 00083 00084 } 00085 00086 namespace internal { 00087 00088 template <int ProductTag> struct product_promote_storage_type<Sparse, PermutationStorage, ProductTag> { typedef Sparse ret; }; 00089 template <int ProductTag> struct product_promote_storage_type<PermutationStorage, Sparse, ProductTag> { typedef Sparse ret; }; 00090 00091 // TODO, the following two overloads are only needed to define the right temporary type through 00092 // typename traits<permutation_sparse_matrix_product<Rhs,Lhs,OnTheRight,false> >::ReturnType 00093 // whereas it should be correctly handled by traits<Product<> >::PlainObject 00094 00095 template<typename Lhs, typename Rhs, int ProductTag> 00096 struct product_evaluator<Product<Lhs, Rhs, AliasFreeProduct>, ProductTag, PermutationShape, SparseShape> 00097 : public evaluator<typename permutation_matrix_product<Rhs,OnTheLeft,false,SparseShape>::ReturnType> 00098 { 00099 typedef Product<Lhs, Rhs, AliasFreeProduct> XprType; 00100 typedef typename permutation_matrix_product<Rhs,OnTheLeft,false,SparseShape>::ReturnType PlainObject; 00101 typedef evaluator<PlainObject> Base; 00102 00103 enum { 00104 Flags = Base::Flags | EvalBeforeNestingBit 00105 }; 00106 00107 explicit product_evaluator(const XprType& xpr) 00108 : m_result(xpr.rows(), xpr.cols()) 00109 { 00110 ::new (static_cast<Base*>(this)) Base(m_result); 00111 generic_product_impl<Lhs, Rhs, PermutationShape, SparseShape, ProductTag>::evalTo(m_result, xpr.lhs(), xpr.rhs()); 00112 } 00113 00114 protected: 00115 PlainObject m_result; 00116 }; 00117 00118 template<typename Lhs, typename Rhs, int ProductTag> 00119 struct product_evaluator<Product<Lhs, Rhs, AliasFreeProduct>, ProductTag, SparseShape, PermutationShape > 00120 : public evaluator<typename permutation_matrix_product<Lhs,OnTheRight,false,SparseShape>::ReturnType> 00121 { 00122 typedef Product<Lhs, Rhs, AliasFreeProduct> XprType; 00123 typedef typename permutation_matrix_product<Lhs,OnTheRight,false,SparseShape>::ReturnType PlainObject; 00124 typedef evaluator<PlainObject> Base; 00125 00126 enum { 00127 Flags = Base::Flags | EvalBeforeNestingBit 00128 }; 00129 00130 explicit product_evaluator(const XprType& xpr) 00131 : m_result(xpr.rows(), xpr.cols()) 00132 { 00133 ::new (static_cast<Base*>(this)) Base(m_result); 00134 generic_product_impl<Lhs, Rhs, SparseShape, PermutationShape, ProductTag>::evalTo(m_result, xpr.lhs(), xpr.rhs()); 00135 } 00136 00137 protected: 00138 PlainObject m_result; 00139 }; 00140 00141 } // end namespace internal 00142 00145 template<typename SparseDerived, typename PermDerived> 00146 inline const Product<SparseDerived, PermDerived, AliasFreeProduct> 00147 operator*(const SparseMatrixBase<SparseDerived>& matrix, const PermutationBase<PermDerived>& perm) 00148 { return Product<SparseDerived, PermDerived, AliasFreeProduct>(matrix.derived(), perm.derived()); } 00149 00152 template<typename SparseDerived, typename PermDerived> 00153 inline const Product<PermDerived, SparseDerived, AliasFreeProduct> 00154 operator*( const PermutationBase<PermDerived>& perm, const SparseMatrixBase<SparseDerived>& matrix) 00155 { return Product<PermDerived, SparseDerived, AliasFreeProduct>(perm.derived(), matrix.derived()); } 00156 00157 00160 template<typename SparseDerived, typename PermutationType> 00161 inline const Product<SparseDerived, Inverse<PermutationType>, AliasFreeProduct> 00162 operator*(const SparseMatrixBase<SparseDerived>& matrix, const InverseImpl<PermutationType, PermutationStorage>& tperm) 00163 { 00164 return Product<SparseDerived, Inverse<PermutationType>, AliasFreeProduct>(matrix.derived(), tperm.derived()); 00165 } 00166 00169 template<typename SparseDerived, typename PermutationType> 00170 inline const Product<Inverse<PermutationType>, SparseDerived, AliasFreeProduct> 00171 operator*(const InverseImpl<PermutationType,PermutationStorage>& tperm, const SparseMatrixBase<SparseDerived>& matrix) 00172 { 00173 return Product<Inverse<PermutationType>, SparseDerived, AliasFreeProduct>(tperm.derived(), matrix.derived()); 00174 } 00175 00176 } // end namespace Eigen 00177 00178 #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H