improve sparse permutations
diff --git a/Eigen/src/SparseCore/SparsePermutation.h b/Eigen/src/SparseCore/SparsePermutation.h index af9a1fe..7e402cc 100644 --- a/Eigen/src/SparseCore/SparsePermutation.h +++ b/Eigen/src/SparseCore/SparsePermutation.h
@@ -18,69 +18,142 @@ namespace internal { +template<typename ExpressionType, typename PlainObjectType, bool NeedEval = !is_same<ExpressionType, PlainObjectType>::value> +struct XprHelper +{ + XprHelper(const ExpressionType& xpr) : m_xpr(xpr) {} + inline const PlainObjectType& xpr() const { return m_xpr; } + // this is a new PlainObjectType initialized by xpr + const PlainObjectType m_xpr; +}; +template<typename ExpressionType, typename PlainObjectType> +struct XprHelper<ExpressionType, PlainObjectType, false> +{ + XprHelper(const ExpressionType& xpr) : m_xpr(xpr) {} + inline const PlainObjectType& xpr() const { return m_xpr; } + // this is a reference to xpr + const PlainObjectType& m_xpr; +}; + +template<typename PermDerived, bool NeedInverseEval> +struct PermHelper +{ + using IndicesType = typename PermDerived::IndicesType; + using PermutationIndex = typename IndicesType::Scalar; + using type = PermutationMatrix<IndicesType::SizeAtCompileTime, IndicesType::MaxSizeAtCompileTime, PermutationIndex>; + PermHelper(const PermDerived& perm) : m_perm(perm.inverse()) {} + inline const type& perm() const { return m_perm; } + // this is a new PermutationMatrix initialized by perm.inverse() + const type m_perm; +}; +template<typename PermDerived> +struct PermHelper<PermDerived, false> +{ + using type = PermDerived; + PermHelper(const PermDerived& perm) : m_perm(perm) {} + inline const type& perm() const { return m_perm; } + // this is a reference to perm + const type& m_perm; +}; + template<typename ExpressionType, int Side, bool Transposed> struct permutation_matrix_product<ExpressionType, Side, Transposed, SparseShape> { - typedef typename nested_eval<ExpressionType, 1>::type MatrixType; - typedef remove_all_t<MatrixType> MatrixTypeCleaned; + using MatrixType = typename nested_eval<ExpressionType, 1>::type; + using MatrixTypeCleaned = remove_all_t<MatrixType>; - typedef typename MatrixTypeCleaned::Scalar Scalar; - typedef typename MatrixTypeCleaned::StorageIndex StorageIndex; + using Scalar = typename MatrixTypeCleaned::Scalar; + using StorageIndex = typename MatrixTypeCleaned::StorageIndex; - enum { - SrcStorageOrder = MatrixTypeCleaned::Flags&RowMajorBit ? RowMajor : ColMajor, - MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight - }; - - typedef std::conditional_t<MoveOuter, - SparseMatrix<Scalar,SrcStorageOrder,StorageIndex>, - SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,StorageIndex> > ReturnType; + // the actual "return type" is `Dest`. this is a temporary type + using ReturnType = SparseMatrix<Scalar, MatrixTypeCleaned::IsRowMajor ? RowMajor : ColMajor, StorageIndex>; + using TmpHelper = XprHelper<ExpressionType, ReturnType>; - template<typename Dest,typename PermutationType> - static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) - { - MatrixType mat(xpr); - if(MoveOuter) - { - SparseMatrix<Scalar,SrcStorageOrder,StorageIndex> tmp(mat.rows(), mat.cols()); - Matrix<StorageIndex,Dynamic,1> sizes(mat.outerSize()); - for(Index j=0; j<mat.outerSize(); ++j) - { + static constexpr bool NeedOuterPermutation = ExpressionType::IsRowMajor ? Side == OnTheLeft : Side == OnTheRight; + static constexpr bool NeedInversePermutation = Transposed ? Side == OnTheLeft : Side == OnTheRight; + + template <typename Dest, typename PermutationType> + static inline void permute_outer(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) { + + // if ExpressionType is not ReturnType, evaluate `xpr` (allocation) + // otherwise, just reference `xpr` + // TODO: handle trivial expressions such as CwiseBinaryOp without temporary + const TmpHelper tmpHelper(xpr); + const ReturnType& tmp = tmpHelper.xpr(); + + ReturnType result(tmp.rows(), tmp.cols()); + + for (Index j = 0; j < tmp.outerSize(); j++) { Index jp = perm.indices().coeff(j); - sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = StorageIndex(mat.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).nonZeros()); + Index jsrc = NeedInversePermutation ? jp : j; + Index jdst = NeedInversePermutation ? j : jp; + Index begin = tmp.outerIndexPtr()[jsrc]; + Index end = tmp.isCompressed() ? tmp.outerIndexPtr()[jsrc + 1] : begin + tmp.innerNonZeroPtr()[jsrc]; + result.outerIndexPtr()[jdst + 1] += end - begin; } - tmp.reserve(sizes); - for(Index j=0; j<mat.outerSize(); ++j) - { - Index jp = perm.indices().coeff(j); - Index jsrc = ((Side==OnTheRight) ^ Transposed) ? jp : j; - Index jdst = ((Side==OnTheLeft) ^ Transposed) ? jp : j; - for(typename MatrixTypeCleaned::InnerIterator it(mat,jsrc); it; ++it) - tmp.insertByOuterInner(jdst,it.index()) = it.value(); - } - dst = tmp; - } - else - { - SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,StorageIndex> tmp(mat.rows(), mat.cols()); - Matrix<StorageIndex,Dynamic,1> sizes(tmp.outerSize()); - sizes.setZero(); - PermutationMatrix<Dynamic,Dynamic,StorageIndex> perm_cpy; - if((Side==OnTheLeft) ^ Transposed) - perm_cpy = perm; - else - perm_cpy = perm.transpose(); - for(Index j=0; j<mat.outerSize(); ++j) - for(typename MatrixTypeCleaned::InnerIterator it(mat,j); it; ++it) - sizes[perm_cpy.indices().coeff(it.index())]++; - tmp.reserve(sizes); - for(Index j=0; j<mat.outerSize(); ++j) - for(typename MatrixTypeCleaned::InnerIterator it(mat,j); it; ++it) - tmp.insertByOuterInner(perm_cpy.indices().coeff(it.index()),j) = it.value(); - dst = tmp; - } + std::partial_sum(result.outerIndexPtr(), result.outerIndexPtr() + result.outerSize() + 1, + result.outerIndexPtr()); + result.resizeNonZeros(result.nonZeros()); + + for (Index j = 0; j < tmp.outerSize(); j++) { + Index jp = perm.indices().coeff(j); + Index jsrc = NeedInversePermutation ? jp : j; + Index jdst = NeedInversePermutation ? j : jp; + Index begin = tmp.outerIndexPtr()[jsrc]; + Index end = tmp.isCompressed() ? tmp.outerIndexPtr()[jsrc + 1] : begin + tmp.innerNonZeroPtr()[jsrc]; + Index target = result.outerIndexPtr()[jdst]; + smart_copy(tmp.innerIndexPtr() + begin, tmp.innerIndexPtr() + end, result.innerIndexPtr() + target); + smart_copy(tmp.valuePtr() + begin, tmp.valuePtr() + end, result.valuePtr() + target); + } + dst = std::move(result); } + + template <typename Dest, typename PermutationType> + static inline void permute_inner(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) { + using InnerPermHelper = PermHelper<PermutationType, NeedInversePermutation>; + using InnerPermType = typename InnerPermHelper::type; + + // if ExpressionType is not ReturnType, evaluate `xpr` (allocation) + // otherwise, just reference `xpr` + // TODO: handle trivial expressions such as CwiseBinaryOp without temporary + const TmpHelper tmpHelper(xpr); + const ReturnType& tmp = tmpHelper.xpr(); + + // if inverse permutation of inner indices is requested, calculate perm.inverse() (allocation) + // otherwise, just reference `perm` + const InnerPermHelper permHelper(perm); + const InnerPermType& innerPerm = permHelper.perm(); + + ReturnType result(tmp.rows(), tmp.cols()); + + for (Index j = 0; j < tmp.outerSize(); j++) { + Index begin = tmp.outerIndexPtr()[j]; + Index end = tmp.isCompressed() ? tmp.outerIndexPtr()[j + 1] : begin + tmp.innerNonZeroPtr()[j]; + result.outerIndexPtr()[j + 1] += end - begin; + } + + std::partial_sum(result.outerIndexPtr(), result.outerIndexPtr() + result.outerSize() + 1, result.outerIndexPtr()); + result.resizeNonZeros(result.nonZeros()); + + for (Index j = 0; j < tmp.outerSize(); j++) { + Index begin = tmp.outerIndexPtr()[j]; + Index end = tmp.isCompressed() ? tmp.outerIndexPtr()[j + 1] : begin + tmp.innerNonZeroPtr()[j]; + Index target = result.outerIndexPtr()[j]; + std::transform(tmp.innerIndexPtr() + begin, tmp.innerIndexPtr() + end, result.innerIndexPtr() + target, + [&innerPerm](StorageIndex i) { return innerPerm.indices().coeff(i); }); + smart_copy(tmp.valuePtr() + begin, tmp.valuePtr() + end, result.valuePtr() + target); + } + // the inner indices were permuted, and must be sorted + result.sortInnerIndices(); + dst = std::move(result); + } + + template <typename Dest, typename PermutationType, bool DoOuter = NeedOuterPermutation, std::enable_if_t<DoOuter, int> = 0> + static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) { permute_outer(dst, perm, xpr); } + + template <typename Dest, typename PermutationType, bool DoOuter = NeedOuterPermutation, std::enable_if_t<!DoOuter, int> = 0> + static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) { permute_inner(dst, perm, xpr); } }; } @@ -143,35 +216,34 @@ } // end namespace internal /** \returns the matrix with the permutation applied to the columns - */ -template<typename SparseDerived, typename PermDerived> -inline const Product<SparseDerived, PermDerived, AliasFreeProduct> -operator*(const SparseMatrixBase<SparseDerived>& matrix, const PermutationBase<PermDerived>& perm) -{ return Product<SparseDerived, PermDerived, AliasFreeProduct>(matrix.derived(), perm.derived()); } + */ +template <typename SparseDerived, typename PermDerived> +inline const Product<SparseDerived, PermDerived, AliasFreeProduct> operator*( + const SparseMatrixBase<SparseDerived>& matrix, const PermutationBase<PermDerived>& perm) { + return Product<SparseDerived, PermDerived, AliasFreeProduct>(matrix.derived(), perm.derived()); +} /** \returns the matrix with the permutation applied to the rows - */ -template<typename SparseDerived, typename PermDerived> -inline const Product<PermDerived, SparseDerived, AliasFreeProduct> -operator*( const PermutationBase<PermDerived>& perm, const SparseMatrixBase<SparseDerived>& matrix) -{ return Product<PermDerived, SparseDerived, AliasFreeProduct>(perm.derived(), matrix.derived()); } - + */ +template <typename SparseDerived, typename PermDerived> +inline const Product<PermDerived, SparseDerived, AliasFreeProduct> operator*( + const PermutationBase<PermDerived>& perm, const SparseMatrixBase<SparseDerived>& matrix) { + return Product<PermDerived, SparseDerived, AliasFreeProduct>(perm.derived(), matrix.derived()); +} /** \returns the matrix with the inverse permutation applied to the columns. - */ -template<typename SparseDerived, typename PermutationType> -inline const Product<SparseDerived, Inverse<PermutationType>, AliasFreeProduct> -operator*(const SparseMatrixBase<SparseDerived>& matrix, const InverseImpl<PermutationType, PermutationStorage>& tperm) -{ + */ +template <typename SparseDerived, typename PermutationType> +inline const Product<SparseDerived, Inverse<PermutationType>, AliasFreeProduct> operator*( + const SparseMatrixBase<SparseDerived>& matrix, const InverseImpl<PermutationType, PermutationStorage>& tperm) { return Product<SparseDerived, Inverse<PermutationType>, AliasFreeProduct>(matrix.derived(), tperm.derived()); } /** \returns the matrix with the inverse permutation applied to the rows. - */ -template<typename SparseDerived, typename PermutationType> -inline const Product<Inverse<PermutationType>, SparseDerived, AliasFreeProduct> -operator*(const InverseImpl<PermutationType,PermutationStorage>& tperm, const SparseMatrixBase<SparseDerived>& matrix) -{ + */ +template <typename SparseDerived, typename PermutationType> +inline const Product<Inverse<PermutationType>, SparseDerived, AliasFreeProduct> operator*( + const InverseImpl<PermutationType, PermutationStorage>& tperm, const SparseMatrixBase<SparseDerived>& matrix) { return Product<Inverse<PermutationType>, SparseDerived, AliasFreeProduct>(tperm.derived(), matrix.derived()); }
diff --git a/test/sparse_permutations.cpp b/test/sparse_permutations.cpp index 5974c74..775c560 100644 --- a/test/sparse_permutations.cpp +++ b/test/sparse_permutations.cpp
@@ -17,6 +17,15 @@ VERIFY( (#XPR) && nb_transposed_copies==N ); \ } +static long int nb_temporaries; +#define EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN {nb_temporaries++;} +#define VERIFY_TEMPORARY_COUNT(XPR,N) {\ + nb_temporaries = 0; \ + XPR; \ + if(nb_temporaries!=N) std::cerr << "nb_temporaries == " << nb_temporaries << "\n"; \ + VERIFY( (#XPR) && nb_temporaries==N ); \ + } + #include "sparse.h" template<typename T> @@ -79,28 +88,52 @@ VERIFY( is_sorted( ::eval(mat*p) )); VERIFY( is_sorted( res = mat*p )); VERIFY_TRANSPOSITION_COUNT( ::eval(mat*p), 0); - //VERIFY_TRANSPOSITION_COUNT( res = mat*p, IsRowMajor ? 1 : 0 ); + VERIFY_TEMPORARY_COUNT( ::eval(mat*p), 1) res_d = mat_d*p; VERIFY(res.isApprox(res_d) && "mat*p"); VERIFY( is_sorted( ::eval(p*mat) )); VERIFY( is_sorted( res = p*mat )); VERIFY_TRANSPOSITION_COUNT( ::eval(p*mat), 0); + VERIFY_TEMPORARY_COUNT( ::eval(p*mat), 1); res_d = p*mat_d; VERIFY(res.isApprox(res_d) && "p*mat"); VERIFY( is_sorted( (mat*p).eval() )); VERIFY( is_sorted( res = mat*p.inverse() )); VERIFY_TRANSPOSITION_COUNT( ::eval(mat*p.inverse()), 0); + VERIFY_TEMPORARY_COUNT( ::eval(mat*p.inverse()), 1); res_d = mat*p.inverse(); VERIFY(res.isApprox(res_d) && "mat*inv(p)"); VERIFY( is_sorted( (p*mat+p*mat).eval() )); VERIFY( is_sorted( res = p.inverse()*mat )); VERIFY_TRANSPOSITION_COUNT( ::eval(p.inverse()*mat), 0); + VERIFY_TEMPORARY_COUNT( ::eval(p.inverse()*mat), 1); res_d = p.inverse()*mat_d; VERIFY(res.isApprox(res_d) && "inv(p)*mat"); + // test non-plaintype expressions that require additional temporary + const Scalar alpha(2.34); + + res_d = p * (alpha * mat_d); + VERIFY_TEMPORARY_COUNT( res = p * (alpha * mat), 2); + VERIFY( res.isApprox(res_d) && "p*(alpha*mat)" ); + + res_d = (alpha * mat_d) * p; + VERIFY_TEMPORARY_COUNT( res = (alpha * mat) * p, 2); + VERIFY( res.isApprox(res_d) && "(alpha*mat)*p" ); + + res_d = p.inverse() * (alpha * mat_d); + VERIFY_TEMPORARY_COUNT( res = p.inverse() * (alpha * mat), 2); + VERIFY( res.isApprox(res_d) && "inv(p)*(alpha*mat)" ); + + res_d = (alpha * mat_d) * p.inverse(); + VERIFY_TEMPORARY_COUNT( res = (alpha * mat) * p.inverse(), 2); + VERIFY( res.isApprox(res_d) && "(alpha*mat)*inv(p)" ); + + // + VERIFY( is_sorted( (p * mat * p.inverse()).eval() )); VERIFY( is_sorted( res = mat.twistedBy(p) )); VERIFY_TRANSPOSITION_COUNT( ::eval(p * mat * p.inverse()), 0);
diff --git a/unsupported/Eigen/src/SparseExtra/MarketIO.h b/unsupported/Eigen/src/SparseExtra/MarketIO.h index cf5828e..5c26133 100644 --- a/unsupported/Eigen/src/SparseExtra/MarketIO.h +++ b/unsupported/Eigen/src/SparseExtra/MarketIO.h
@@ -189,6 +189,7 @@ readsizes = true; mat.resize(M,N); mat.reserve(NNZ); + elements.reserve(NNZ); } } else