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/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..5f714d71037 --- /dev/null +++ b/Sofa/framework/LinearAlgebra/Testing/src/Sofa.LinearAlgebra.Testing/SparseMatrixProduct_test.h @@ -0,0 +1,150 @@ +/****************************************************************************** +* 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 SparseMatrixProductInit +{ + static void init(T& product) + { + SOFA_UNUSED(product); + }; + + static void cleanup(T& product) + { + SOFA_UNUSED(product); + } +}; + +/** + * 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::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; + + 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 + + T product(&A, &B); + SparseMatrixProductInit::init(product); + 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())); + + SparseMatrixProductInit::cleanup(product); + + return true; + } +}; + +TYPED_TEST_SUITE_P(TestSparseMatrixProduct); + +TYPED_TEST_P(TestSparseMatrixProduct, squareMatrix ) +{ + EXPECT_TRUE( this->checkMatrix( 5, 5, 5, 1. / 5. ) ); + EXPECT_TRUE( this->checkMatrix( 5, 5, 5, 3. / 5. ) ); + + 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. ) ); + + EXPECT_TRUE( this->checkMatrix( 20, 20, 20, 1. ) ); +} + +TYPED_TEST_P(TestSparseMatrixProduct, rectangularMatrix ) +{ + EXPECT_TRUE( this->checkMatrix( 5, 10, 7, 1. / 5. ) ); + EXPECT_TRUE( this->checkMatrix( 5, 10, 7, 3. / 5. ) ); + + EXPECT_TRUE( this->checkMatrix( 10, 5, 7, 1. / 5. ) ); + EXPECT_TRUE( this->checkMatrix( 10, 5, 7, 3. / 5. ) ); + + EXPECT_TRUE( this->checkMatrix( 1000, 3000, 2000, 1. / 1000. ) ); + EXPECT_TRUE( this->checkMatrix( 1000, 3000, 2000, 20. / 1000. ) ); + + EXPECT_TRUE( this->checkMatrix( 20, 30, 10, 1. ) ); +} + +REGISTER_TYPED_TEST_SUITE_P(TestSparseMatrixProduct, + squareMatrix,rectangularMatrix); + +} 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..e65aaf023c4 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,81 +44,63 @@ 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* lhs, Rhs* rhs) : m_lhs(lhs), m_rhs(rhs) {} SparseMatrixProduct() = default; + virtual ~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 LHS * RHS bool m_hasComputedIntersection { false }; void computeIntersection(); - void computeProductFromIntersection(); + virtual void computeProductFromIntersection(); Intersection m_intersectionAB; }; -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..98cbc623c54 --- /dev/null +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct.inl @@ -0,0 +1,316 @@ +/****************************************************************************** +* 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::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 +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, + EigenOptions_v, + typename Lhs::StorageIndex + >; + + using LocalRhs = Eigen::SparseMatrix< + sparsematrixproduct::IndexedValue, + EigenOptions_v, + 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 + //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(); +} + +template +void SparseMatrixProduct::computeProductFromIntersection() +{ + assert(intersection.intersection.size() == product->nonZeros()); + + auto* lhs_ptr = m_lhs->valuePtr(); + auto* rhs_ptr = m_rhs->valuePtr(); + auto* product_ptr = m_productResult.valuePtr(); + + const auto lhsNonZeros = m_lhs->nonZeros(); + const auto rhsNonZeros = m_rhs->nonZeros(); + + for (const auto& pairs : m_intersectionAB.intersection) + { + auto& value = *product_ptr++; + value = 0; + for (const auto& [lhsIndex, rhsIndex] : pairs) + { + assert(lhsIndex < lhsNonZeros); + assert(rhsIndex < rhsNonZeros); + value += lhs_ptr[lhsIndex] * rhs_ptr[rhsIndex]; + } + } +} + +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/test/SparseMatrixProduct_test.cpp b/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp index e10a1eaff6b..2c2e208325e 100644 --- a/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp +++ b/Sofa/framework/LinearAlgebra/test/SparseMatrixProduct_test.cpp @@ -19,119 +19,42 @@ * * * Contact information: contact@sofa-framework.org * ******************************************************************************/ -#include -#include -#include -#include +#include +#include -#include -#include -#include - -template -struct TestSparseMatrixProductTraits -{ - using Matrix = TMatrix; - 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 Matrix = typename T::Matrix; - 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) - { - Eigen::SparseMatrix eigen_a, eigen_b; - - generateRandomSparseMatrix(eigen_a, nbRowsA, nbColsA, sparsity); - generateRandomSparseMatrix(eigen_b, nbColsA, nbColsB, sparsity); - - Matrix A, 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())); +using namespace sofa::linearalgebra::testing; - // 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.matrixA = &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; +#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, float>, - TestSparseMatrixProductTraits, double>, - TestSparseMatrixProductTraits, float>, - TestSparseMatrixProductTraits, double> - // TestSparseMatrixProductTraits > + 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) >; -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. ) ); +#undef DEFINE_TEST_FOR_STORAGE +#undef DEFINE_TEST_FOR_TYPE - ASSERT_TRUE( this->checkMatrix( 100, 1000, 1000, 1. / 1000. ) ); - ASSERT_TRUE( this->checkMatrix( 1000, 1000, 1000, 20. / 1000. ) ); +INSTANTIATE_TYPED_TEST_SUITE_P( + TestSparseMatrixProduct, + TestSparseMatrixProduct, + TestSparseMatrixProductImplementations +); - 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. ) ); - - ASSERT_TRUE( this->checkMatrix( 20, 30, 10, 1. ) ); -} \ No newline at end of file diff --git a/Sofa/framework/Simulation/Core/CMakeLists.txt b/Sofa/framework/Simulation/Core/CMakeLists.txt index 5a08deb6aa8..8dbfba77981 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/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[EigenSparseMatrix].h b/Sofa/framework/Simulation/Core/src/sofa/simulation/ParallelSparseMatrixProduct.h similarity index 52% rename from Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[EigenSparseMatrix].h rename to Sofa/framework/Simulation/Core/src/sofa/simulation/ParallelSparseMatrixProduct.h index e9a7d446839..87000e6951c 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrixProduct[EigenSparseMatrix].h +++ b/Sofa/framework/Simulation/Core/src/sofa/simulation/ParallelSparseMatrixProduct.h @@ -20,24 +20,49 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #pragma once +#include +#include +#include -#include -#include -namespace sofa::linearalgebra +namespace sofa::simulation { -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(); +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(); -#if !defined(SOFA_LINEARAGEBRA_SPARSEMATRIXPRODUCT_EIGENSPARSEMATRIX_CPP) - extern template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; - extern template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; + 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; - extern template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; - extern template class SOFA_LINEARALGEBRA_API SparseMatrixProduct >; -#endif + 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]; + } + } + }); + } +}; -} //namespace sofa::linearalgebra \ No newline at end of file +} 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 +); + +}