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; + } +} + + +}