From fa0bd2c34ec7cb8590fa157a209d2002c8116f03 Mon Sep 17 00:00:00 2001 From: Charles Schlosser Date: Sun, 15 Jan 2023 03:21:25 +0000 Subject: [PATCH] improve sparse permutations --- Eigen/src/SparseCore/SparsePermutation.h | 220 ++++++++++++------- test/sparse_permutations.cpp | 35 ++- unsupported/Eigen/src/SparseExtra/MarketIO.h | 1 + 3 files changed, 181 insertions(+), 75 deletions(-) diff --git a/Eigen/src/SparseCore/SparsePermutation.h b/Eigen/src/SparseCore/SparsePermutation.h index af9a1fea6..7e402cced 100644 --- a/Eigen/src/SparseCore/SparsePermutation.h +++ b/Eigen/src/SparseCore/SparsePermutation.h @@ -18,69 +18,142 @@ namespace Eigen { namespace internal { +template::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 +struct XprHelper +{ + 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 +struct PermHelper +{ + using IndicesType = typename PermDerived::IndicesType; + using PermutationIndex = typename IndicesType::Scalar; + using type = PermutationMatrix; + 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 +struct PermHelper +{ + 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 struct permutation_matrix_product { - typedef typename nested_eval::type MatrixType; - typedef remove_all_t MatrixTypeCleaned; + using MatrixType = typename nested_eval::type; + using MatrixTypeCleaned = remove_all_t; - 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, - SparseMatrix > ReturnType; + // the actual "return type" is `Dest`. this is a temporary type + using ReturnType = SparseMatrix; + using TmpHelper = XprHelper; - template - static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) - { - MatrixType mat(xpr); - if(MoveOuter) - { - SparseMatrix tmp(mat.rows(), mat.cols()); - Matrix sizes(mat.outerSize()); - for(Index j=0; j + 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 tmp(mat.rows(), mat.cols()); - Matrix sizes(tmp.outerSize()); - sizes.setZero(); - PermutationMatrix perm_cpy; - if((Side==OnTheLeft) ^ Transposed) - perm_cpy = perm; - else - perm_cpy = perm.transpose(); - for(Index j=0; j + static inline void permute_inner(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) { + using InnerPermHelper = PermHelper; + 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 = 0> + static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) { permute_outer(dst, perm, xpr); } + + template = 0> + static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) { permute_inner(dst, perm, xpr); } }; } @@ -143,35 +216,34 @@ protected: } // end namespace internal /** \returns the matrix with the permutation applied to the columns - */ -template -inline const Product -operator*(const SparseMatrixBase& matrix, const PermutationBase& perm) -{ return Product(matrix.derived(), perm.derived()); } + */ +template +inline const Product operator*( + const SparseMatrixBase& matrix, const PermutationBase& perm) { + return Product(matrix.derived(), perm.derived()); +} /** \returns the matrix with the permutation applied to the rows - */ -template -inline const Product -operator*( const PermutationBase& perm, const SparseMatrixBase& matrix) -{ return Product(perm.derived(), matrix.derived()); } - + */ +template +inline const Product operator*( + const PermutationBase& perm, const SparseMatrixBase& matrix) { + return Product(perm.derived(), matrix.derived()); +} /** \returns the matrix with the inverse permutation applied to the columns. - */ -template -inline const Product, AliasFreeProduct> -operator*(const SparseMatrixBase& matrix, const InverseImpl& tperm) -{ + */ +template +inline const Product, AliasFreeProduct> operator*( + const SparseMatrixBase& matrix, const InverseImpl& tperm) { return Product, AliasFreeProduct>(matrix.derived(), tperm.derived()); } /** \returns the matrix with the inverse permutation applied to the rows. - */ -template -inline const Product, SparseDerived, AliasFreeProduct> -operator*(const InverseImpl& tperm, const SparseMatrixBase& matrix) -{ + */ +template +inline const Product, SparseDerived, AliasFreeProduct> operator*( + const InverseImpl& tperm, const SparseMatrixBase& matrix) { return Product, SparseDerived, AliasFreeProduct>(tperm.derived(), matrix.derived()); } diff --git a/test/sparse_permutations.cpp b/test/sparse_permutations.cpp index 5974c74f8..775c5608b 100644 --- a/test/sparse_permutations.cpp +++ b/test/sparse_permutations.cpp @@ -17,6 +17,15 @@ static long int nb_transposed_copies; 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 @@ -79,28 +88,52 @@ template void sparse_permutations(c 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 cf5828e18..5c2613333 100644 --- a/unsupported/Eigen/src/SparseExtra/MarketIO.h +++ b/unsupported/Eigen/src/SparseExtra/MarketIO.h @@ -189,6 +189,7 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename) readsizes = true; mat.resize(M,N); mat.reserve(NNZ); + elements.reserve(NNZ); } } else