From d517c89a4d6ce139444ad752b67d4be0d326effe Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 26 Feb 2024 21:09:21 +0100 Subject: [PATCH 1/6] [LinearAlgebra] Refactor sparse matrix product --- Sofa/framework/LinearAlgebra/CMakeLists.txt | 6 +- .../SparseMatrixTest.h | 57 ++-- ...parseMatrix].h => SparseMatrixProduct.cpp} | 13 +- .../sofa/linearalgebra/SparseMatrixProduct.h | 71 ++--- .../linearalgebra/SparseMatrixProduct.inl | 279 ++++++++++++++++++ ...trixProduct[CompressedRowSparseMatrix].cpp | 43 --- ...SparseMatrixProduct[EigenSparseMatrix].cpp | 262 ---------------- .../SparseMatrixProduct[EigenSparseMatrix].h | 43 --- .../test/SparseMatrixProduct_test.cpp | 42 +-- 9 files changed, 351 insertions(+), 465 deletions(-) rename Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/{SparseMatrixProduct[CompressedRowSparseMatrix].h => SparseMatrixProduct.cpp} (78%) create mode 100644 Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.inl delete mode 100644 Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[CompressedRowSparseMatrix].cpp delete mode 100644 Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[EigenSparseMatrix].cpp delete mode 100644 Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[EigenSparseMatrix].h diff --git a/Sofa/framework/LinearAlgebra/CMakeLists.txt b/Sofa/framework/LinearAlgebra/CMakeLists.txt index 34f87cc30cd..f97ec21cbed 100644 --- a/Sofa/framework/LinearAlgebra/CMakeLists.txt +++ b/Sofa/framework/LinearAlgebra/CMakeLists.txt @@ -37,8 +37,7 @@ set(HEADER_FILES ${SOFALINEARALGEBRASRC_ROOT}/RotationMatrix.h ${SOFALINEARALGEBRASRC_ROOT}/SparseMatrix.h ${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixProduct.h - ${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixProduct[CompressedRowSparseMatrix].h - ${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixProduct[EigenSparseMatrix].h + ${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixProduct.inl ${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixStorageOrder.h ${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixStorageOrder[EigenSparseMatrix].h ${SOFALINEARALGEBRASRC_ROOT}/TriangularSystemSolver.h @@ -59,8 +58,7 @@ set(SOURCE_FILES ${SOFALINEARALGEBRASRC_ROOT}/FullMatrix.cpp ${SOFALINEARALGEBRASRC_ROOT}/FullVector.cpp ${SOFALINEARALGEBRASRC_ROOT}/RotationMatrix.cpp - ${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixProduct[CompressedRowSparseMatrix].cpp - ${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixProduct[EigenSparseMatrix].cpp + ${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixProduct.cpp ${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixStorageOrder[EigenSparseMatrix].cpp ${SOFALINEARALGEBRASRC_ROOT}/init.cpp ) diff --git a/Sofa/framework/LinearAlgebra/Testing/src/Sofa.LinearAlgebra.Testing/SparseMatrixTest.h b/Sofa/framework/LinearAlgebra/Testing/src/Sofa.LinearAlgebra.Testing/SparseMatrixTest.h index 7d92c5bf313..e654875d821 100644 --- a/Sofa/framework/LinearAlgebra/Testing/src/Sofa.LinearAlgebra.Testing/SparseMatrixTest.h +++ b/Sofa/framework/LinearAlgebra/Testing/src/Sofa.LinearAlgebra.Testing/SparseMatrixTest.h @@ -68,12 +68,11 @@ struct SparseMatrixTest : public virtual NumericTest eigenMatrix.setFromTriplets(first, last); } - static void copyFromEigen(Eigen::SparseMatrix& dst, const Eigen::SparseMatrix& src) - { - dst = src; - } - - static void copyFromEigen(Eigen::SparseMatrix& dst, const Eigen::SparseMatrix& src) + template< + typename _DstScalar, int _DstOptions, typename _DstStorageIndex, + typename _SrcScalar, int _SrcOptions, typename _SrcStorageIndex + > + static void copyFromEigen(Eigen::SparseMatrix<_DstScalar, _DstOptions, _DstStorageIndex>& dst, const Eigen::SparseMatrix<_SrcScalar, _SrcOptions, _SrcStorageIndex>& src) { dst = src; } @@ -93,48 +92,26 @@ struct SparseMatrixTest : public virtual NumericTest } } - static bool compareSparseMatrix(const Eigen::SparseMatrix& A, const Eigen::SparseMatrix& B) + template< + typename _AScalar, int _AOptions, typename _AStorageIndex, + typename _BScalar, int _BOptions, typename _BStorageIndex + > + static bool compareSparseMatrix(const Eigen::SparseMatrix<_AScalar, _AOptions, _AStorageIndex>& A, const Eigen::SparseMatrix<_BScalar, _BOptions, _BStorageIndex>& B) { return compareEigenSparseMatrix(A, B); } - static bool compareEigenSparseMatrix(const Eigen::SparseMatrix& A, const Eigen::SparseMatrix& B) + template< + typename _AScalar, int _AOptions, typename _AStorageIndex, + typename _BScalar, int _BOptions, typename _BStorageIndex + > + static bool compareEigenSparseMatrix(const Eigen::SparseMatrix<_AScalar, _AOptions, _AStorageIndex>& A, const Eigen::SparseMatrix<_BScalar, _BOptions, _BStorageIndex>& B) { - if (A.outerSize() != B.outerSize()) - return false; - - for (int k = 0; k < A.outerSize(); ++k) - { - sofa::type::vector > triplets_a, triplets_b; - for (typename Eigen::SparseMatrix::InnerIterator it(A, k); it; ++it) - { - triplets_a.emplace_back(it.row(), it.col(), it.value()); - } - for (typename Eigen::SparseMatrix::InnerIterator it(B, k); it; ++it) - { - triplets_b.emplace_back(it.row(), it.col(), it.value()); - } - - if (triplets_a.size() != triplets_b.size()) - return false; - - for (size_t i = 0 ; i < triplets_a.size(); ++i) - { - const auto& a = triplets_a[i]; - const auto& b = triplets_b[i]; - - if (a.row() != b.row() || a.col() != b.col() || !NumericTest::isSmall(a.value() - b.value(), 1)) - { - // ret = false; - return false; - } - } - } - return true; + return A.isApprox(B); } }; -} // namespace sofa::testing \ No newline at end of file +} // namespace sofa::testing diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[CompressedRowSparseMatrix].h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.cpp similarity index 78% rename from Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[CompressedRowSparseMatrix].h rename to Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.cpp index 57f1965da4e..7f50530ceef 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[CompressedRowSparseMatrix].h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.cpp @@ -19,17 +19,8 @@ * * * Contact information: contact@sofa-framework.org * ******************************************************************************/ -#pragma once - -#include -#include +#define SOFA_LINEARALGEBRA_SPARSEMATRIXPRODUCT_DEFINITION namespace sofa::linearalgebra { - -#if !defined(SOFA_LINEARAGEBRA_SPARSEMATRIXPRODUCT_COMPRESSEDROWSPARSEMATRIX_CPP) - extern template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; - extern template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; -#endif - -} //namespace sofa::linearalgebra \ No newline at end of file +} diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.h index b74528466ee..e9b86afb8e8 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.h @@ -23,8 +23,8 @@ #include #include - -#include +#include +#include namespace sofa::linearalgebra { @@ -44,40 +44,54 @@ namespace sofa::linearalgebra * and * Saupin, G., 2008. Vers la simulation interactive réaliste de corps déformables virtuels (Doctoral dissertation, Lille 1). */ -template +template class SparseMatrixProduct { public: + + using LhsCleaned = std::decay_t; + using RhsCleaned = std::decay_t; + using ResultCleaned = std::decay_t; + + using LhsScalar = typename LhsCleaned::Scalar; + using RhsScalar = typename RhsCleaned::Scalar; + using ResultScalar = typename ResultCleaned::Scalar; + /// Left side of the product A*B - const TMatrix* matrixA { nullptr }; + const LhsCleaned* m_lhs { nullptr }; /// Right side of the product A*B - const TMatrix* matrixB { nullptr }; + const RhsCleaned* m_rhs { nullptr }; + + using Index = Eigen::Index; + + using ProductResult = ResultType; + void computeProduct(bool forceComputeIntersection = false); void computeRegularProduct(); - const TMatrix& getProductResult() const { return matrixC; } + [[nodiscard]] const ResultType& getProductResult() const { return m_productResult; } void invalidateIntersection(); - SparseMatrixProduct(TMatrix* a, TMatrix* b) : matrixA(a), matrixB(b) {} + SparseMatrixProduct(Lhs* a, Rhs* b) : m_lhs(a), m_rhs(b) {} SparseMatrixProduct() = default; struct Intersection { // Two indices: the first for the values vector of the matrix A, the second for the values vector of the matrix B - using PairIndex = std::pair; + using PairIndex = std::pair; // A list of pairs of indices - using ListPairIndex = type::vector; + using ListPairIndex = sofa::type::vector; /// Each element of this vector gives the list of values from matrix A and B to multiply together and accumulate /// them into the matrix C at the same location in the values vector - using ValuesIntersection = type::vector; + using ValuesIntersection = sofa::type::vector; ValuesIntersection intersection; }; protected: - TMatrix matrixC; /// Result of A*B + ProductResult m_productResult; /// Result of A*B bool m_hasComputedIntersection { false }; void computeIntersection(); @@ -87,38 +101,5 @@ class SparseMatrixProduct }; -template -void SparseMatrixProduct::computeProduct(bool forceComputeIntersection) -{ - if (forceComputeIntersection) - { - m_hasComputedIntersection = false; - } - - if (m_hasComputedIntersection == false) - { - computeIntersection(); - m_hasComputedIntersection = true; - } - computeProductFromIntersection(); -} - -template -void SparseMatrixProduct::computeIntersection() -{ - -} - -template -void SparseMatrixProduct::computeProductFromIntersection() -{ - matrixC = (*matrixA) * (*matrixB); -} - -template -void SparseMatrixProduct::invalidateIntersection() -{ - m_hasComputedIntersection = false; -} -}// sofa::linearalgebra \ No newline at end of file +}// sofa::linearalgebra diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.inl b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.inl new file mode 100644 index 00000000000..070619075d0 --- /dev/null +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.inl @@ -0,0 +1,279 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include +#include +#include +#include + + +namespace sofa::linearalgebra::sparsematrixproduct +{ + +/** + * Represent a scalar and its index in an array of scalars + */ +template +struct IndexedValue +{ + Eigen::Index index {}; + Scalar value; + + IndexedValue() = default; + + template > > + IndexedValue(AnyScalar s) : value(s) {} + + IndexedValue(const IndexedValue& other) = default; + + operator Scalar() const + { + return value; + } +}; + +template +std::ostream& operator<<(std::ostream& o, IndexedValue& p) +{ + o << "(" << p.value << ", " << p.index << ")"; + return o; +} + + +/** + * Represent a sum of scalar products. It stores: + * - a value for the result + * - a list of pairs of indices to know what scalars were used for the computation + */ +template +class IndexValueProduct +{ +private: + using IndexLHS = Eigen::Index; + using IndexRHS = Eigen::Index; + + using ScalarProduct = std::pair; + sofa::type::vector m_indices {}; + + Scalar value {}; + +public: + + [[nodiscard]] const sofa::type::vector& getIndices() const + { + return m_indices; + } + IndexValueProduct() = default; + + template > > + IndexValueProduct(AnyScalar s) : value(s) {} + + operator Scalar() const + { + return value; + } + + template + IndexValueProduct(const IndexValueProduct& other) + : m_indices(other.indices) + , value(static_cast(other.value)) + {} + + template + void operator+=(const IndexValueProduct& other) + { + m_indices.insert(m_indices.end(), other.m_indices.begin(), other.m_indices.end()); + value += static_cast(other.value); + } + + template + friend IndexValueProduct + operator*(const IndexedValue& lhs, const IndexedValue& rhs); +}; + +template +std::ostream& operator<<(std::ostream& o, IndexValueProduct& p) +{ + o << "(" << p.value << ", [" << p.indices << "])"; + return o; +} + +template +IndexValueProduct +operator*(const IndexedValue& lhs, const IndexedValue& rhs) +{ + IndexValueProduct product; + product.m_indices.resize(1, {lhs.index, rhs.index}); + product.value = lhs.value * rhs.value; + return product; +} + +} + +//this is to inform Eigen that the product of two IndexedValue is a IndexValueProduct +#define DEFINE_PRODUCT_OP_FOR_TYPES(lhs, rhs) \ +template<> \ +struct Eigen::ScalarBinaryOpTraits< \ + sofa::linearalgebra::sparsematrixproduct::IndexedValue, \ + sofa::linearalgebra::sparsematrixproduct::IndexedValue, \ + Eigen::internal::scalar_product_op< \ + sofa::linearalgebra::sparsematrixproduct::IndexedValue, \ + sofa::linearalgebra::sparsematrixproduct::IndexedValue \ + > \ +> \ +{ \ + typedef sofa::linearalgebra::sparsematrixproduct::IndexValueProduct ReturnType; \ +}; + +DEFINE_PRODUCT_OP_FOR_TYPES(float, float) +DEFINE_PRODUCT_OP_FOR_TYPES(double, float) +DEFINE_PRODUCT_OP_FOR_TYPES(float, double) +DEFINE_PRODUCT_OP_FOR_TYPES(double, double) + +namespace sofa::linearalgebra +{ +template +void SparseMatrixProduct::computeProduct(bool forceComputeIntersection) +{ + if (forceComputeIntersection) + { + m_hasComputedIntersection = false; + } + + if (m_hasComputedIntersection == false) + { + computeIntersection(); + m_hasComputedIntersection = true; + } + else + { + computeProductFromIntersection(); + } +} + +template +void SparseMatrixProduct::computeRegularProduct() +{ + m_productResult = *m_lhs * *m_rhs; +} + +template +void flagValueIndices(Eigen::SparseMatrix, _Options, _StorageIndex>& matrix) +{ + for (Eigen::Index i = 0; i < matrix.nonZeros(); ++i) + { + matrix.valuePtr()[i].index = i; + } +} + +template +void SparseMatrixProduct::computeIntersection() +{ + using LocalLhs = Eigen::SparseMatrix< + sparsematrixproduct::IndexedValue, + Lhs::Options, + typename Lhs::StorageIndex + >; + + using LocalRhs = Eigen::SparseMatrix< + sparsematrixproduct::IndexedValue, + Rhs::Options, + typename Rhs::StorageIndex + >; + + //copy the input matrices in an intermediate matrix with the same properties + //except that the type of values is IndexedValue + LocalLhs lhs = m_lhs->template cast>(); + LocalRhs rhs = m_rhs->template cast>(); + + flagValueIndices(lhs); + flagValueIndices(rhs); + + using LocalResult = Eigen::SparseMatrix< + sparsematrixproduct::IndexValueProduct, + ResultType::Options, + typename ResultType::StorageIndex + >; + + const LocalResult product = lhs * rhs; + + const auto productNonZeros = product.nonZeros(); + m_intersectionAB.intersection.clear(); + m_intersectionAB.intersection.reserve(productNonZeros); + + for (Eigen::Index i = 0; i < productNonZeros; ++i) + { + m_intersectionAB.intersection.push_back(product.valuePtr()[i].getIndices()); + + //depending on the storage scheme, Eigen can change the order of the lhs and rhs + if constexpr (Lhs::IsRowMajor || Rhs::IsRowMajor) + { + for (auto& [lhsIndex, rhsIndex] : m_intersectionAB.intersection.back()) + { + std::swap(lhsIndex, rhsIndex); + } + } + } + + m_productResult = product.template cast(); +} + +template +void __computeProductFromIntersection(const Lhs* lhs, const Rhs* rhs, ResultType* product, + const typename SparseMatrixProduct::Intersection& intersection) +{ + assert(intersection.intersection.size() == product->nonZeros()); + + auto* lhs_ptr = lhs->valuePtr(); + auto* rhs_ptr = rhs->valuePtr(); + auto* product_ptr = product->valuePtr(); + + for (const auto& pairs : intersection.intersection) + { + auto& value = *product_ptr++; + value = 0; + for (auto [lhsIndex, rhsIndex] : pairs) + { + value += lhs_ptr[lhsIndex] * rhs_ptr[rhsIndex]; + } + } +} + +template +void SparseMatrixProduct::computeProductFromIntersection() +{ + if (m_hasComputedIntersection == false) + { + msg_error("SparseMatrixProduct") << "Intersection computation is required before computing the product"; + return; + } + __computeProductFromIntersection(m_lhs, m_rhs, &m_productResult, m_intersectionAB); +} + +template +void SparseMatrixProduct::invalidateIntersection() +{ + m_hasComputedIntersection = false; +} + +} diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[CompressedRowSparseMatrix].cpp b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[CompressedRowSparseMatrix].cpp deleted file mode 100644 index 15b59f51580..00000000000 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[CompressedRowSparseMatrix].cpp +++ /dev/null @@ -1,43 +0,0 @@ -/****************************************************************************** -* SOFA, Simulation Open-Framework Architecture * -* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * -* * -* This program is free software; you can redistribute it and/or modify it * -* under the terms of the GNU Lesser General Public License as published by * -* the Free Software Foundation; either version 2.1 of the License, or (at * -* your option) any later version. * -* * -* This program is distributed in the hope that it will be useful, but WITHOUT * -* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * -* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * -* for more details. * -* * -* You should have received a copy of the GNU Lesser General Public License * -* along with this program. If not, see . * -******************************************************************************* -* Authors: The SOFA Team and external contributors (see Authors.txt) * -* * -* Contact information: contact@sofa-framework.org * -******************************************************************************/ -#define SOFA_LINEARAGEBRA_SPARSEMATRIXPRODUCT_COMPRESSEDROWSPARSEMATRIX_CPP -#include - -namespace sofa::linearalgebra -{ - -template <> -void SparseMatrixProduct >::computeProductFromIntersection() -{ - matrixA->mul(matrixC, *matrixB); -} - -template <> -void SparseMatrixProduct >::computeProductFromIntersection() -{ - matrixA->mul(matrixC, *matrixB); -} - -template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; -template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; - -} //namespace sofa::linearalgebra \ No newline at end of file diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[EigenSparseMatrix].cpp b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[EigenSparseMatrix].cpp deleted file mode 100644 index 6f42508d962..00000000000 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[EigenSparseMatrix].cpp +++ /dev/null @@ -1,262 +0,0 @@ -/****************************************************************************** -* SOFA, Simulation Open-Framework Architecture * -* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * -* * -* This program is free software; you can redistribute it and/or modify it * -* under the terms of the GNU Lesser General Public License as published by * -* the Free Software Foundation; either version 2.1 of the License, or (at * -* your option) any later version. * -* * -* This program is distributed in the hope that it will be useful, but WITHOUT * -* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * -* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * -* for more details. * -* * -* You should have received a copy of the GNU Lesser General Public License * -* along with this program. If not, see . * -******************************************************************************* -* Authors: The SOFA Team and external contributors (see Authors.txt) * -* * -* Contact information: contact@sofa-framework.org * -******************************************************************************/ -#define SOFA_LINEARAGEBRA_SPARSEMATRIXPRODUCT_EIGENSPARSEMATRIX_CPP -#include -#include -#include - -#include - -namespace sofa::linearalgebra -{ - -template <> -void SparseMatrixProduct >::computeRegularProduct() -{ - matrixC = (*matrixA) * (*matrixB); -} - -template <> -void SparseMatrixProduct >::computeRegularProduct() -{ - matrixC = (*matrixA) * (*matrixB); -} - -template <> -void SparseMatrixProduct >::computeRegularProduct() -{ - matrixC = (*matrixA) * (*matrixB); -} - -template <> -void SparseMatrixProduct >::computeRegularProduct() -{ - matrixC = (*matrixA) * (*matrixB); -} - -template -void __computeIntersectionColumnMajor(const TMatrix* A, const TMatrix* B, TMatrix* C, typename SparseMatrixProduct::Intersection& intersection) -{ - C->resize(A->rows(), B->cols()); - - *C = (*A) * (*B); - - const SparseMatrixStorageOrder transpose(A); - - const auto& outerStarts = transpose.getOuterStarts(); - const auto& innerIndices = transpose.getInnerIndices(); - const auto& perm = transpose.getPermutations(); - - if (outerStarts.empty()) - return; - - intersection.intersection.clear(); - intersection.intersection.reserve(C->nonZeros()); - - for (Eigen::Index c = 0; c < B->outerSize(); ++c) - { - const auto beginB = B->outerIndexPtr()[c]; - const auto endB = B->outerIndexPtr()[c + 1]; - - for (std::size_t r = 0; r < outerStarts.size() - 1; ++r) - { - const auto beginA = outerStarts[r]; - const auto endA = outerStarts[r + 1]; - - auto iA = beginA; - auto iB = beginB; - typename SparseMatrixProduct::Intersection::ListPairIndex listPairs; - - while( iA < endA && iB < endB) - { - const auto inner_A = innerIndices[iA]; - const auto inner_B = B->innerIndexPtr()[iB]; - if (inner_A == inner_B) //intersection - { - listPairs.emplace_back(perm[iA], iB); - } - - if (inner_A < inner_B) - ++iA; - else - ++iB; - } - - if (!listPairs.empty()) - intersection.intersection.push_back(listPairs); - - } - - } -} - -template -void __computeIntersectionRowMajor(const TMatrix* A, const TMatrix* B, TMatrix* C, typename SparseMatrixProduct::Intersection& intersection) -{ - C->resize(A->rows(), B->cols()); - - *C = (*A) * (*B); - - const SparseMatrixStorageOrder transpose(B); - - const auto& outerStarts = transpose.getOuterStarts(); - const auto& innerIndices = transpose.getInnerIndices(); - const auto& perm = transpose.getPermutations(); - - if (outerStarts.empty()) - return; - - intersection.intersection.clear(); - intersection.intersection.reserve(C->nonZeros()); - - for (Eigen::Index r = 0; r < A->outerSize(); ++r) - { - const auto beginA = A->outerIndexPtr()[r]; - const auto endA = A->outerIndexPtr()[r + 1]; - - for (std::size_t c = 0; c < outerStarts.size() - 1; ++c) - { - const auto beginB = outerStarts[c]; - const auto endB = outerStarts[c + 1]; - - auto iA = beginA; - auto iB = beginB; - typename SparseMatrixProduct::Intersection::ListPairIndex listPairs; - - while( iA < endA && iB < endB) - { - const auto inner_A = A->innerIndexPtr()[iA]; - const auto inner_B = innerIndices[iB]; - if (inner_A == inner_B) //intersection - { - listPairs.emplace_back(iA, perm[iB]); - } - - if (inner_A < inner_B) - ++iA; - else - ++iB; - } - - if (!listPairs.empty()) - intersection.intersection.push_back(listPairs); - - } - - } -} - -template <> -void SparseMatrixProduct >::computeIntersection() -{ - __computeIntersectionColumnMajor(matrixA, matrixB, &matrixC, m_intersectionAB); -} - -template <> -void SparseMatrixProduct >::computeIntersection() -{ - __computeIntersectionColumnMajor(matrixA, matrixB, &matrixC, m_intersectionAB); -} - -template <> -void SparseMatrixProduct >::computeIntersection() -{ - __computeIntersectionRowMajor(matrixA, matrixB, &matrixC, m_intersectionAB); -} - -template <> -void SparseMatrixProduct >::computeIntersection() -{ - __computeIntersectionRowMajor(matrixA, matrixB, &matrixC, m_intersectionAB); -} - -template -void __computeProductFromIntersection(const TMatrix* A, const TMatrix* B, TMatrix* C, const typename SparseMatrixProduct::Intersection& intersection) -{ - assert(intersection.intersection.size() == C->nonZeros()); - - auto* a_ptr = A->valuePtr(); - auto* b_ptr = B->valuePtr(); - auto* c_ptr = C->valuePtr(); - - for (const auto& pairs : intersection.intersection) - { - auto& value = *c_ptr++; - value = 0; - for (const auto& p : pairs) - { - value += a_ptr[p.first] * b_ptr[p.second]; - } - } -} - -template <> -void SparseMatrixProduct >::computeProductFromIntersection() -{ - if (m_hasComputedIntersection == false) - { - msg_error("SparseMatrixProduct") << "Intersection computation is required before computing the product"; - return; - } - __computeProductFromIntersection(matrixA, matrixB, &matrixC, m_intersectionAB); -} - -template <> -void SparseMatrixProduct >::computeProductFromIntersection() -{ - if (m_hasComputedIntersection == false) - { - msg_error("SparseMatrixProduct") << "Intersection computation is required before computing the product"; - return; - } - __computeProductFromIntersection(matrixA, matrixB, &matrixC, m_intersectionAB); -} - -template <> -void SparseMatrixProduct >::computeProductFromIntersection() -{ - if (m_hasComputedIntersection == false) - { - msg_error("SparseMatrixProduct") << "Intersection computation is required before computing the product"; - return; - } - __computeProductFromIntersection(matrixA, matrixB, &matrixC, m_intersectionAB); -} - -template <> -void SparseMatrixProduct >::computeProductFromIntersection() -{ - if (m_hasComputedIntersection == false) - { - msg_error("SparseMatrixProduct") << "Intersection computation is required before computing the product"; - return; - } - __computeProductFromIntersection(matrixA, matrixB, &matrixC, m_intersectionAB); -} - -template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; -template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; - -template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; -template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; - -} //namespace sofa::linearalgebra \ No newline at end of file diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[EigenSparseMatrix].h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[EigenSparseMatrix].h deleted file mode 100644 index e9a7d446839..00000000000 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[EigenSparseMatrix].h +++ /dev/null @@ -1,43 +0,0 @@ -/****************************************************************************** -* SOFA, Simulation Open-Framework Architecture * -* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * -* * -* This program is free software; you can redistribute it and/or modify it * -* under the terms of the GNU Lesser General Public License as published by * -* the Free Software Foundation; either version 2.1 of the License, or (at * -* your option) any later version. * -* * -* This program is distributed in the hope that it will be useful, but WITHOUT * -* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * -* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * -* for more details. * -* * -* You should have received a copy of the GNU Lesser General Public License * -* along with this program. If not, see . * -******************************************************************************* -* Authors: The SOFA Team and external contributors (see Authors.txt) * -* * -* Contact information: contact@sofa-framework.org * -******************************************************************************/ -#pragma once - -#include -#include - -namespace sofa::linearalgebra -{ - -template<> void SOFA_LINEARALGEBRA_API SparseMatrixProduct >::computeRegularProduct(); -template<> void SOFA_LINEARALGEBRA_API SparseMatrixProduct >::computeRegularProduct(); -template<> void SOFA_LINEARALGEBRA_API SparseMatrixProduct >::computeRegularProduct(); -template<> void SOFA_LINEARALGEBRA_API SparseMatrixProduct >::computeRegularProduct(); - -#if !defined(SOFA_LINEARAGEBRA_SPARSEMATRIXPRODUCT_EIGENSPARSEMATRIX_CPP) - extern template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; - extern template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; - - extern template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; - extern template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; -#endif - -} //namespace sofa::linearalgebra \ No newline at end of file diff --git a/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp b/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp index e10a1eaff6b..78579e5a414 100644 --- a/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp +++ b/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp @@ -19,19 +19,21 @@ * * * Contact information: contact@sofa-framework.org * ******************************************************************************/ -#include -#include +#include #include #include #include #include #include +#include -template + +template struct TestSparseMatrixProductTraits { - using Matrix = TMatrix; + using LHSMatrix = TLHSMatrix; + using RHSMatrix = TRHSMatrix; using Real = TReal; }; @@ -43,34 +45,37 @@ struct TestSparseMatrixProductTraits template struct TestSparseMatrixProduct : public sofa::testing::SparseMatrixTest { - using Matrix = typename T::Matrix; + using LHSMatrix = typename T::LHSMatrix; + using RHSMatrix = typename T::RHSMatrix; using Real = typename T::Real; using Base = sofa::testing::SparseMatrixTest; using Base::generateRandomSparseMatrix; using Base::copyFromEigen; using Base::compareSparseMatrix; - bool checkMatrix(typename Matrix::Index nbRowsA, typename Matrix::Index nbColsA, typename Matrix::Index nbColsB, Real sparsity) + bool checkMatrix(typename LHSMatrix::Index nbRowsA, typename LHSMatrix::Index nbColsA, typename RHSMatrix::Index nbColsB, Real sparsity) { - Eigen::SparseMatrix eigen_a, eigen_b; + Eigen::SparseMatrix eigen_a; + Eigen::SparseMatrix eigen_b; generateRandomSparseMatrix(eigen_a, nbRowsA, nbColsA, sparsity); generateRandomSparseMatrix(eigen_b, nbColsA, nbColsB, sparsity); - Matrix A, B; + LHSMatrix A; + RHSMatrix B; copyFromEigen(A, eigen_a); copyFromEigen(B, eigen_b); EXPECT_GT(eigen_a.outerSize(), 0); EXPECT_GT(eigen_b.outerSize(), 0); - Eigen::SparseMatrix eigen_c = eigen_a * eigen_b; + Eigen::SparseMatrix eigen_c = eigen_a * eigen_b; EXPECT_EQ(eigen_c.rows(), nbRowsA); EXPECT_EQ(eigen_c.cols(), nbColsB); EXPECT_GT(eigen_c.outerSize(), 0); //to make sure that there are non-zero values in the result matrix - sofa::linearalgebra::SparseMatrixProduct product(&A, &B); + sofa::linearalgebra::SparseMatrixProduct > product(&A, &B); product.computeProduct(); EXPECT_TRUE(compareSparseMatrix(eigen_c, product.getProductResult())); @@ -92,7 +97,7 @@ struct TestSparseMatrixProduct : public sofa::testing::SparseMatrixTest; using TestSparseMatrixProductImplementations = ::testing::Types< - TestSparseMatrixProductTraits, float>, - TestSparseMatrixProductTraits, double>, - TestSparseMatrixProductTraits, float>, - TestSparseMatrixProductTraits, double> - // TestSparseMatrixProductTraits > + TestSparseMatrixProductTraits, Eigen::SparseMatrix, float>, + TestSparseMatrixProductTraits, Eigen::SparseMatrix, double>, + TestSparseMatrixProductTraits, Eigen::SparseMatrix, float>, + TestSparseMatrixProductTraits, Eigen::SparseMatrix, double>, + TestSparseMatrixProductTraits, Eigen::SparseMatrix, float>, + TestSparseMatrixProductTraits, Eigen::SparseMatrix, double>, + TestSparseMatrixProductTraits, Eigen::SparseMatrix, float>, + TestSparseMatrixProductTraits, Eigen::SparseMatrix, double> >; TYPED_TEST_SUITE(TestSparseMatrixProduct, TestSparseMatrixProductImplementations); @@ -134,4 +142,4 @@ TYPED_TEST(TestSparseMatrixProduct, rectangularMatrix ) ASSERT_TRUE( this->checkMatrix( 1000, 3000, 2000, 20. / 1000. ) ); ASSERT_TRUE( this->checkMatrix( 20, 30, 10, 1. ) ); -} \ No newline at end of file +} From 1ca9dc3b88580076ce3dc97be74d36a0d6206f0f Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Tue, 27 Feb 2024 12:37:57 +0100 Subject: [PATCH 2/6] Move tests as a template test --- .../LinearAlgebra/Testing/CMakeLists.txt | 1 + .../SparseMatrixProduct_test.h | 139 ++++++++++++++++++ .../sofa/linearalgebra/SparseMatrixProduct.h | 7 +- .../test/SparseMatrixProduct_test.cpp | 114 +------------- 4 files changed, 152 insertions(+), 109 deletions(-) create mode 100644 Sofa/framework/LinearAlgebra/Testing/src/Sofa.LinearAlgebra.Testing/SparseMatrixProduct_test.h diff --git a/Sofa/framework/LinearAlgebra/Testing/CMakeLists.txt b/Sofa/framework/LinearAlgebra/Testing/CMakeLists.txt index 64dc968e5d8..f2f4f157164 100644 --- a/Sofa/framework/LinearAlgebra/Testing/CMakeLists.txt +++ b/Sofa/framework/LinearAlgebra/Testing/CMakeLists.txt @@ -4,6 +4,7 @@ project(Sofa.LinearAlgebra.Testing LANGUAGES CXX) set(HEADER_FILES src/Sofa.LinearAlgebra.Testing/BaseMatrix_test.h src/Sofa.LinearAlgebra.Testing/SparseMatrixTest.h + src/Sofa.LinearAlgebra.Testing/SparseMatrixProduct_test.h ) add_library(${PROJECT_NAME} INTERFACE) diff --git a/Sofa/framework/LinearAlgebra/Testing/src/Sofa.LinearAlgebra.Testing/SparseMatrixProduct_test.h b/Sofa/framework/LinearAlgebra/Testing/src/Sofa.LinearAlgebra.Testing/SparseMatrixProduct_test.h new file mode 100644 index 00000000000..6d50153cd9c --- /dev/null +++ b/Sofa/framework/LinearAlgebra/Testing/src/Sofa.LinearAlgebra.Testing/SparseMatrixProduct_test.h @@ -0,0 +1,139 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include +#include + +#include +#include +#include + +namespace sofa::linearalgebra::testing +{ +template +struct TestSparseMatrixProductTraits +{ + using LHSMatrix = TLHSMatrix; + using RHSMatrix = TRHSMatrix; + using Real = TReal; +}; + +/** + * Test the class SparseMatrixProduct + * The class is designed to use the templates defined in TestSparseMatrixProductTraits + * The type of matrix can be any of the types supported by SparseMatrixProduct. + */ +template +struct TestSparseMatrixProduct : public sofa::testing::SparseMatrixTest +{ + using LHSMatrix = typename T::LHSMatrix; + using RHSMatrix = typename T::RHSMatrix; + using Real = typename T::Real; + using Base = sofa::testing::SparseMatrixTest; + using Base::generateRandomSparseMatrix; + using Base::copyFromEigen; + using Base::compareSparseMatrix; + + bool checkMatrix(typename LHSMatrix::Index nbRowsA, typename LHSMatrix::Index nbColsA, typename RHSMatrix::Index nbColsB, Real sparsity) + { + Eigen::SparseMatrix eigen_a; + Eigen::SparseMatrix eigen_b; + + generateRandomSparseMatrix(eigen_a, nbRowsA, nbColsA, sparsity); + generateRandomSparseMatrix(eigen_b, nbColsA, nbColsB, sparsity); + + LHSMatrix A; + RHSMatrix B; + copyFromEigen(A, eigen_a); + copyFromEigen(B, eigen_b); + + EXPECT_GT(eigen_a.outerSize(), 0); + EXPECT_GT(eigen_b.outerSize(), 0); + + Eigen::SparseMatrix eigen_c = eigen_a * eigen_b; + + EXPECT_EQ(eigen_c.rows(), nbRowsA); + EXPECT_EQ(eigen_c.cols(), nbColsB); + EXPECT_GT(eigen_c.outerSize(), 0); //to make sure that there are non-zero values in the result matrix + + sofa::linearalgebra::SparseMatrixProduct > product(&A, &B); + product.computeProduct(); + + EXPECT_TRUE(compareSparseMatrix(eigen_c, product.getProductResult())); + + //the second time computeProduct is called uses the faster algorithm + product.computeProduct(); + EXPECT_TRUE(compareSparseMatrix(eigen_c, product.getProductResult())); + + // force re-computing the intersection + product.computeProduct(true); + EXPECT_TRUE(compareSparseMatrix(eigen_c, product.getProductResult())); + + //modify the values of A, but not its pattern + for (int i = 0; i < eigen_a.nonZeros(); ++i) + { + eigen_a.valuePtr()[i] = static_cast(sofa::helper::drand(1)); + } + + eigen_c = eigen_a * eigen_b; //result is updated using the regular matrix product + copyFromEigen(A, eigen_a); + + product.m_lhs = &A; + product.computeProduct(); //intersection is already computed: uses the faster algorithm + EXPECT_TRUE(compareSparseMatrix(eigen_c, product.getProductResult())); + + return true; + } +}; + +TYPED_TEST_SUITE_P(TestSparseMatrixProduct); + +TYPED_TEST_P(TestSparseMatrixProduct, squareMatrix ) +{ + ASSERT_TRUE( this->checkMatrix( 5, 5, 5, 1. / 5. ) ); + ASSERT_TRUE( this->checkMatrix( 5, 5, 5, 3. / 5. ) ); + + ASSERT_TRUE( this->checkMatrix( 100, 1000, 1000, 1. / 1000. ) ); + ASSERT_TRUE( this->checkMatrix( 1000, 1000, 1000, 20. / 1000. ) ); + + ASSERT_TRUE( this->checkMatrix( 20, 20, 20, 1. ) ); +} + +TYPED_TEST_P(TestSparseMatrixProduct, rectangularMatrix ) +{ + ASSERT_TRUE( this->checkMatrix( 5, 10, 7, 1. / 5. ) ); + ASSERT_TRUE( this->checkMatrix( 5, 10, 7, 3. / 5. ) ); + + ASSERT_TRUE( this->checkMatrix( 10, 5, 7, 1. / 5. ) ); + ASSERT_TRUE( this->checkMatrix( 10, 5, 7, 3. / 5. ) ); + + ASSERT_TRUE( this->checkMatrix( 1000, 3000, 2000, 1. / 1000. ) ); + ASSERT_TRUE( this->checkMatrix( 1000, 3000, 2000, 20. / 1000. ) ); + + ASSERT_TRUE( this->checkMatrix( 20, 30, 10, 1. ) ); +} + +REGISTER_TYPED_TEST_SUITE_P(TestSparseMatrixProduct, + squareMatrix,rectangularMatrix); + +} diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.h index e9b86afb8e8..e65aaf023c4 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.h @@ -74,8 +74,9 @@ class SparseMatrixProduct void invalidateIntersection(); - SparseMatrixProduct(Lhs* a, Rhs* b) : m_lhs(a), m_rhs(b) {} + SparseMatrixProduct(Lhs* lhs, Rhs* rhs) : m_lhs(lhs), m_rhs(rhs) {} SparseMatrixProduct() = default; + virtual ~SparseMatrixProduct() = default; struct Intersection { @@ -91,11 +92,11 @@ class SparseMatrixProduct }; protected: - ProductResult m_productResult; /// Result of A*B + ProductResult m_productResult; /// Result of LHS * RHS bool m_hasComputedIntersection { false }; void computeIntersection(); - void computeProductFromIntersection(); + virtual void computeProductFromIntersection(); Intersection m_intersectionAB; diff --git a/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp b/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp index 78579e5a414..6a4804e1481 100644 --- a/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp +++ b/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp @@ -19,93 +19,13 @@ * * * Contact information: contact@sofa-framework.org * ******************************************************************************/ +#include #include -#include -#include -#include -#include -#include -#include - - -template -struct TestSparseMatrixProductTraits -{ - using LHSMatrix = TLHSMatrix; - using RHSMatrix = TRHSMatrix; - using Real = TReal; -}; - -/** - * Test the class SparseMatrixProduct - * The class is designed to use the templates defined in TestSparseMatrixProductTraits - * The type of matrix can be any of the types supported by SparseMatrixProduct. - */ -template -struct TestSparseMatrixProduct : public sofa::testing::SparseMatrixTest +namespace sofa { - using LHSMatrix = typename T::LHSMatrix; - using RHSMatrix = typename T::RHSMatrix; - using Real = typename T::Real; - using Base = sofa::testing::SparseMatrixTest; - using Base::generateRandomSparseMatrix; - using Base::copyFromEigen; - using Base::compareSparseMatrix; - - bool checkMatrix(typename LHSMatrix::Index nbRowsA, typename LHSMatrix::Index nbColsA, typename RHSMatrix::Index nbColsB, Real sparsity) - { - Eigen::SparseMatrix eigen_a; - Eigen::SparseMatrix eigen_b; - - generateRandomSparseMatrix(eigen_a, nbRowsA, nbColsA, sparsity); - generateRandomSparseMatrix(eigen_b, nbColsA, nbColsB, sparsity); - - LHSMatrix A; - RHSMatrix B; - copyFromEigen(A, eigen_a); - copyFromEigen(B, eigen_b); - - EXPECT_GT(eigen_a.outerSize(), 0); - EXPECT_GT(eigen_b.outerSize(), 0); - - Eigen::SparseMatrix eigen_c = eigen_a * eigen_b; - - EXPECT_EQ(eigen_c.rows(), nbRowsA); - EXPECT_EQ(eigen_c.cols(), nbColsB); - EXPECT_GT(eigen_c.outerSize(), 0); //to make sure that there are non-zero values in the result matrix - - sofa::linearalgebra::SparseMatrixProduct > product(&A, &B); - product.computeProduct(); - - EXPECT_TRUE(compareSparseMatrix(eigen_c, product.getProductResult())); - - //the second time computeProduct is called uses the faster algorithm - product.computeProduct(); - EXPECT_TRUE(compareSparseMatrix(eigen_c, product.getProductResult())); - // force re-computing the intersection - product.computeProduct(true); - EXPECT_TRUE(compareSparseMatrix(eigen_c, product.getProductResult())); - - //modify the values of A, but not its pattern - for (int i = 0; i < eigen_a.nonZeros(); ++i) - { - eigen_a.valuePtr()[i] = static_cast(sofa::helper::drand(1)); - } - - eigen_c = eigen_a * eigen_b; //result is updated using the regular matrix product - copyFromEigen(A, eigen_a); - - product.m_lhs = &A; - product.computeProduct(); //intersection is already computed: uses the faster algorithm - EXPECT_TRUE(compareSparseMatrix(eigen_c, product.getProductResult())); - - return true; - } -}; - -using CRSMatrixScalar = sofa::linearalgebra::CompressedRowSparseMatrix; +using namespace sofa::linearalgebra::testing; using TestSparseMatrixProductImplementations = ::testing::Types< TestSparseMatrixProductTraits, Eigen::SparseMatrix, float>, @@ -117,29 +37,11 @@ using TestSparseMatrixProductImplementations = ::testing::Types< TestSparseMatrixProductTraits, Eigen::SparseMatrix, float>, TestSparseMatrixProductTraits, Eigen::SparseMatrix, double> >; -TYPED_TEST_SUITE(TestSparseMatrixProduct, TestSparseMatrixProductImplementations); - -TYPED_TEST(TestSparseMatrixProduct, squareMatrix ) -{ - ASSERT_TRUE( this->checkMatrix( 5, 5, 5, 1. / 5. ) ); - ASSERT_TRUE( this->checkMatrix( 5, 5, 5, 3. / 5. ) ); - - ASSERT_TRUE( this->checkMatrix( 100, 1000, 1000, 1. / 1000. ) ); - ASSERT_TRUE( this->checkMatrix( 1000, 1000, 1000, 20. / 1000. ) ); - - ASSERT_TRUE( this->checkMatrix( 20, 20, 20, 1. ) ); -} - -TYPED_TEST(TestSparseMatrixProduct, rectangularMatrix ) -{ - ASSERT_TRUE( this->checkMatrix( 5, 10, 7, 1. / 5. ) ); - ASSERT_TRUE( this->checkMatrix( 5, 10, 7, 3. / 5. ) ); - - ASSERT_TRUE( this->checkMatrix( 10, 5, 7, 1. / 5. ) ); - ASSERT_TRUE( this->checkMatrix( 10, 5, 7, 3. / 5. ) ); - ASSERT_TRUE( this->checkMatrix( 1000, 3000, 2000, 1. / 1000. ) ); - ASSERT_TRUE( this->checkMatrix( 1000, 3000, 2000, 20. / 1000. ) ); +INSTANTIATE_TYPED_TEST_SUITE_P( + TestSparseMatrixProduct, + TestSparseMatrixProduct, + TestSparseMatrixProductImplementations +); - ASSERT_TRUE( this->checkMatrix( 20, 30, 10, 1. ) ); } From 9e725ec58ea4ed3c1ca84c7aecf446179be27c8f Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 29 Feb 2024 16:28:36 +0100 Subject: [PATCH 3/6] Fix all combinations of storage --- .../SparseMatrixProduct_test.h | 57 +++++++++++-------- .../linearalgebra/SparseMatrixProduct.inl | 29 +++++++--- .../test/SparseMatrixProduct_test.cpp | 29 +++++++--- 3 files changed, 75 insertions(+), 40 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/Testing/src/Sofa.LinearAlgebra.Testing/SparseMatrixProduct_test.h b/Sofa/framework/LinearAlgebra/Testing/src/Sofa.LinearAlgebra.Testing/SparseMatrixProduct_test.h index 6d50153cd9c..5f714d71037 100644 --- a/Sofa/framework/LinearAlgebra/Testing/src/Sofa.LinearAlgebra.Testing/SparseMatrixProduct_test.h +++ b/Sofa/framework/LinearAlgebra/Testing/src/Sofa.LinearAlgebra.Testing/SparseMatrixProduct_test.h @@ -30,12 +30,19 @@ namespace sofa::linearalgebra::testing { -template -struct TestSparseMatrixProductTraits + +template +struct SparseMatrixProductInit { - using LHSMatrix = TLHSMatrix; - using RHSMatrix = TRHSMatrix; - using Real = TReal; + static void init(T& product) + { + SOFA_UNUSED(product); + }; + + static void cleanup(T& product) + { + SOFA_UNUSED(product); + } }; /** @@ -44,12 +51,12 @@ struct TestSparseMatrixProductTraits * The type of matrix can be any of the types supported by SparseMatrixProduct. */ template -struct TestSparseMatrixProduct : public sofa::testing::SparseMatrixTest +struct TestSparseMatrixProduct : public sofa::testing::SparseMatrixTest { - using LHSMatrix = typename T::LHSMatrix; - using RHSMatrix = typename T::RHSMatrix; - using Real = typename T::Real; - using Base = sofa::testing::SparseMatrixTest; + using LHSMatrix = typename T::LhsCleaned; + using RHSMatrix = typename T::RhsCleaned; + using Real = typename T::LhsScalar; + using Base = sofa::testing::SparseMatrixTest; using Base::generateRandomSparseMatrix; using Base::copyFromEigen; using Base::compareSparseMatrix; @@ -76,7 +83,8 @@ struct TestSparseMatrixProduct : public sofa::testing::SparseMatrixTest > product(&A, &B); + T product(&A, &B); + SparseMatrixProductInit::init(product); product.computeProduct(); EXPECT_TRUE(compareSparseMatrix(eigen_c, product.getProductResult())); @@ -102,6 +110,8 @@ struct TestSparseMatrixProduct : public sofa::testing::SparseMatrixTest::cleanup(product); + return true; } }; @@ -110,27 +120,28 @@ TYPED_TEST_SUITE_P(TestSparseMatrixProduct); TYPED_TEST_P(TestSparseMatrixProduct, squareMatrix ) { - ASSERT_TRUE( this->checkMatrix( 5, 5, 5, 1. / 5. ) ); - ASSERT_TRUE( this->checkMatrix( 5, 5, 5, 3. / 5. ) ); + EXPECT_TRUE( this->checkMatrix( 5, 5, 5, 1. / 5. ) ); + EXPECT_TRUE( this->checkMatrix( 5, 5, 5, 3. / 5. ) ); - ASSERT_TRUE( this->checkMatrix( 100, 1000, 1000, 1. / 1000. ) ); - ASSERT_TRUE( this->checkMatrix( 1000, 1000, 1000, 20. / 1000. ) ); + EXPECT_TRUE( this->checkMatrix( 100, 1000, 1000, 1. / 1000. ) ); + EXPECT_TRUE( this->checkMatrix( 1000, 1000, 1000, 20. / 1000. ) ); + // EXPECT_TRUE( this->checkMatrix( 1000, 1000, 1000, 150. / 1000. ) ); - ASSERT_TRUE( this->checkMatrix( 20, 20, 20, 1. ) ); + EXPECT_TRUE( this->checkMatrix( 20, 20, 20, 1. ) ); } TYPED_TEST_P(TestSparseMatrixProduct, rectangularMatrix ) { - ASSERT_TRUE( this->checkMatrix( 5, 10, 7, 1. / 5. ) ); - ASSERT_TRUE( this->checkMatrix( 5, 10, 7, 3. / 5. ) ); + EXPECT_TRUE( this->checkMatrix( 5, 10, 7, 1. / 5. ) ); + EXPECT_TRUE( this->checkMatrix( 5, 10, 7, 3. / 5. ) ); - ASSERT_TRUE( this->checkMatrix( 10, 5, 7, 1. / 5. ) ); - ASSERT_TRUE( this->checkMatrix( 10, 5, 7, 3. / 5. ) ); + EXPECT_TRUE( this->checkMatrix( 10, 5, 7, 1. / 5. ) ); + EXPECT_TRUE( this->checkMatrix( 10, 5, 7, 3. / 5. ) ); - ASSERT_TRUE( this->checkMatrix( 1000, 3000, 2000, 1. / 1000. ) ); - ASSERT_TRUE( this->checkMatrix( 1000, 3000, 2000, 20. / 1000. ) ); + EXPECT_TRUE( this->checkMatrix( 1000, 3000, 2000, 1. / 1000. ) ); + EXPECT_TRUE( this->checkMatrix( 1000, 3000, 2000, 20. / 1000. ) ); - ASSERT_TRUE( this->checkMatrix( 20, 30, 10, 1. ) ); + EXPECT_TRUE( this->checkMatrix( 20, 30, 10, 1. ) ); } REGISTER_TYPED_TEST_SUITE_P(TestSparseMatrixProduct, diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.inl b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.inl index 070619075d0..ac01c957822 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.inl +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.inl @@ -23,8 +23,6 @@ #include #include #include -#include -#include namespace sofa::linearalgebra::sparsematrixproduct @@ -226,13 +224,26 @@ void SparseMatrixProduct::computeIntersection() m_intersectionAB.intersection.push_back(product.valuePtr()[i].getIndices()); //depending on the storage scheme, Eigen can change the order of the lhs and rhs - if constexpr (Lhs::IsRowMajor || Rhs::IsRowMajor) + //Note: the condition has been determined empirically, using unit tests + //testing all possible combinations = 2^3 = 8 + if constexpr ((Lhs::IsRowMajor && Rhs::IsRowMajor && ResultType::IsRowMajor) + || ((Lhs::IsRowMajor || Rhs::IsRowMajor) && !ResultType::IsRowMajor)) { for (auto& [lhsIndex, rhsIndex] : m_intersectionAB.intersection.back()) { std::swap(lhsIndex, rhsIndex); } } + +#if !defined(NDEBUG) + const auto lhsNonZeros = m_lhs->nonZeros(); + const auto rhsNonZeros = m_rhs->nonZeros(); + for (const auto& [lhsIndex, rhsIndex] : m_intersectionAB.intersection.back()) + { + assert(lhsIndex < lhsNonZeros); + assert(rhsIndex < rhsNonZeros); + } +#endif } m_productResult = product.template cast(); @@ -248,12 +259,17 @@ void __computeProductFromIntersection(const Lhs* lhs, const Rhs* rhs, ResultType auto* rhs_ptr = rhs->valuePtr(); auto* product_ptr = product->valuePtr(); + const auto lhsNonZeros = lhs->nonZeros(); + const auto rhsNonZeros = rhs->nonZeros(); + for (const auto& pairs : intersection.intersection) { auto& value = *product_ptr++; value = 0; - for (auto [lhsIndex, rhsIndex] : pairs) + for (const auto& [lhsIndex, rhsIndex] : pairs) { + assert(lhsIndex < lhsNonZeros); + assert(rhsIndex < rhsNonZeros); value += lhs_ptr[lhsIndex] * rhs_ptr[rhsIndex]; } } @@ -262,11 +278,6 @@ void __computeProductFromIntersection(const Lhs* lhs, const Rhs* rhs, ResultType template void SparseMatrixProduct::computeProductFromIntersection() { - if (m_hasComputedIntersection == false) - { - msg_error("SparseMatrixProduct") << "Intersection computation is required before computing the product"; - return; - } __computeProductFromIntersection(m_lhs, m_rhs, &m_productResult, m_intersectionAB); } diff --git a/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp b/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp index 6a4804e1481..2c2e208325e 100644 --- a/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp +++ b/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp @@ -27,17 +27,30 @@ namespace sofa using namespace sofa::linearalgebra::testing; +#define DEFINE_TEST_FOR_TYPE(scalar, StorageLHS, StorageRHS, StorageResult)\ + sofa::linearalgebra::SparseMatrixProduct<\ + Eigen::SparseMatrix,\ + Eigen::SparseMatrix,\ + Eigen::SparseMatrix\ + > +#define DEFINE_TEST_FOR_STORAGE(StorageLHS, StorageRHS, StorageResult)\ + DEFINE_TEST_FOR_TYPE(float, StorageLHS, StorageRHS, StorageResult),\ + DEFINE_TEST_FOR_TYPE(double, StorageLHS, StorageRHS, StorageResult) + using TestSparseMatrixProductImplementations = ::testing::Types< - TestSparseMatrixProductTraits, Eigen::SparseMatrix, float>, - TestSparseMatrixProductTraits, Eigen::SparseMatrix, double>, - TestSparseMatrixProductTraits, Eigen::SparseMatrix, float>, - TestSparseMatrixProductTraits, Eigen::SparseMatrix, double>, - TestSparseMatrixProductTraits, Eigen::SparseMatrix, float>, - TestSparseMatrixProductTraits, Eigen::SparseMatrix, double>, - TestSparseMatrixProductTraits, Eigen::SparseMatrix, float>, - TestSparseMatrixProductTraits, Eigen::SparseMatrix, double> + DEFINE_TEST_FOR_STORAGE(Eigen::ColMajor, Eigen::ColMajor, Eigen::ColMajor), + DEFINE_TEST_FOR_STORAGE(Eigen::RowMajor, Eigen::ColMajor, Eigen::ColMajor), + DEFINE_TEST_FOR_STORAGE(Eigen::ColMajor, Eigen::RowMajor, Eigen::ColMajor), + DEFINE_TEST_FOR_STORAGE(Eigen::RowMajor, Eigen::RowMajor, Eigen::ColMajor), + DEFINE_TEST_FOR_STORAGE(Eigen::ColMajor, Eigen::ColMajor, Eigen::RowMajor), + DEFINE_TEST_FOR_STORAGE(Eigen::RowMajor, Eigen::ColMajor, Eigen::RowMajor), + DEFINE_TEST_FOR_STORAGE(Eigen::ColMajor, Eigen::RowMajor, Eigen::RowMajor), + DEFINE_TEST_FOR_STORAGE(Eigen::RowMajor, Eigen::RowMajor, Eigen::RowMajor) >; +#undef DEFINE_TEST_FOR_STORAGE +#undef DEFINE_TEST_FOR_TYPE + INSTANTIATE_TYPED_TEST_SUITE_P( TestSparseMatrixProduct, TestSparseMatrixProduct, From 787c14154c200ac214d82659ca6b8ce17ee90950 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 29 Feb 2024 16:28:49 +0100 Subject: [PATCH 4/6] Introduce parallel computation --- Sofa/framework/Simulation/Core/CMakeLists.txt | 1 + .../simulation/ParallelSparseMatrixProduct.h | 68 ++++++++++++++++ .../Simulation/Core/test/CMakeLists.txt | 3 +- .../test/ParallelSparseMatrixProduct_test.cpp | 78 +++++++++++++++++++ 4 files changed, 149 insertions(+), 1 deletion(-) create mode 100644 Sofa/framework/Simulation/Core/src/sofa/simulation/ParallelSparseMatrixProduct.h create mode 100644 Sofa/framework/Simulation/Core/test/ParallelSparseMatrixProduct_test.cpp diff --git a/Sofa/framework/Simulation/Core/CMakeLists.txt b/Sofa/framework/Simulation/Core/CMakeLists.txt index c0d6c8d5da7..3e3f94b3980 100644 --- a/Sofa/framework/Simulation/Core/CMakeLists.txt +++ b/Sofa/framework/Simulation/Core/CMakeLists.txt @@ -38,6 +38,7 @@ set(HEADER_FILES ${SRC_ROOT}/Node.h ${SRC_ROOT}/Node.inl ${SRC_ROOT}/ParallelForEach.h + ${SRC_ROOT}/ParallelSparseMatrixProduct.h ${SRC_ROOT}/ParallelVisitorScheduler.h ${SRC_ROOT}/PauseEvent.h ${SRC_ROOT}/PipelineImpl.h diff --git a/Sofa/framework/Simulation/Core/src/sofa/simulation/ParallelSparseMatrixProduct.h b/Sofa/framework/Simulation/Core/src/sofa/simulation/ParallelSparseMatrixProduct.h new file mode 100644 index 00000000000..87000e6951c --- /dev/null +++ b/Sofa/framework/Simulation/Core/src/sofa/simulation/ParallelSparseMatrixProduct.h @@ -0,0 +1,68 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include +#include + + +namespace sofa::simulation +{ + +template +class ParallelSparseMatrixProduct + : public linearalgebra::SparseMatrixProduct +{ +public: + using linearalgebra::SparseMatrixProduct::SparseMatrixProduct; + TaskScheduler* taskScheduler { nullptr }; + + void computeProductFromIntersection() override + { + assert(this->m_intersectionAB.intersection.size() == this->m_productResult.nonZeros()); + assert(taskScheduler); + + auto* lhs_ptr = this->m_lhs->valuePtr(); + auto* rhs_ptr = this->m_rhs->valuePtr(); + auto* product_ptr = this->m_productResult.valuePtr(); + + parallelForEachRange(*taskScheduler, + this->m_intersectionAB.intersection.begin(), this->m_intersectionAB.intersection.end(), + [lhs_ptr, rhs_ptr, product_ptr, this](const auto& range) + { + auto i = std::distance(this->m_intersectionAB.intersection.begin(), range.start); + auto* p = product_ptr + i; + + for (auto it = range.start; it != range.end; ++it) + { + auto& value = *p++; + value = 0; + for (const auto& [lhsIndex, rhsIndex] : *it) + { + value += lhs_ptr[lhsIndex] * rhs_ptr[rhsIndex]; + } + } + }); + } +}; + +} diff --git a/Sofa/framework/Simulation/Core/test/CMakeLists.txt b/Sofa/framework/Simulation/Core/test/CMakeLists.txt index 09724fd448b..fe9d7599be8 100644 --- a/Sofa/framework/Simulation/Core/test/CMakeLists.txt +++ b/Sofa/framework/Simulation/Core/test/CMakeLists.txt @@ -7,6 +7,7 @@ set(SOURCE_FILES RequiredPlugin_test.cpp SceneCheckRegistry_test.cpp Simulation_test.cpp + ParallelSparseMatrixProduct_test.cpp TaskSchedulerFactory_test.cpp TaskSchedulerTestTasks.cpp TaskSchedulerTestTasks.h @@ -14,6 +15,6 @@ set(SOURCE_FILES ) add_executable(${PROJECT_NAME} ${SOURCE_FILES}) -target_link_libraries(${PROJECT_NAME} Sofa.Testing Sofa.Simulation.Core) +target_link_libraries(${PROJECT_NAME} Sofa.Testing Sofa.Simulation.Core Sofa.LinearAlgebra.Testing) add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME}) diff --git a/Sofa/framework/Simulation/Core/test/ParallelSparseMatrixProduct_test.cpp b/Sofa/framework/Simulation/Core/test/ParallelSparseMatrixProduct_test.cpp new file mode 100644 index 00000000000..3b8c2de2e42 --- /dev/null +++ b/Sofa/framework/Simulation/Core/test/ParallelSparseMatrixProduct_test.cpp @@ -0,0 +1,78 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#include +#include +#include +#include + +namespace sofa +{ + +using namespace sofa::linearalgebra::testing; + +template +struct SparseMatrixProductInit< + sofa::simulation::ParallelSparseMatrixProduct> +{ + static void init(sofa::simulation::ParallelSparseMatrixProduct& product) + { + product.taskScheduler = simulation::MainTaskSchedulerFactory::createInRegistry(); + product.taskScheduler->init(); + }; + + static void cleanup(sofa::simulation::ParallelSparseMatrixProduct& product) + { + // simulation::MainTaskSchedulerRegistry::clear(); + } +}; + +#define DEFINE_TEST_FOR_TYPE(scalar, StorageLHS, StorageRHS, StorageResult)\ + sofa::simulation::ParallelSparseMatrixProduct<\ + Eigen::SparseMatrix,\ + Eigen::SparseMatrix,\ + Eigen::SparseMatrix\ + > +#define DEFINE_TEST_FOR_STORAGE(StorageLHS, StorageRHS, StorageResult)\ + DEFINE_TEST_FOR_TYPE(float, StorageLHS, StorageRHS, StorageResult),\ + DEFINE_TEST_FOR_TYPE(double, StorageLHS, StorageRHS, StorageResult) + +using TestSparseMatrixProductImplementations = ::testing::Types< + DEFINE_TEST_FOR_STORAGE(Eigen::ColMajor, Eigen::ColMajor, Eigen::ColMajor), + DEFINE_TEST_FOR_STORAGE(Eigen::RowMajor, Eigen::ColMajor, Eigen::ColMajor), + DEFINE_TEST_FOR_STORAGE(Eigen::ColMajor, Eigen::RowMajor, Eigen::ColMajor), + DEFINE_TEST_FOR_STORAGE(Eigen::RowMajor, Eigen::RowMajor, Eigen::ColMajor), + DEFINE_TEST_FOR_STORAGE(Eigen::ColMajor, Eigen::ColMajor, Eigen::RowMajor), + DEFINE_TEST_FOR_STORAGE(Eigen::RowMajor, Eigen::ColMajor, Eigen::RowMajor), + DEFINE_TEST_FOR_STORAGE(Eigen::ColMajor, Eigen::RowMajor, Eigen::RowMajor), + DEFINE_TEST_FOR_STORAGE(Eigen::RowMajor, Eigen::RowMajor, Eigen::RowMajor) +>; + +#undef DEFINE_TEST_FOR_STORAGE +#undef DEFINE_TEST_FOR_TYPE + +INSTANTIATE_TYPED_TEST_SUITE_P( + TestParallelSparseMatrixProduct, + TestSparseMatrixProduct, + TestSparseMatrixProductImplementations +); + +} From 64f4ff2a326f67ae4baae50c42dca54137d04df6 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 1 Mar 2024 14:05:29 +0100 Subject: [PATCH 5/6] Manage options for more types --- .../linearalgebra/SparseMatrixProduct.inl | 58 ++++++++++++++----- 1 file changed, 42 insertions(+), 16 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.inl b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.inl index ac01c957822..98cbc623c54 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.inl +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.inl @@ -184,18 +184,51 @@ void flagValueIndices(Eigen::SparseMatrix +struct EigenOptions +{ + static constexpr auto value = T::Options; +}; + +template +static constexpr auto EigenOptions_v = EigenOptions::value; + +template +struct EigenOptions> +{ + static constexpr auto value = EigenOptions_v; +}; + +template +struct EigenOptions> +{ + static constexpr auto value = EigenOptions_v; +}; + +template +struct EigenOptions> +{ + static constexpr auto value = (EigenOptions_v == Eigen::RowMajor) ? Eigen::ColMajor : Eigen::RowMajor; +}; + +template +struct EigenOptions> +{ + static constexpr auto value = (EigenOptions_v == Eigen::RowMajor) ? Eigen::ColMajor : Eigen::RowMajor; +}; + template void SparseMatrixProduct::computeIntersection() { using LocalLhs = Eigen::SparseMatrix< sparsematrixproduct::IndexedValue, - Lhs::Options, + EigenOptions_v, typename Lhs::StorageIndex >; using LocalRhs = Eigen::SparseMatrix< sparsematrixproduct::IndexedValue, - Rhs::Options, + EigenOptions_v, typename Rhs::StorageIndex >; @@ -250,19 +283,18 @@ void SparseMatrixProduct::computeIntersection() } template -void __computeProductFromIntersection(const Lhs* lhs, const Rhs* rhs, ResultType* product, - const typename SparseMatrixProduct::Intersection& intersection) +void SparseMatrixProduct::computeProductFromIntersection() { assert(intersection.intersection.size() == product->nonZeros()); - auto* lhs_ptr = lhs->valuePtr(); - auto* rhs_ptr = rhs->valuePtr(); - auto* product_ptr = product->valuePtr(); + auto* lhs_ptr = m_lhs->valuePtr(); + auto* rhs_ptr = m_rhs->valuePtr(); + auto* product_ptr = m_productResult.valuePtr(); - const auto lhsNonZeros = lhs->nonZeros(); - const auto rhsNonZeros = rhs->nonZeros(); + const auto lhsNonZeros = m_lhs->nonZeros(); + const auto rhsNonZeros = m_rhs->nonZeros(); - for (const auto& pairs : intersection.intersection) + for (const auto& pairs : m_intersectionAB.intersection) { auto& value = *product_ptr++; value = 0; @@ -275,12 +307,6 @@ void __computeProductFromIntersection(const Lhs* lhs, const Rhs* rhs, ResultType } } -template -void SparseMatrixProduct::computeProductFromIntersection() -{ - __computeProductFromIntersection(m_lhs, m_rhs, &m_productResult, m_intersectionAB); -} - template void SparseMatrixProduct::invalidateIntersection() { From 663a2b7cb792f65979e1d3ee4f5c49d0847aa6ad Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Fri, 1 Mar 2024 14:06:23 +0100 Subject: [PATCH 6/6] [LinearSystem] Introduce ConstantSparsityProjectionMethod --- Sofa/Component/LinearSystem/CMakeLists.txt | 3 + .../ConstantSparsityProjectionMethod.cpp | 33 +++++ .../ConstantSparsityProjectionMethod.h | 68 ++++++++++ .../ConstantSparsityProjectionMethod.inl | 124 ++++++++++++++++++ 4 files changed, 228 insertions(+) create mode 100644 Sofa/Component/LinearSystem/src/sofa/component/linearsystem/ConstantSparsityProjectionMethod.cpp create mode 100644 Sofa/Component/LinearSystem/src/sofa/component/linearsystem/ConstantSparsityProjectionMethod.h create mode 100644 Sofa/Component/LinearSystem/src/sofa/component/linearsystem/ConstantSparsityProjectionMethod.inl diff --git a/Sofa/Component/LinearSystem/CMakeLists.txt b/Sofa/Component/LinearSystem/CMakeLists.txt index 3fee230b0fd..c9de4c3fa38 100644 --- a/Sofa/Component/LinearSystem/CMakeLists.txt +++ b/Sofa/Component/LinearSystem/CMakeLists.txt @@ -11,6 +11,8 @@ set(HEADER_FILES ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/CompositeLinearSystem.inl ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/ConstantSparsityPatternSystem.h ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/ConstantSparsityPatternSystem.inl + ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/ConstantSparsityProjectionMethod.h + ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/ConstantSparsityProjectionMethod.inl ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/CreateMatrixDispatcher.h ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/LinearSystemData.h ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/MappingGraph.h @@ -40,6 +42,7 @@ set(SOURCE_FILES ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/BaseMatrixProjectionMethod.cpp ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/CompositeLinearSystem.cpp ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/ConstantSparsityPatternSystem.cpp + ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/ConstantSparsityProjectionMethod.cpp ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/MappingGraph.cpp ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/MatrixLinearSystem.cpp ${SOFACOMPONENTLINEARSOLVERLINEARSYSTEM_SOURCE_DIR}/MatrixProjectionMethod.cpp diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/ConstantSparsityProjectionMethod.cpp b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/ConstantSparsityProjectionMethod.cpp new file mode 100644 index 00000000000..346c5d2d0dd --- /dev/null +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/ConstantSparsityProjectionMethod.cpp @@ -0,0 +1,33 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define SOFA_COMPONENT_LINEARSYSTEM_CONSTANTSPARSITYPROJECTIONMETHOD_CPP +#include +#include + +namespace sofa::component::linearsystem +{ +template class SOFA_COMPONENT_LINEARSYSTEM_API ConstantSparsityProjectionMethod >; + +int ConstantSparsityProjectionMethodClass = core::RegisterObject("Matrix mapping computing the matrix projection taking advantage of the constant sparsity pattern") + .add< ConstantSparsityProjectionMethod > >(true) + ; +} diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/ConstantSparsityProjectionMethod.h b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/ConstantSparsityProjectionMethod.h new file mode 100644 index 00000000000..13b656525ab --- /dev/null +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/ConstantSparsityProjectionMethod.h @@ -0,0 +1,68 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include + +namespace sofa::component::linearsystem +{ + +/** + * Matrix prjection method computing the matrix projection taking advantage of the constant sparsity pattern + */ +template +class ConstantSparsityProjectionMethod : public MatrixProjectionMethod +{ +public: + SOFA_CLASS(SOFA_TEMPLATE(ConstantSparsityProjectionMethod, TMatrix), SOFA_TEMPLATE(MatrixProjectionMethod, TMatrix)); + using PairMechanicalStates = typename BaseMatrixProjectionMethod::PairMechanicalStates; + using Block = typename Inherit1::Block; + + ConstantSparsityProjectionMethod(); + ~ConstantSparsityProjectionMethod() override; + + Data d_parallelProduct; + + void init() override; + +protected: + void computeProjection( + const Eigen::Map > KMap, + const sofa::type::fixed_array, 2> J, + Eigen::SparseMatrix& JT_K_J) override; + + using K_Type = Eigen::Map >; + using J_Type = Eigen::Map >; + using KJ_Type = Eigen::SparseMatrix; + using JT_Type = const Eigen::Transpose > >; + using JTKJ_Type = Eigen::SparseMatrix; + + std::unique_ptr > m_matrixProductKJ; + std::unique_ptr > m_matrixProductJTKJ; +}; + +#if !defined(SOFA_COMPONENT_LINEARSYSTEM_CONSTANTSPARSITYPROJECTIONMETHOD_CPP) +extern template class SOFA_COMPONENT_LINEARSYSTEM_API ConstantSparsityProjectionMethod >; +#endif + +} + diff --git a/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/ConstantSparsityProjectionMethod.inl b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/ConstantSparsityProjectionMethod.inl new file mode 100644 index 00000000000..a36766c2a64 --- /dev/null +++ b/Sofa/Component/LinearSystem/src/sofa/component/linearsystem/ConstantSparsityProjectionMethod.inl @@ -0,0 +1,124 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include + +#include +#include +#include +#include + +namespace sofa::component::linearsystem +{ + +template +ConstantSparsityProjectionMethod::ConstantSparsityProjectionMethod() + : d_parallelProduct(initData(&d_parallelProduct, true, "parallelProduct", "Compute the matrix product in parallel")) +{} + +template +ConstantSparsityProjectionMethod:: +~ConstantSparsityProjectionMethod() = default; + +template +void ConstantSparsityProjectionMethod::init() +{ + Inherit1::init(); + + if (d_parallelProduct.getValue()) + { + auto* taskScheduler = simulation::MainTaskSchedulerFactory::createInRegistry(); + taskScheduler->init(); + + auto* matrixPrductKJ = new sofa::simulation::ParallelSparseMatrixProduct< + K_Type, J_Type, KJ_Type>(); + + m_matrixProductKJ = std::unique_ptr>(matrixPrductKJ); + + matrixPrductKJ->taskScheduler = taskScheduler; + + + auto* matrixPrductJTKJ = new sofa::simulation::ParallelSparseMatrixProduct< + JT_Type, KJ_Type, JTKJ_Type>(); + + m_matrixProductJTKJ = std::unique_ptr>(matrixPrductJTKJ); + + matrixPrductJTKJ->taskScheduler = taskScheduler; + } + else + { + m_matrixProductKJ = std::make_unique>(); + + m_matrixProductJTKJ = std::make_unique>(); + } +} + +template +void ConstantSparsityProjectionMethod::computeProjection( + const Eigen::Map> KMap, + const sofa::type::fixed_array, 2> J, + Eigen::SparseMatrix& JT_K_J) +{ + if (J[0] && J[1]) + { + const auto JMap0 = this->makeEigenMap(*J[0]); + const auto JMap1 = this->makeEigenMap(*J[1]); + + m_matrixProductKJ->m_lhs = &KMap; + m_matrixProductKJ->m_rhs = &JMap1; + m_matrixProductKJ->computeProduct(); + + const JT_Type JMap0T = JMap0.transpose(); + m_matrixProductJTKJ->m_lhs = &JMap0T; + m_matrixProductJTKJ->m_rhs = &m_matrixProductKJ->getProductResult(); + m_matrixProductJTKJ->computeProduct(); + + JT_K_J = m_matrixProductJTKJ->getProductResult(); + } + else if (J[0] && !J[1]) + { + const auto JMap0 = this->makeEigenMap(*J[0]); + JT_K_J = JMap0.transpose() * KMap; + } + else if (!J[0] && J[1]) + { + const auto JMap1 = this->makeEigenMap(*J[1]); + + m_matrixProductKJ->m_lhs = &KMap; + m_matrixProductKJ->m_rhs = &JMap1; + + m_matrixProductKJ->computeProduct(); + + JT_K_J = m_matrixProductKJ->getProductResult(); + } + else + { + JT_K_J = KMap; + } +} + + +}