From 97f4b1477c5afe8bf500468e4d08a09d463bf902 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Tue, 7 Apr 2026 15:57:49 +0900 Subject: [PATCH 01/10] implment EigenCholmodSupernodalLLT --- .../LinearSolver/Direct/CMakeLists.txt | 14 ++++ .../direct/EigenCholmodSupernodalLLT.cpp | 53 +++++++++++++ .../direct/EigenCholmodSupernodalLLT.h | 65 ++++++++++++++++ .../linearsolver/direct/EigenSolverFactory.h | 23 ++++++ .../component/linearsolver/direct/init.cpp | 31 ++++++++ cmake/Modules/FindCHOLMOD.cmake | 75 +++++++++++++++++++ 6 files changed, 261 insertions(+) create mode 100644 Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenCholmodSupernodalLLT.cpp create mode 100644 Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenCholmodSupernodalLLT.h create mode 100644 cmake/Modules/FindCHOLMOD.cmake diff --git a/Sofa/Component/LinearSolver/Direct/CMakeLists.txt b/Sofa/Component/LinearSolver/Direct/CMakeLists.txt index 5ee3287128b..301cabb3c1d 100644 --- a/Sofa/Component/LinearSolver/Direct/CMakeLists.txt +++ b/Sofa/Component/LinearSolver/Direct/CMakeLists.txt @@ -7,6 +7,7 @@ set(HEADER_FILES ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/config.h.in ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/init.h ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/AsyncSparseLDLSolver.h + ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenCholmodSupernodalLLT.h ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/AsyncSparseLDLSolver.inl ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/BTDLinearSolver.h ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/BTDLinearSolver.inl @@ -38,6 +39,7 @@ set(HEADER_FILES set(SOURCE_FILES ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/init.cpp ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/AsyncSparseLDLSolver.cpp + ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenCholmodSupernodalLLT.cpp ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/BTDLinearSolver.cpp ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/CholeskySolver.cpp ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenSimplicialLDLT.cpp @@ -59,11 +61,23 @@ sofa_find_package(Sofa.Simulation.Core REQUIRED) sofa_find_package(Sofa.Component.LinearSolver.Iterative REQUIRED) # Only for MatrixLinearSolver which will move to core sofa_find_package(Sofa.Component.LinearSolver.Ordering REQUIRED) +# Optional CHOLMOD support from SuiteSparse for supernodal Cholesky factorization +sofa_find_package(CHOLMOD QUIET) + add_library(${PROJECT_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES} ${WRAPPER_FILES}) target_link_libraries(${PROJECT_NAME} PUBLIC Sofa.Simulation.Core Sofa.Component.LinearSolver.Iterative Sofa.Component.LinearSolver.Ordering) target_link_libraries(${PROJECT_NAME} PUBLIC Threads::Threads) +if(CHOLMOD_FOUND) + target_compile_definitions(${PROJECT_NAME} PUBLIC SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD) + target_link_libraries(${PROJECT_NAME} PUBLIC ${CHOLMOD_LIBRARIES}) + target_include_directories(${PROJECT_NAME} SYSTEM PUBLIC ${CHOLMOD_INCLUDE_DIRS}) + message(STATUS "CHOLMOD found: EigenCholmodSupernodalLLT solver enabled") +else() + message(STATUS "CHOLMOD not found: EigenCholmodSupernodalLLT solver disabled") +endif() + sofa_create_package_with_targets( PACKAGE_NAME ${PROJECT_NAME} PACKAGE_VERSION ${Sofa_VERSION} diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenCholmodSupernodalLLT.cpp b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenCholmodSupernodalLLT.cpp new file mode 100644 index 00000000000..e0bc2207cbc --- /dev/null +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenCholmodSupernodalLLT.cpp @@ -0,0 +1,53 @@ +/****************************************************************************** +* 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 + +#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD + +#define SOFA_COMPONENT_LINEARSOLVER_DIRECT_EIGENCHOLMODSUPERNODALLLT_CPP +#include +#include +#include + +#include + +namespace sofa::component::linearsolver::direct +{ + +// CHOLMOD only supports double precision +static_assert(std::is_same_v, "EigenCholmodSupernodalLLT requires double precision (SReal must be double)"); + +template class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API EigenCholmodSupernodalLLT< SReal >; +template class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API EigenCholmodSupernodalLLT< sofa::type::Mat<3,3,SReal> >; + +MainCholmodSupernodalLLTFactory::~MainCholmodSupernodalLLTFactory() = default; + +void registerEigenCholmodSupernodalLLT(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(core::ObjectRegistrationData("Direct linear solver using a supernodal sparse LL^T Cholesky factorization from CHOLMOD (SuiteSparse).") + .add< EigenCholmodSupernodalLLT< SReal > >() + .add< EigenCholmodSupernodalLLT< sofa::type::Mat<3, 3, SReal> > >()); +} + +} // namespace sofa::component::linearsolver::direct + +#endif // SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenCholmodSupernodalLLT.h b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenCholmodSupernodalLLT.h new file mode 100644 index 00000000000..46cfd950b6c --- /dev/null +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenCholmodSupernodalLLT.h @@ -0,0 +1,65 @@ +/****************************************************************************** +* 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 + +#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD + +#include +#include + +namespace sofa::component::linearsolver::direct +{ + +/** + * Direct linear solver based on a supernodal sparse LL^T Cholesky factorization + * using CHOLMOD from the SuiteSparse library (via Eigen's CholmodSupport). + * + * This solver is significantly faster than simplicial solvers (SparseLDLSolver, + * EigenSimplicialLDLT) for medium-to-large systems (n > ~2000) thanks to + * supernodal BLAS3-based factorization. CHOLMOD automatically selects the best + * fill-reducing ordering (AMD, METIS, or NESDIS) regardless of the linked + * OrderingMethod component. + */ +template +class EigenCholmodSupernodalLLT + : public EigenDirectSparseSolver< + TBlockType, + MainCholmodSupernodalLLTFactory + > +{ +public: + typedef sofa::linearalgebra::CompressedRowSparseMatrix Matrix; + using Real = typename Matrix::Real; + typedef sofa::linearalgebra::FullVector Vector; + + SOFA_CLASS(SOFA_TEMPLATE(EigenCholmodSupernodalLLT, TBlockType), SOFA_TEMPLATE2(EigenDirectSparseSolver, TBlockType, MainCholmodSupernodalLLTFactory)); +}; + +#ifndef SOFA_COMPONENT_LINEARSOLVER_DIRECT_EIGENCHOLMODSUPERNODALLLT_CPP + extern template class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API EigenCholmodSupernodalLLT< SReal >; + extern template class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API EigenCholmodSupernodalLLT< sofa::type::Mat<3,3,SReal> >; +#endif + +} // namespace sofa::component::linearsolver::direct + +#endif // SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSolverFactory.h b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSolverFactory.h index ee994822382..b299dc93351 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSolverFactory.h +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSolverFactory.h @@ -28,6 +28,7 @@ #include #include #include + #include #include #include @@ -331,5 +332,27 @@ class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API MainLUFactory : public BaseMainEige } }; +#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD +/** + * Singleton factory dedicated to CHOLMOD supernodal LLT solvers. + * + * CHOLMOD manages its own fill-reducing ordering internally (AMD, METIS, NESDIS), + * so the OrderingMethodType template parameter is intentionally ignored. + * The ordering method name is only used as a key for the factory lookup. + * + * Unlike other factories, registerSolver is not defined in this header because + * it requires which we don't want to pull into every + * translation unit. It is defined in init.cpp instead. + */ +class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API MainCholmodSupernodalLLTFactory : public BaseMainEigenSolverFactory +{ +public: + ~MainCholmodSupernodalLLTFactory(); + + template + static void registerSolver(const std::string& orderingMethodName); +}; +#endif + } diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/init.cpp b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/init.cpp index 2515edac27b..c336258dffb 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/init.cpp +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/init.cpp @@ -20,6 +20,9 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #include +#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD +#include +#endif #include #include #include @@ -43,6 +46,9 @@ extern void registerEigenSparseQR(sofa::core::ObjectFactory* factory); extern void registerPrecomputedLinearSolver(sofa::core::ObjectFactory* factory); extern void registerSparseLDLSolver(sofa::core::ObjectFactory* factory); extern void registerSVDLinearSolver(sofa::core::ObjectFactory* factory); +#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD +extern void registerEigenCholmodSupernodalLLT(sofa::core::ObjectFactory* factory); +#endif extern "C" { SOFA_EXPORT_DYNAMIC_LIBRARY void initExternalModule(); @@ -83,6 +89,25 @@ void registerOrderingMethods() registerOrderingMethods(); } +#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD +// Definition of MainCholmodSupernodalLLTFactory::registerSolver. +// Kept here (not in the header) to avoid pulling everywhere. +template +void MainCholmodSupernodalLLTFactory::registerSolver(const std::string& orderingMethodName) +{ + std::lock_guard lock(s_mutex); + // CHOLMOD uses its own internal ordering; the OrderingMethodType is unused + getFactory().registerType< + Eigen::CholmodSupernodalLLT> >(orderingMethodName); +} + +void registerCholmodOrderingMethods() +{ + // CHOLMOD only supports double precision + registerOrderingMethods(); +} +#endif + void registerObjects(sofa::core::ObjectFactory* factory) { registerAsyncSparseLDLSolver(factory); @@ -97,6 +122,9 @@ void registerObjects(sofa::core::ObjectFactory* factory) registerSparseLDLSolver(factory); registerSVDLinearSolver(factory); linearsystem::registerTypedMatrixLinearSystemBTDMatrix(factory); +#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD + registerEigenCholmodSupernodalLLT(factory); +#endif } void init() @@ -109,6 +137,9 @@ void init() registerOrderingMethods(); registerOrderingMethods(); +#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD + registerCholmodOrderingMethods(); +#endif first = false; } diff --git a/cmake/Modules/FindCHOLMOD.cmake b/cmake/Modules/FindCHOLMOD.cmake new file mode 100644 index 00000000000..d394257c328 --- /dev/null +++ b/cmake/Modules/FindCHOLMOD.cmake @@ -0,0 +1,75 @@ +# Find the CHOLMOD library from SuiteSparse +# +# Behavior: first tries to find SuiteSparse CMake config files (SuiteSparse >= 7), +# then falls back to manual header/library search. +# +# Defines: +# CHOLMOD_FOUND : True if CHOLMOD is found +# CHOLMOD_INCLUDE_DIRS : CHOLMOD include directories +# CHOLMOD_LIBRARIES : CHOLMOD libraries to link against +# +# Provides target SuiteSparse::CHOLMOD when using SuiteSparse config files, +# or a manual CHOLMOD::CHOLMOD imported target otherwise. + +# Try SuiteSparse >= 7 CMake config first +find_package(SuiteSparse CONFIG QUIET COMPONENTS CHOLMOD) +if(TARGET SuiteSparse::CHOLMOD) + set(CHOLMOD_FOUND TRUE) + get_target_property(CHOLMOD_INCLUDE_DIRS SuiteSparse::CHOLMOD INTERFACE_INCLUDE_DIRECTORIES) + set(CHOLMOD_LIBRARIES SuiteSparse::CHOLMOD) + return() +endif() + +# Try CHOLMOD standalone config +find_package(CHOLMOD CONFIG QUIET) +if(TARGET CHOLMOD::CHOLMOD) + set(CHOLMOD_FOUND TRUE) + get_target_property(CHOLMOD_INCLUDE_DIRS CHOLMOD::CHOLMOD INTERFACE_INCLUDE_DIRECTORIES) + set(CHOLMOD_LIBRARIES CHOLMOD::CHOLMOD) + return() +endif() + +# Fallback: manual search +find_path(CHOLMOD_INCLUDE_DIR + NAMES cholmod.h + PATH_SUFFIXES suitesparse SuiteSparse +) + +find_library(CHOLMOD_LIBRARY + NAMES cholmod +) + +# CHOLMOD also needs SuiteSparse_config and AMD at minimum +find_library(SUITESPARSE_CONFIG_LIBRARY + NAMES suitesparseconfig +) + +find_library(AMD_LIBRARY + NAMES amd +) + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(CHOLMOD + REQUIRED_VARS CHOLMOD_LIBRARY CHOLMOD_INCLUDE_DIR +) + +if(CHOLMOD_FOUND) + set(CHOLMOD_INCLUDE_DIRS ${CHOLMOD_INCLUDE_DIR}) + set(CHOLMOD_LIBRARIES ${CHOLMOD_LIBRARY}) + if(SUITESPARSE_CONFIG_LIBRARY) + list(APPEND CHOLMOD_LIBRARIES ${SUITESPARSE_CONFIG_LIBRARY}) + endif() + if(AMD_LIBRARY) + list(APPEND CHOLMOD_LIBRARIES ${AMD_LIBRARY}) + endif() + + if(NOT TARGET CHOLMOD::CHOLMOD) + add_library(CHOLMOD::CHOLMOD UNKNOWN IMPORTED) + set_target_properties(CHOLMOD::CHOLMOD PROPERTIES + IMPORTED_LOCATION "${CHOLMOD_LIBRARY}" + INTERFACE_INCLUDE_DIRECTORIES "${CHOLMOD_INCLUDE_DIR}" + ) + endif() +endif() + +mark_as_advanced(CHOLMOD_INCLUDE_DIR CHOLMOD_LIBRARY SUITESPARSE_CONFIG_LIBRARY AMD_LIBRARY) From c2cb687042b96bf2d3625f2ac905eea2badbe347 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Wed, 8 Apr 2026 08:08:42 +0900 Subject: [PATCH 02/10] add example (and increase numbers of nodes) --- .../Direct/FEMBAR_EigenCholmodSupernodalLLT.scn | 9 +++++++++ examples/Component/LinearSolver/FEMBAR-common.xml | 2 +- 2 files changed, 10 insertions(+), 1 deletion(-) create mode 100644 examples/Component/LinearSolver/Direct/FEMBAR_EigenCholmodSupernodalLLT.scn diff --git a/examples/Component/LinearSolver/Direct/FEMBAR_EigenCholmodSupernodalLLT.scn b/examples/Component/LinearSolver/Direct/FEMBAR_EigenCholmodSupernodalLLT.scn new file mode 100644 index 00000000000..42c96b570de --- /dev/null +++ b/examples/Component/LinearSolver/Direct/FEMBAR_EigenCholmodSupernodalLLT.scn @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/examples/Component/LinearSolver/FEMBAR-common.xml b/examples/Component/LinearSolver/FEMBAR-common.xml index 36a7ea2f858..c6c9948dc87 100644 --- a/examples/Component/LinearSolver/FEMBAR-common.xml +++ b/examples/Component/LinearSolver/FEMBAR-common.xml @@ -24,7 +24,7 @@ - + From 32e042d1e979fc83b3fe3e483f6160d8c2757bb5 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 10 Jul 2026 12:19:02 +0900 Subject: [PATCH 03/10] Set EigenCholmodSupernodalLLT in its own plugin, to isolate dependency to SuiteSparse/CHOLMOD --- .../LinearSolver/Direct/CMakeLists.txt | 14 --- .../linearsolver/direct/EigenSolverFactory.h | 23 ---- .../component/linearsolver/direct/init.cpp | 31 ------ applications/plugins/CMakeLists.txt | 1 + .../plugins/SofaCHOLMOD/CMakeLists.txt | 39 +++++++ applications/plugins/SofaCHOLMOD/README.md | 0 .../SofaCHOLMOD/SofaCHOLMODConfig.cmake.in | 12 +++ .../SofaCHOLMOD/examples/FEMBAR-common.xml | 30 ++++++ .../FEMBAR_EigenCholmodSupernodalLLT.scn | 11 ++ .../EigenCholmodSupernodalLLT.cpp | 19 ++-- .../SofaCHOLMOD}/EigenCholmodSupernodalLLT.h | 41 ++++--- .../SofaCHOLMOD/src/SofaCHOLMOD/config.h.in | 37 +++++++ .../SofaCHOLMOD/src/SofaCHOLMOD/init.cpp | 102 ++++++++++++++++++ .../SofaCHOLMOD/src/SofaCHOLMOD/init.h | 29 +++++ 14 files changed, 298 insertions(+), 91 deletions(-) create mode 100644 applications/plugins/SofaCHOLMOD/CMakeLists.txt create mode 100644 applications/plugins/SofaCHOLMOD/README.md create mode 100644 applications/plugins/SofaCHOLMOD/SofaCHOLMODConfig.cmake.in create mode 100644 applications/plugins/SofaCHOLMOD/examples/FEMBAR-common.xml create mode 100644 applications/plugins/SofaCHOLMOD/examples/FEMBAR_EigenCholmodSupernodalLLT.scn rename {Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct => applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD}/EigenCholmodSupernodalLLT.cpp (74%) rename {Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct => applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD}/EigenCholmodSupernodalLLT.h (66%) create mode 100644 applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/config.h.in create mode 100644 applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/init.cpp create mode 100644 applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/init.h diff --git a/Sofa/Component/LinearSolver/Direct/CMakeLists.txt b/Sofa/Component/LinearSolver/Direct/CMakeLists.txt index 301cabb3c1d..5ee3287128b 100644 --- a/Sofa/Component/LinearSolver/Direct/CMakeLists.txt +++ b/Sofa/Component/LinearSolver/Direct/CMakeLists.txt @@ -7,7 +7,6 @@ set(HEADER_FILES ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/config.h.in ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/init.h ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/AsyncSparseLDLSolver.h - ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenCholmodSupernodalLLT.h ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/AsyncSparseLDLSolver.inl ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/BTDLinearSolver.h ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/BTDLinearSolver.inl @@ -39,7 +38,6 @@ set(HEADER_FILES set(SOURCE_FILES ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/init.cpp ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/AsyncSparseLDLSolver.cpp - ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenCholmodSupernodalLLT.cpp ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/BTDLinearSolver.cpp ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/CholeskySolver.cpp ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenSimplicialLDLT.cpp @@ -61,23 +59,11 @@ sofa_find_package(Sofa.Simulation.Core REQUIRED) sofa_find_package(Sofa.Component.LinearSolver.Iterative REQUIRED) # Only for MatrixLinearSolver which will move to core sofa_find_package(Sofa.Component.LinearSolver.Ordering REQUIRED) -# Optional CHOLMOD support from SuiteSparse for supernodal Cholesky factorization -sofa_find_package(CHOLMOD QUIET) - add_library(${PROJECT_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES} ${WRAPPER_FILES}) target_link_libraries(${PROJECT_NAME} PUBLIC Sofa.Simulation.Core Sofa.Component.LinearSolver.Iterative Sofa.Component.LinearSolver.Ordering) target_link_libraries(${PROJECT_NAME} PUBLIC Threads::Threads) -if(CHOLMOD_FOUND) - target_compile_definitions(${PROJECT_NAME} PUBLIC SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD) - target_link_libraries(${PROJECT_NAME} PUBLIC ${CHOLMOD_LIBRARIES}) - target_include_directories(${PROJECT_NAME} SYSTEM PUBLIC ${CHOLMOD_INCLUDE_DIRS}) - message(STATUS "CHOLMOD found: EigenCholmodSupernodalLLT solver enabled") -else() - message(STATUS "CHOLMOD not found: EigenCholmodSupernodalLLT solver disabled") -endif() - sofa_create_package_with_targets( PACKAGE_NAME ${PROJECT_NAME} PACKAGE_VERSION ${Sofa_VERSION} diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSolverFactory.h b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSolverFactory.h index b299dc93351..c754c35257d 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSolverFactory.h +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSolverFactory.h @@ -332,27 +332,4 @@ class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API MainLUFactory : public BaseMainEige } }; -#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD -/** - * Singleton factory dedicated to CHOLMOD supernodal LLT solvers. - * - * CHOLMOD manages its own fill-reducing ordering internally (AMD, METIS, NESDIS), - * so the OrderingMethodType template parameter is intentionally ignored. - * The ordering method name is only used as a key for the factory lookup. - * - * Unlike other factories, registerSolver is not defined in this header because - * it requires which we don't want to pull into every - * translation unit. It is defined in init.cpp instead. - */ -class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API MainCholmodSupernodalLLTFactory : public BaseMainEigenSolverFactory -{ -public: - ~MainCholmodSupernodalLLTFactory(); - - template - static void registerSolver(const std::string& orderingMethodName); -}; -#endif - - } diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/init.cpp b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/init.cpp index c336258dffb..2515edac27b 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/init.cpp +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/init.cpp @@ -20,9 +20,6 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #include -#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD -#include -#endif #include #include #include @@ -46,9 +43,6 @@ extern void registerEigenSparseQR(sofa::core::ObjectFactory* factory); extern void registerPrecomputedLinearSolver(sofa::core::ObjectFactory* factory); extern void registerSparseLDLSolver(sofa::core::ObjectFactory* factory); extern void registerSVDLinearSolver(sofa::core::ObjectFactory* factory); -#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD -extern void registerEigenCholmodSupernodalLLT(sofa::core::ObjectFactory* factory); -#endif extern "C" { SOFA_EXPORT_DYNAMIC_LIBRARY void initExternalModule(); @@ -89,25 +83,6 @@ void registerOrderingMethods() registerOrderingMethods(); } -#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD -// Definition of MainCholmodSupernodalLLTFactory::registerSolver. -// Kept here (not in the header) to avoid pulling everywhere. -template -void MainCholmodSupernodalLLTFactory::registerSolver(const std::string& orderingMethodName) -{ - std::lock_guard lock(s_mutex); - // CHOLMOD uses its own internal ordering; the OrderingMethodType is unused - getFactory().registerType< - Eigen::CholmodSupernodalLLT> >(orderingMethodName); -} - -void registerCholmodOrderingMethods() -{ - // CHOLMOD only supports double precision - registerOrderingMethods(); -} -#endif - void registerObjects(sofa::core::ObjectFactory* factory) { registerAsyncSparseLDLSolver(factory); @@ -122,9 +97,6 @@ void registerObjects(sofa::core::ObjectFactory* factory) registerSparseLDLSolver(factory); registerSVDLinearSolver(factory); linearsystem::registerTypedMatrixLinearSystemBTDMatrix(factory); -#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD - registerEigenCholmodSupernodalLLT(factory); -#endif } void init() @@ -137,9 +109,6 @@ void init() registerOrderingMethods(); registerOrderingMethods(); -#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD - registerCholmodOrderingMethods(); -#endif first = false; } diff --git a/applications/plugins/CMakeLists.txt b/applications/plugins/CMakeLists.txt index 736fbf72423..ba69551ca83 100644 --- a/applications/plugins/CMakeLists.txt +++ b/applications/plugins/CMakeLists.txt @@ -40,6 +40,7 @@ sofa_add_subdirectory(plugin SofaHAPI SofaHAPI) sofa_add_subdirectory(plugin SofaCarving SofaCarving) sofa_add_subdirectory(plugin LeapMotion LeapMotion) sofa_add_subdirectory(plugin Geomagic Geomagic) +sofa_add_subdirectory(plugin SofaCHOLMOD SofaCHOLMOD) sofa_add_external(plugin SofaAssimp GIT_REF master GIT_REPOSITORY https://www.github.com/sofa-framework/SofaAssimp.git) # ColladaSceneLoader Depends on Flexible and image sofa_add_subdirectory(plugin SofaMatrix SofaMatrix) # Depends on image, CImgPlugin sofa_add_external(plugin BeamAdapter GIT_REF master GIT_REPOSITORY https://www.github.com/sofa-framework/BeamAdapter.git) diff --git a/applications/plugins/SofaCHOLMOD/CMakeLists.txt b/applications/plugins/SofaCHOLMOD/CMakeLists.txt new file mode 100644 index 00000000000..3621355b952 --- /dev/null +++ b/applications/plugins/SofaCHOLMOD/CMakeLists.txt @@ -0,0 +1,39 @@ +cmake_minimum_required(VERSION 3.22) +project(SofaCHOLMOD LANGUAGES CXX) + +find_package(Sofa.Config REQUIRED) +sofa_find_package(Sofa.Simulation.Core REQUIRED) +sofa_find_package(Sofa.Component.LinearSolver.Direct REQUIRED) +sofa_find_package(Eigen3 REQUIRED) +sofa_find_package(CHOLMOD REQUIRED) + +get_target_property(SOFACHOLMOD_VERSION Sofa.Config SOFA_VERSION) + +set(SOFACHOLMOD_SOURCE_DIR "src/${PROJECT_NAME}") + +set(HEADER_FILES + ${SOFACHOLMOD_SOURCE_DIR}/config.h.in + ${SOFACHOLMOD_SOURCE_DIR}/init.h + ${SOFACHOLMOD_SOURCE_DIR}/EigenCholmodSupernodalLLT.h +) + +set(SOURCE_FILES + ${SOFACHOLMOD_SOURCE_DIR}/init.cpp + ${SOFACHOLMOD_SOURCE_DIR}/EigenCholmodSupernodalLLT.cpp +) + +set(TARGET_NAME ${PROJECT_NAME}) + +add_library(${TARGET_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES}) +target_link_libraries(${TARGET_NAME} PUBLIC Sofa.Config Sofa.Simulation.Core) +target_link_libraries(${TARGET_NAME} PUBLIC Sofa.Component.LinearSolver.Direct) +target_link_libraries(${PROJECT_NAME} PUBLIC ${CHOLMOD_LIBRARIES}) +target_include_directories(${PROJECT_NAME} SYSTEM PUBLIC ${CHOLMOD_INCLUDE_DIRS}) + +sofa_create_package_with_targets( + PACKAGE_NAME ${PROJECT_NAME} + PACKAGE_VERSION ${SOFACHOLMOD_VERSION} + TARGETS ${PROJECT_NAME} AUTO_SET_TARGET_PROPERTIES + INCLUDE_SOURCE_DIR "src" + INCLUDE_INSTALL_DIR "${PROJECT_NAME}" +) diff --git a/applications/plugins/SofaCHOLMOD/README.md b/applications/plugins/SofaCHOLMOD/README.md new file mode 100644 index 00000000000..e69de29bb2d diff --git a/applications/plugins/SofaCHOLMOD/SofaCHOLMODConfig.cmake.in b/applications/plugins/SofaCHOLMOD/SofaCHOLMODConfig.cmake.in new file mode 100644 index 00000000000..03ebbc1df89 --- /dev/null +++ b/applications/plugins/SofaCHOLMOD/SofaCHOLMODConfig.cmake.in @@ -0,0 +1,12 @@ +# CMake package configuration file for the SofaCHOLMOD plugin + +@PACKAGE_GUARD@ +@PACKAGE_INIT@ + +find_package(Sofa.Config QUIET REQUIRED) + +sofa_find_package(Sofa.Simulation.Core QUIET REQUIRED) +sofa_find_package(Sofa.Compoment.LinearSolver.Direct QUIET REQUIRED) +sofa_find_package(CHOLMOD QUIET REQUIRED) + +check_required_components(SofaCHOLMOD) diff --git a/applications/plugins/SofaCHOLMOD/examples/FEMBAR-common.xml b/applications/plugins/SofaCHOLMOD/examples/FEMBAR-common.xml new file mode 100644 index 00000000000..c6c9948dc87 --- /dev/null +++ b/applications/plugins/SofaCHOLMOD/examples/FEMBAR-common.xml @@ -0,0 +1,30 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/applications/plugins/SofaCHOLMOD/examples/FEMBAR_EigenCholmodSupernodalLLT.scn b/applications/plugins/SofaCHOLMOD/examples/FEMBAR_EigenCholmodSupernodalLLT.scn new file mode 100644 index 00000000000..a3450bbe72d --- /dev/null +++ b/applications/plugins/SofaCHOLMOD/examples/FEMBAR_EigenCholmodSupernodalLLT.scn @@ -0,0 +1,11 @@ + + + + + + + + + + + diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenCholmodSupernodalLLT.cpp b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp similarity index 74% rename from Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenCholmodSupernodalLLT.cpp rename to applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp index e0bc2207cbc..a6fd1a6db24 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenCholmodSupernodalLLT.cpp +++ b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp @@ -19,35 +19,32 @@ * * * Contact information: contact@sofa-framework.org * ******************************************************************************/ -#include +#define SOFACHOLMOD_EIGENCHOLMODSUPERNODALLLT_CPP -#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD +#include -#define SOFA_COMPONENT_LINEARSOLVER_DIRECT_EIGENCHOLMODSUPERNODALLLT_CPP #include -#include +#include #include #include -namespace sofa::component::linearsolver::direct +namespace sofacholmod { // CHOLMOD only supports double precision static_assert(std::is_same_v, "EigenCholmodSupernodalLLT requires double precision (SReal must be double)"); -template class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API EigenCholmodSupernodalLLT< SReal >; -template class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API EigenCholmodSupernodalLLT< sofa::type::Mat<3,3,SReal> >; +template class SOFACHOLMOD_API EigenCholmodSupernodalLLT< SReal >; +template class SOFACHOLMOD_API EigenCholmodSupernodalLLT< sofa::type::Mat<3,3,SReal> >; MainCholmodSupernodalLLTFactory::~MainCholmodSupernodalLLTFactory() = default; void registerEigenCholmodSupernodalLLT(sofa::core::ObjectFactory* factory) { - factory->registerObjects(core::ObjectRegistrationData("Direct linear solver using a supernodal sparse LL^T Cholesky factorization from CHOLMOD (SuiteSparse).") + factory->registerObjects(sofa::core::ObjectRegistrationData("Direct linear solver using a supernodal sparse LL^T Cholesky factorization from CHOLMOD (SuiteSparse).") .add< EigenCholmodSupernodalLLT< SReal > >() .add< EigenCholmodSupernodalLLT< sofa::type::Mat<3, 3, SReal> > >()); } -} // namespace sofa::component::linearsolver::direct - -#endif // SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD +} // namespace sofacholmod diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenCholmodSupernodalLLT.h b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.h similarity index 66% rename from Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenCholmodSupernodalLLT.h rename to applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.h index 46cfd950b6c..8cb1d557230 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenCholmodSupernodalLLT.h +++ b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.h @@ -20,16 +20,35 @@ * Contact information: contact@sofa-framework.org * ******************************************************************************/ #pragma once -#include - -#ifdef SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD +#include #include #include -namespace sofa::component::linearsolver::direct +namespace sofacholmod { +/** + * Singleton factory dedicated to CHOLMOD supernodal LLT solvers. + * + * CHOLMOD manages its own fill-reducing ordering internally (AMD, METIS, NESDIS), + * so the OrderingMethodType template parameter is intentionally ignored. + * The ordering method name is only used as a key for the factory lookup. + * + * Unlike other factories, registerSolver is not defined in this header because + * it requires which we don't want to pull into every + * translation unit. It is defined in init.cpp instead. + */ +class SOFACHOLMOD_API MainCholmodSupernodalLLTFactory : public sofa::component::linearsolver::direct::BaseMainEigenSolverFactory +{ + public: + ~MainCholmodSupernodalLLTFactory(); + + template + static void registerSolver(const std::string& orderingMethodName); +}; + + /** * Direct linear solver based on a supernodal sparse LL^T Cholesky factorization * using CHOLMOD from the SuiteSparse library (via Eigen's CholmodSupport). @@ -42,7 +61,7 @@ namespace sofa::component::linearsolver::direct */ template class EigenCholmodSupernodalLLT - : public EigenDirectSparseSolver< + : public sofa::component::linearsolver::direct::EigenDirectSparseSolver< TBlockType, MainCholmodSupernodalLLTFactory > @@ -52,14 +71,12 @@ class EigenCholmodSupernodalLLT using Real = typename Matrix::Real; typedef sofa::linearalgebra::FullVector Vector; - SOFA_CLASS(SOFA_TEMPLATE(EigenCholmodSupernodalLLT, TBlockType), SOFA_TEMPLATE2(EigenDirectSparseSolver, TBlockType, MainCholmodSupernodalLLTFactory)); + SOFA_CLASS(SOFA_TEMPLATE(EigenCholmodSupernodalLLT, TBlockType), SOFA_TEMPLATE2(sofa::component::linearsolver::direct::EigenDirectSparseSolver, TBlockType, MainCholmodSupernodalLLTFactory)); }; -#ifndef SOFA_COMPONENT_LINEARSOLVER_DIRECT_EIGENCHOLMODSUPERNODALLLT_CPP - extern template class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API EigenCholmodSupernodalLLT< SReal >; - extern template class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API EigenCholmodSupernodalLLT< sofa::type::Mat<3,3,SReal> >; +#ifndef SOFACHOLMOD_EIGENCHOLMODSUPERNODALLLT_CPP +extern template class SOFACHOLMOD_API EigenCholmodSupernodalLLT< SReal >; +extern template class SOFACHOLMOD_API EigenCholmodSupernodalLLT< sofa::type::Mat<3,3,SReal> >; #endif -} // namespace sofa::component::linearsolver::direct - -#endif // SOFA_COMPONENT_LINEARSOLVER_DIRECT_HAVE_CHOLMOD +} // namespace sofacholmod diff --git a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/config.h.in b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/config.h.in new file mode 100644 index 00000000000..19ab0a677cc --- /dev/null +++ b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/config.h.in @@ -0,0 +1,37 @@ +/****************************************************************************** +* 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 + +#ifdef SOFA_BUILD_SOFACHOLMOD +# define SOFACHOLMOD_API SOFA_EXPORT_DYNAMIC_LIBRARY +#else +# define SOFACHOLMOD_API SOFA_IMPORT_DYNAMIC_LIBRARY +#endif + +namespace sofacholmod +{ + constexpr const char* MODULE_NAME = "@PROJECT_NAME@"; + constexpr const char* MODULE_VERSION = "@PROJECT_VERSION@"; +} // namespace sofacholmod diff --git a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/init.cpp b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/init.cpp new file mode 100644 index 00000000000..211c522908b --- /dev/null +++ b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/init.cpp @@ -0,0 +1,102 @@ +/****************************************************************************** +* 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 +#include + +namespace sofacholmod +{ + +// Definition of MainCholmodSupernodalLLTFactory::registerSolver. +// Kept here (not in the header) to avoid pulling everywhere. +template +void MainCholmodSupernodalLLTFactory::registerSolver(const std::string& orderingMethodName) +{ + std::lock_guard lock(s_mutex); + // CHOLMOD uses its own internal ordering; the OrderingMethodType is unused + getFactory().registerType< + Eigen::CholmodSupernodalLLT> >(orderingMethodName); +} + + +template +void registerOrderingMethods() +{ + EigenSolverFactory::template registerSolver, Scalar >("AMD"); + EigenSolverFactory::template registerSolver, Scalar >("COLAMD"); + EigenSolverFactory::template registerSolver, Scalar >("Natural"); +} + +void registerCholmodOrderingMethods() +{ + // CHOLMOD only supports double precision + registerOrderingMethods(); +} + +extern void registerEigenCholmodSupernodalLLT(sofa::core::ObjectFactory* factory); + +extern "C" { + SOFA_EXPORT_DYNAMIC_LIBRARY void initExternalModule(); + SOFA_EXPORT_DYNAMIC_LIBRARY const char* getModuleName(); + SOFA_EXPORT_DYNAMIC_LIBRARY const char* getModuleVersion(); + SOFA_EXPORT_DYNAMIC_LIBRARY void registerObjects(sofa::core::ObjectFactory* factory); +} + +void initExternalModule() +{ + init(); +} + +const char* getModuleName() +{ + return MODULE_NAME; +} + +const char* getModuleVersion() +{ + return MODULE_VERSION; +} + +void registerObjects(sofa::core::ObjectFactory* factory) +{ + registerEigenCholmodSupernodalLLT(factory); +} + +void init() +{ + static bool first = true; + if (first) + { + // make sure that this plugin is registered into the PluginManager + sofa::helper::system::PluginManager::getInstance().registerPlugin(MODULE_NAME); + + registerCholmodOrderingMethods(); + + first = false; + } +} + +} // namespace sofacholmod diff --git a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/init.h b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/init.h new file mode 100644 index 00000000000..967a023939e --- /dev/null +++ b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/init.h @@ -0,0 +1,29 @@ +/****************************************************************************** +* 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 + +namespace sofacholmod +{ + SOFACHOLMOD_API void init(); +} // namespace sofacholmod From 08a1e09efc59e37f43fca03c6f66677f9c8fa7f2 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 10 Jul 2026 16:12:29 +0900 Subject: [PATCH 04/10] add Data to control the nb of threads for BLAS (if applicable) --- .../plugins/SofaCHOLMOD/CMakeLists.txt | 4 + .../FEMBAR_EigenCholmodSupernodalLLT.scn | 2 +- .../SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp | 119 ++++++++++++++++++ .../SofaCHOLMOD/EigenCholmodSupernodalLLT.h | 22 ++++ 4 files changed, 146 insertions(+), 1 deletion(-) diff --git a/applications/plugins/SofaCHOLMOD/CMakeLists.txt b/applications/plugins/SofaCHOLMOD/CMakeLists.txt index 3621355b952..eb8bfc59361 100644 --- a/applications/plugins/SofaCHOLMOD/CMakeLists.txt +++ b/applications/plugins/SofaCHOLMOD/CMakeLists.txt @@ -30,6 +30,10 @@ target_link_libraries(${TARGET_NAME} PUBLIC Sofa.Component.LinearSolver.Direct) target_link_libraries(${PROJECT_NAME} PUBLIC ${CHOLMOD_LIBRARIES}) target_include_directories(${PROJECT_NAME} SYSTEM PUBLIC ${CHOLMOD_INCLUDE_DIRS}) +# Needed to resolve the BLAS backend's runtime thread-count setter via dlsym +# (empty on Windows, where GetProcAddress from kernel32 is used instead). +target_link_libraries(${PROJECT_NAME} PRIVATE ${CMAKE_DL_LIBS}) + sofa_create_package_with_targets( PACKAGE_NAME ${PROJECT_NAME} PACKAGE_VERSION ${SOFACHOLMOD_VERSION} diff --git a/applications/plugins/SofaCHOLMOD/examples/FEMBAR_EigenCholmodSupernodalLLT.scn b/applications/plugins/SofaCHOLMOD/examples/FEMBAR_EigenCholmodSupernodalLLT.scn index a3450bbe72d..28ca143592d 100644 --- a/applications/plugins/SofaCHOLMOD/examples/FEMBAR_EigenCholmodSupernodalLLT.scn +++ b/applications/plugins/SofaCHOLMOD/examples/FEMBAR_EigenCholmodSupernodalLLT.scn @@ -5,7 +5,7 @@ - + diff --git a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp index a6fd1a6db24..5e772ef1259 100644 --- a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp +++ b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp @@ -29,12 +29,131 @@ #include +#if defined(_WIN32) +#include +#else +#include +#endif + namespace sofacholmod { // CHOLMOD only supports double precision static_assert(std::is_same_v, "EigenCholmodSupernodalLLT requires double precision (SReal must be double)"); +namespace +{ + +/// Try to set the number of threads used by the BLAS backend at runtime. +/// +/// CHOLMOD's supernodal factorization delegates to dense BLAS3 kernels, so the +/// thread count that matters is the BLAS one. It cannot be changed through +/// environment variables once the process is running (OpenBLAS reads them when +/// the library is loaded, before this component is created), hence we resolve +/// the backend's runtime setter dynamically and call it. +/// +/// Returns true if a supported BLAS thread-control entry point was found. +bool trySetBlasNumThreads(int numThreads) +{ + using SetNumThreadsFn = void (*)(int); + + // Value-based setters exported by the common optimized BLAS backends. + static const char* const candidates[] = { + "openblas_set_num_threads", // OpenBLAS + "goto_set_num_threads", // GotoBLAS / older OpenBLAS + "MKL_Set_Num_Threads", // Intel MKL + }; + +#if defined(_WIN32) + // The BLAS DLL is already loaded (pulled in by CHOLMOD); look it up by name. + static const char* const moduleNames[] = { + "libopenblas", "openblas", "libblas", "mkl_rt" + }; + for (const char* moduleName : moduleNames) + { + HMODULE module = GetModuleHandleA(moduleName); + if (!module) + { + continue; + } + for (const char* name : candidates) + { + if (auto fn = reinterpret_cast( + reinterpret_cast(GetProcAddress(module, name)))) + { + fn(numThreads); + return true; + } + } + } + return false; +#else + // The BLAS symbols are already loaded in the process (via CHOLMOD), + // so search the global symbol table. + for (const char* name : candidates) + { + if (auto fn = reinterpret_cast(dlsym(RTLD_DEFAULT, name))) + { + fn(numThreads); + return true; + } + } + return false; +#endif +} + +} // namespace + +template +EigenCholmodSupernodalLLT::EigenCholmodSupernodalLLT() + : d_numThreads(initData(&d_numThreads, 0, + "numThreads", + "Number of threads used by the BLAS backend of CHOLMOD's supernodal " + "factorization. A value <= 0 (default) keeps the BLAS default " + "(OPENBLAS_NUM_THREADS / OMP_NUM_THREADS). A small value (1-4) is often " + "faster on medium-sized systems by avoiding thread oversubscription. " + "Only effective with OpenBLAS or MKL.")) +{ +} + +template +void EigenCholmodSupernodalLLT::init() +{ + Inherit1::init(); + applyBlasNumThreads(); +} + +template +void EigenCholmodSupernodalLLT::reinit() +{ + Inherit1::reinit(); + applyBlasNumThreads(); +} + +template +void EigenCholmodSupernodalLLT::applyBlasNumThreads() +{ + const int numThreads = d_numThreads.getValue(); + if (numThreads <= 0) + { + // Keep the BLAS default (controlled by the environment). + return; + } + + if (trySetBlasNumThreads(numThreads)) + { + msg_info() << "BLAS backend set to use " << numThreads + << " thread(s) for CHOLMOD factorization."; + } + else + { + msg_warning() << "Could not set the number of BLAS threads at runtime: no supported " + "BLAS thread-control API (OpenBLAS/MKL) was found. The value of '" + << d_numThreads.getName() << "' is ignored on this platform. Control the " + "thread count via environment variables instead (e.g. OPENBLAS_NUM_THREADS)."; + } +} + template class SOFACHOLMOD_API EigenCholmodSupernodalLLT< SReal >; template class SOFACHOLMOD_API EigenCholmodSupernodalLLT< sofa::type::Mat<3,3,SReal> >; diff --git a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.h b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.h index 8cb1d557230..a418e9e4cfd 100644 --- a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.h +++ b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.h @@ -24,6 +24,7 @@ #include #include +#include namespace sofacholmod { @@ -72,6 +73,27 @@ class EigenCholmodSupernodalLLT typedef sofa::linearalgebra::FullVector Vector; SOFA_CLASS(SOFA_TEMPLATE(EigenCholmodSupernodalLLT, TBlockType), SOFA_TEMPLATE2(sofa::component::linearsolver::direct::EigenDirectSparseSolver, TBlockType, MainCholmodSupernodalLLTFactory)); + + void init() override; + void reinit() override; + + /// Number of threads the underlying BLAS (used by CHOLMOD's supernodal + /// factorization) is allowed to use. + /// A value <= 0 (default) leaves the BLAS default untouched, i.e. controlled + /// by the OPENBLAS_NUM_THREADS / OMP_NUM_THREADS environment variables. + /// For the medium-sized systems typically solved here, a small value (1-4) + /// is often faster than the default (= number of cores) because it avoids + /// thread oversubscription. Only effective with OpenBLAS or MKL; ignored + /// with BLAS backends that expose no runtime thread-control API (e.g. Apple + /// Accelerate), where the environment variables must be used instead. + sofa::core::objectmodel::Data d_numThreads; + +protected: + EigenCholmodSupernodalLLT(); + + /// Apply d_numThreads to the underlying BLAS backend if it exposes a runtime + /// thread-count API (OpenBLAS/MKL). No-op if d_numThreads <= 0. + void applyBlasNumThreads(); }; #ifndef SOFACHOLMOD_EIGENCHOLMODSUPERNODALLLT_CPP From 4e8deef0e83e16ac2af42ab8e3e486d57814e95c Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 10 Jul 2026 17:17:48 +0900 Subject: [PATCH 05/10] add README --- applications/plugins/SofaCHOLMOD/README.md | 132 +++++++++++++++++++++ 1 file changed, 132 insertions(+) diff --git a/applications/plugins/SofaCHOLMOD/README.md b/applications/plugins/SofaCHOLMOD/README.md index e69de29bb2d..2688417bf62 100644 --- a/applications/plugins/SofaCHOLMOD/README.md +++ b/applications/plugins/SofaCHOLMOD/README.md @@ -0,0 +1,132 @@ +# SofaCHOLMOD + +A SOFA plugin providing a direct linear solver based on the **supernodal sparse +Cholesky factorization** of [CHOLMOD](https://github.com/DrTimothyAldenDavis/SuiteSparse) +(from SuiteSparse), wrapped through [Eigen's CholmodSupport](https://eigen.tuxfamily.org/dox/group__CholmodSupport__Module.html). + +The plugin is kept separate from the core SOFA solvers so that the dependency on +SuiteSparse / CHOLMOD stays isolated to this optional plugin. + +## Component + +### `EigenCholmodSupernodalLLT` + +A direct solver computing a supernodal `L·Lᵀ` Cholesky factorization of a +symmetric positive-definite system. + +Unlike simplicial solvers (`SparseLDLSolver`, `EigenSimplicialLDLT`), the +supernodal factorization groups columns with a similar sparsity pattern into +dense blocks and factorizes them with **dense BLAS3 kernels** (`dgemm`, +`dsyrk`, `dpotrf`). For medium-to-large systems this is dramatically faster — +we measured up to **~8-11x** faster than `SparseLDLSolver` on a 15k-DOF FEM beam +(see [Performance](#performance)). + +CHOLMOD manages its own fill-reducing ordering internally (AMD / METIS / +NESDIS). An `OrderingMethod` component is still required in the scene (it is part +of the common linear-solver API), but its choice is **ignored** by this solver. + +> **Note:** CHOLMOD only supports **double** precision. The plugin will not build +> if SOFA is configured with `SReal = float`. + +#### Templates + +| `template` value | Block type | +|--------------------------------------|-----------------------| +| `CompressedRowSparseMatrix` | scalar (`SReal`) | +| `CompressedRowSparseMatrixMat3x3` | `Mat<3,3,SReal>` | + +`CompressedRowSparseMatrixMat3x3` is recommended for 3D mechanics (3-DOF nodes). + +#### Data + +| Data | Type | Default | Description | +|--------------|-------|---------|-------------| +| `numThreads` | `int` | `0` | Number of threads the underlying BLAS backend may use for the factorization. `<= 0` leaves the BLAS default untouched (controlled by the `OPENBLAS_NUM_THREADS` / `OMP_NUM_THREADS` environment variables). A small value (`1`-`4`) is often faster on medium-sized systems by avoiding thread oversubscription. Only effective with OpenBLAS or MKL (see [Performance](#performance)). | + +## Usage + +```xml + + + + + + +``` + +A complete example is provided in +[`examples/FEMBAR_EigenCholmodSupernodalLLT.scn`](examples/FEMBAR_EigenCholmodSupernodalLLT.scn). + +## Dependencies + +- **SuiteSparse / CHOLMOD** (with its `SuiteSparse_config` and `AMD` modules). +- An **optimized BLAS/LAPACK** implementation backing CHOLMOD (OpenBLAS, Intel + MKL, Apple Accelerate, …). This is the single most important factor for + performance — see below. + +Installing the dependency: + +| Platform | Command | +|----------|---------| +| Ubuntu / Debian | `sudo apt install libsuitesparse-dev libopenblas-dev` | +| macOS (Homebrew) | `brew install suite-sparse` (links Apple Accelerate) | +| vcpkg | `vcpkg install suitesparse` | + +The plugin locates CHOLMOD via [`cmake/Modules/FindCHOLMOD.cmake`](../../../cmake/Modules/FindCHOLMOD.cmake), +which first tries the SuiteSparse CMake config (SuiteSparse >= 7) and falls back +to a manual header/library search. + +## Building + +Enable the plugin when configuring SOFA: + +```bash +cmake -DPLUGIN_SOFACHOLMOD=ON +``` + +## Performance + +CHOLMOD's supernodal factorization delegates almost all of its work to dense +BLAS3 kernels, so **its speed is dictated by whichever BLAS is linked behind +CHOLMOD**, not by SOFA. Keep this in mind when benchmarking: + +- On **macOS**, SuiteSparse links Apple **Accelerate** (fast), so the speedup is + always present. +- On **Linux**, the BLAS backing `libcholmod` is selected by + `update-alternatives`. If it resolves to the reference *netlib* BLAS instead of + an optimized one (OpenBLAS/MKL), the speedup largely disappears. Check and fix + with: + + ```bash + # what BLAS does CHOLMOD use? + ldd $(ldconfig -p | awk '/libcholmod.so/{print $NF; exit}') | grep -i blas + + # prefer OpenBLAS system-wide + sudo update-alternatives --config libblas.so.3-x86_64-linux-gnu + sudo update-alternatives --config liblapack.so.3-x86_64-linux-gnu + ``` + +### Threading + +With an optimized multithreaded BLAS (e.g. OpenBLAS-pthread), the default +behavior is to use *all* cores. For the medium-sized systems typical in SOFA, +this **oversubscribes** threads and the launch/synchronization overhead can make +the solver *slower* than with a handful of threads. Tune it with the +`numThreads` Data (or the `OPENBLAS_NUM_THREADS` environment variable). + +Illustrative measurement — `examples/FEMBAR_EigenCholmodSupernodalLLT.scn` +(10×10×50 grid, ~15k DOF), 24-core Linux box, OpenBLAS backend: + +| Solver / configuration | ms / step | Speedup | +|-----------------------------------------------|-----------|---------| +| `SparseLDLSolver` | ~762 | 1.0x | +| `EigenCholmodSupernodalLLT`, reference BLAS | ~346 | 2.2x | +| `EigenCholmodSupernodalLLT`, OpenBLAS, all cores | ~92 | 8.3x | +| `EigenCholmodSupernodalLLT`, OpenBLAS, `numThreads="1"` | ~70 | 10.9x | + +> The optimal thread count is problem- and machine-dependent; small values +> (`1`-`4`) are a good starting point for medium systems. + +## License + +LGPL 2.1+ — see the SOFA license. From 5f5072da3a21bb1d208377b3c11deebcb081eff5 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 10 Jul 2026 17:25:37 +0900 Subject: [PATCH 06/10] Move check on sofa_floating_point=double in cmake (not at compilation time) --- applications/plugins/SofaCHOLMOD/CMakeLists.txt | 11 +++++++++++ .../src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp | 4 ++-- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/applications/plugins/SofaCHOLMOD/CMakeLists.txt b/applications/plugins/SofaCHOLMOD/CMakeLists.txt index eb8bfc59361..62f7466f734 100644 --- a/applications/plugins/SofaCHOLMOD/CMakeLists.txt +++ b/applications/plugins/SofaCHOLMOD/CMakeLists.txt @@ -2,6 +2,17 @@ cmake_minimum_required(VERSION 3.22) project(SofaCHOLMOD LANGUAGES CXX) find_package(Sofa.Config REQUIRED) + +# CHOLMOD only supports double precision, so this plugin requires SOFA to be +# configured with SOFA_FLOATING_POINT_TYPE=double. +if(NOT "${SOFA_FLOATING_POINT_TYPE}" STREQUAL "double") + message(FATAL_ERROR + "SofaCHOLMOD requires SOFA_FLOATING_POINT_TYPE=double (CHOLMOD only " + "supports double precision), but it is currently set to " + "'${SOFA_FLOATING_POINT_TYPE}'. Set SOFA_FLOATING_POINT_TYPE=double or " + "disable this plugin (PLUGIN_SOFACHOLMOD=OFF).") +endif() + sofa_find_package(Sofa.Simulation.Core REQUIRED) sofa_find_package(Sofa.Component.LinearSolver.Direct REQUIRED) sofa_find_package(Eigen3 REQUIRED) diff --git a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp index 5e772ef1259..487e325081d 100644 --- a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp +++ b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp @@ -38,8 +38,8 @@ namespace sofacholmod { -// CHOLMOD only supports double precision -static_assert(std::is_same_v, "EigenCholmodSupernodalLLT requires double precision (SReal must be double)"); +// CHOLMOD only supports double precision. This is enforced at configure time by +// the plugin's CMakeLists.txt, which requires SOFA_FLOATING_POINT_TYPE=double. namespace { From 11d9f02132a9cfe1a3c4f0ad51581dafac9a2f77 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 10 Jul 2026 17:35:17 +0900 Subject: [PATCH 07/10] add comparison with SparseLDL --- .../plugins/SofaCHOLMOD/examples/FEMBAR_SparseLDL.scn | 9 +++++++++ examples/Component/LinearSolver/FEMBAR-common.xml | 2 +- 2 files changed, 10 insertions(+), 1 deletion(-) create mode 100644 applications/plugins/SofaCHOLMOD/examples/FEMBAR_SparseLDL.scn diff --git a/applications/plugins/SofaCHOLMOD/examples/FEMBAR_SparseLDL.scn b/applications/plugins/SofaCHOLMOD/examples/FEMBAR_SparseLDL.scn new file mode 100644 index 00000000000..e659e6f5e02 --- /dev/null +++ b/applications/plugins/SofaCHOLMOD/examples/FEMBAR_SparseLDL.scn @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/examples/Component/LinearSolver/FEMBAR-common.xml b/examples/Component/LinearSolver/FEMBAR-common.xml index c6c9948dc87..36a7ea2f858 100644 --- a/examples/Component/LinearSolver/FEMBAR-common.xml +++ b/examples/Component/LinearSolver/FEMBAR-common.xml @@ -24,7 +24,7 @@ - + From 4ca369c7b22261a23c53ec75a3f43492460c5c13 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 10 Jul 2026 17:42:29 +0900 Subject: [PATCH 08/10] remove exameple in sofa codebase --- .../Direct/FEMBAR_EigenCholmodSupernodalLLT.scn | 9 --------- 1 file changed, 9 deletions(-) delete mode 100644 examples/Component/LinearSolver/Direct/FEMBAR_EigenCholmodSupernodalLLT.scn diff --git a/examples/Component/LinearSolver/Direct/FEMBAR_EigenCholmodSupernodalLLT.scn b/examples/Component/LinearSolver/Direct/FEMBAR_EigenCholmodSupernodalLLT.scn deleted file mode 100644 index 42c96b570de..00000000000 --- a/examples/Component/LinearSolver/Direct/FEMBAR_EigenCholmodSupernodalLLT.scn +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - From 42e7d4bed0eeb5be8ec74eba3015664f11add94d Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 10 Jul 2026 22:26:44 +0900 Subject: [PATCH 09/10] add addJMInvJtLocal using dsyrk --- .../linearsolver/direct/EigenSolverFactory.h | 10 + applications/plugins/SofaCHOLMOD/README.md | 68 ++- .../SofaCHOLMOD/examples/TorusFall.scn | 443 ++++++++++++++++++ .../SofaCHOLMOD/examples/mesh/torus1-16x8.obj | 385 +++++++++++++++ .../src/SofaCHOLMOD/CholmodSolverProxy.h | 151 ++++++ .../SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp | 214 +++++++-- .../SofaCHOLMOD/EigenCholmodSupernodalLLT.h | 29 +- .../SofaCHOLMOD/src/SofaCHOLMOD/init.cpp | 9 +- 8 files changed, 1251 insertions(+), 58 deletions(-) create mode 100644 applications/plugins/SofaCHOLMOD/examples/TorusFall.scn create mode 100644 applications/plugins/SofaCHOLMOD/examples/mesh/torus1-16x8.obj create mode 100644 applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/CholmodSolverProxy.h diff --git a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSolverFactory.h b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSolverFactory.h index c754c35257d..4867161a74a 100644 --- a/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSolverFactory.h +++ b/Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSolverFactory.h @@ -181,6 +181,16 @@ class SOFA_COMPONENT_LINEARSOLVER_DIRECT_API EigenSolverFactory m_registeredTypes[join(orderingMethodName)] = [](){ return new EigenSolverWrapper(); }; } + /// Register a solver proxy that already derives from BaseEigenSolverProxy, + /// instead of wrapping a raw Eigen solver in an EigenSolverWrapper. This lets + /// a solver expose backend-specific functionality (e.g. CHOLMOD's raw factor + /// for optimized compliance-matrix computation) through its own proxy type. + template + void registerProxyType(const std::string& orderingMethodName) + { + m_registeredTypes[join(orderingMethodName)] = [](){ return new ProxyClass(); }; + } + template BaseEigenSolverProxy* getObject(const std::string& orderingMethodName ) { diff --git a/applications/plugins/SofaCHOLMOD/README.md b/applications/plugins/SofaCHOLMOD/README.md index 2688417bf62..bf5311a2313 100644 --- a/applications/plugins/SofaCHOLMOD/README.md +++ b/applications/plugins/SofaCHOLMOD/README.md @@ -41,7 +41,7 @@ of the common linear-solver API), but its choice is **ignored** by this solver. | Data | Type | Default | Description | |--------------|-------|---------|-------------| -| `numThreads` | `int` | `0` | Number of threads the underlying BLAS backend may use for the factorization. `<= 0` leaves the BLAS default untouched (controlled by the `OPENBLAS_NUM_THREADS` / `OMP_NUM_THREADS` environment variables). A small value (`1`-`4`) is often faster on medium-sized systems by avoiding thread oversubscription. Only effective with OpenBLAS or MKL (see [Performance](#performance)). | +| `numThreads` | `int` | `1` | Number of threads the underlying BLAS backend may use for the factorization. The default of `1` is the fastest setting for the vast majority of SOFA scenes and avoids thread oversubscription (see [Threading](#threading)). Increase it only for a single, very large standalone system. A value `<= 0` leaves the BLAS default untouched (controlled by the `OPENBLAS_NUM_THREADS` / `OMP_NUM_THREADS` environment variables). Only effective with OpenBLAS or MKL. | ## Usage @@ -108,24 +108,62 @@ CHOLMOD**, not by SOFA. Keep this in mind when benchmarking: ### Threading -With an optimized multithreaded BLAS (e.g. OpenBLAS-pthread), the default -behavior is to use *all* cores. For the medium-sized systems typical in SOFA, -this **oversubscribes** threads and the launch/synchronization overhead can make -the solver *slower* than with a handful of threads. Tune it with the -`numThreads` Data (or the `OPENBLAS_NUM_THREADS` environment variable). +An optimized multithreaded BLAS (e.g. OpenBLAS-pthread) uses *all* cores by +default. For the medium-sized systems typical in SOFA this **oversubscribes** +threads: the launch/synchronization overhead makes the solver *slower* than with +a single thread, and the problem becomes catastrophic when several solvers are +factorized in parallel (e.g. multiple objects with `parallelODESolving="true"`), +where each solver would otherwise spawn one BLAS pool per core. -Illustrative measurement — `examples/FEMBAR_EigenCholmodSupernodalLLT.scn` -(10×10×50 grid, ~15k DOF), 24-core Linux box, OpenBLAS backend: +For these reasons `numThreads` **defaults to `1`**, which is the fastest setting +in nearly all cases. Increase it only for a single, very large standalone +system. -| Solver / configuration | ms / step | Speedup | -|-----------------------------------------------|-----------|---------| -| `SparseLDLSolver` | ~762 | 1.0x | -| `EigenCholmodSupernodalLLT`, reference BLAS | ~346 | 2.2x | -| `EigenCholmodSupernodalLLT`, OpenBLAS, all cores | ~92 | 8.3x | +Illustrative measurements on a 24-core Linux box, OpenBLAS backend: + +Single medium system — `examples/FEMBAR_EigenCholmodSupernodalLLT.scn` +(10×10×50 grid, ~15k DOF): + +| Solver / configuration | ms / step | Speedup | +|-------------------------------------------------|-----------|---------| +| `SparseLDLSolver` | ~762 | 1.0x | +| `EigenCholmodSupernodalLLT`, reference BLAS | ~346 | 2.2x | +| `EigenCholmodSupernodalLLT`, OpenBLAS, all cores| ~92 | 8.3x | | `EigenCholmodSupernodalLLT`, OpenBLAS, `numThreads="1"` | ~70 | 10.9x | -> The optimal thread count is problem- and machine-dependent; small values -> (`1`-`4`) are a good starting point for medium systems. +Many small systems in parallel — `examples/TorusFall.scn` (10 tori, +`parallelODESolving="true"`): + +| Solver / configuration | FPS | Speedup | +|-------------------------------------------------|-------|---------| +| `EigenCholmodSupernodalLLT`, OpenBLAS, all cores| ~20 | 0.5x | +| `SparseLDLSolver` (`parallelInverseProduct`) | ~40 | 1.0x | +| `EigenCholmodSupernodalLLT`, `numThreads="1"` | ~108 | 2.7x | + +> Multithreaded BLAS only pays off for a single, very large system; whenever +> several solvers run concurrently, keep `numThreads="1"`. + +## Constraint solving (compliance matrix) + +For Lagrangian constraints (e.g. contacts), `addJMInvJtLocal` computes the +compliance block `W = J·A⁻¹·Jᵀ`. Since CHOLMOD factorizes `P·A·Pᵀ = L·Lᵀ`, this +is `W = Zᵀ·Z` with `Z = L⁻¹·P·Jᵀ` — one triangular forward-solve (all constraint +columns at once) plus one symmetric product, instead of a full solve per row. + +Implementation notes: + +- The `Z = L⁻¹·P·Jᵀ` solve uses **`cholmod_solve2`** with workspace buffers held + on the solver's `CholmodSolverProxy` and **reused across calls**, avoiding a + per-step allocate/free of the dense right-hand side (relevant for scenes with + many small systems solved every step). +- `W = Zᵀ·Z` is computed with the BLAS symmetric rank-k update **`dsyrk`** (only + the triangle that is needed, half the flops of a full product), threaded per + `numThreads`; it falls back to an Eigen rank update if `dsyrk` is unavailable. + +Performance vs `SparseLDLSolver` is size-dependent: CHOLMOD's compliance is +faster for systems above ~1000 DOF per object (and the margin grows with size), +while for very small systems the fixed per-call overhead makes the two solvers +comparable. ## License diff --git a/applications/plugins/SofaCHOLMOD/examples/TorusFall.scn b/applications/plugins/SofaCHOLMOD/examples/TorusFall.scn new file mode 100644 index 00000000000..4b9f3d5f26f --- /dev/null +++ b/applications/plugins/SofaCHOLMOD/examples/TorusFall.scn @@ -0,0 +1,443 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/applications/plugins/SofaCHOLMOD/examples/mesh/torus1-16x8.obj b/applications/plugins/SofaCHOLMOD/examples/mesh/torus1-16x8.obj new file mode 100644 index 00000000000..dbfe8da3718 --- /dev/null +++ b/applications/plugins/SofaCHOLMOD/examples/mesh/torus1-16x8.obj @@ -0,0 +1,385 @@ +# AXIS=Y R1=1 R2=3 NL=16 NC=8 php torus.obj.php +v 0 0 3 +v 0 0.70710678118655 2.7071067811865 +v 0 1 2 +v 0 0.70710678118655 1.2928932188135 +v 0 1.2246467991474E-16 1 +v 0 -0.70710678118655 1.2928932188135 +v 0 -1 2 +v 0 -0.70710678118655 2.7071067811865 +v 1.1480502970953 0 2.7716385975339 +v 1.0359649148033 0.70710678118655 2.5010405474608 +v 0.76536686473018 1 1.8477590650226 +v 0.49476881465708 0.70710678118655 1.1944775825844 +v 0.38268343236509 1.2246467991474E-16 0.92387953251129 +v 0.49476881465708 -0.70710678118655 1.1944775825844 +v 0.76536686473018 -1 1.8477590650226 +v 1.0359649148033 -0.70710678118655 2.5010405474608 +v 2.1213203435596 0 2.1213203435596 +v 1.9142135623731 0.70710678118655 1.9142135623731 +v 1.4142135623731 1 1.4142135623731 +v 0.9142135623731 0.70710678118655 0.9142135623731 +v 0.70710678118655 1.2246467991474E-16 0.70710678118655 +v 0.91421356237309 -0.70710678118655 0.9142135623731 +v 1.4142135623731 -1 1.4142135623731 +v 1.9142135623731 -0.70710678118655 1.9142135623731 +v 2.7716385975339 0 1.1480502970953 +v 2.5010405474608 0.70710678118655 1.0359649148033 +v 1.8477590650226 1 0.76536686473018 +v 1.1944775825844 0.70710678118655 0.49476881465708 +v 0.92387953251129 1.2246467991474E-16 0.38268343236509 +v 1.1944775825844 -0.70710678118655 0.49476881465708 +v 1.8477590650226 -1 0.76536686473018 +v 2.5010405474608 -0.70710678118655 1.0359649148033 +v 3 0 1.836970198721E-16 +v 2.7071067811865 0.70710678118655 1.6576248272651E-16 +v 2 1 1.2246467991474E-16 +v 1.2928932188135 0.70710678118655 7.9166877102961E-17 +v 1 1.2246467991474E-16 6.1232339957368E-17 +v 1.2928932188135 -0.70710678118655 7.9166877102961E-17 +v 2 -1 1.2246467991474E-16 +v 2.7071067811865 -0.70710678118655 1.6576248272651E-16 +v 2.7716385975339 0 -1.1480502970953 +v 2.5010405474608 0.70710678118655 -1.0359649148033 +v 1.8477590650226 1 -0.76536686473018 +v 1.1944775825844 0.70710678118655 -0.49476881465708 +v 0.92387953251129 1.2246467991474E-16 -0.38268343236509 +v 1.1944775825844 -0.70710678118655 -0.49476881465708 +v 1.8477590650226 -1 -0.76536686473018 +v 2.5010405474608 -0.70710678118655 -1.0359649148033 +v 2.1213203435596 0 -2.1213203435596 +v 1.9142135623731 0.70710678118655 -1.9142135623731 +v 1.4142135623731 1 -1.4142135623731 +v 0.9142135623731 0.70710678118655 -0.9142135623731 +v 0.70710678118655 1.2246467991474E-16 -0.70710678118655 +v 0.9142135623731 -0.70710678118655 -0.91421356237309 +v 1.4142135623731 -1 -1.4142135623731 +v 1.9142135623731 -0.70710678118655 -1.9142135623731 +v 1.1480502970953 0 -2.7716385975339 +v 1.0359649148033 0.70710678118655 -2.5010405474608 +v 0.76536686473018 1 -1.8477590650226 +v 0.49476881465708 0.70710678118655 -1.1944775825844 +v 0.38268343236509 1.2246467991474E-16 -0.92387953251129 +v 0.49476881465708 -0.70710678118655 -1.1944775825844 +v 0.76536686473018 -1 -1.8477590650226 +v 1.0359649148033 -0.70710678118655 -2.5010405474608 +v 3.6739403974421E-16 0 -3 +v 3.3152496545302E-16 0.70710678118655 -2.7071067811865 +v 2.4492935982947E-16 1 -2 +v 1.5833375420592E-16 0.70710678118655 -1.2928932188135 +v 1.2246467991474E-16 1.2246467991474E-16 -1 +v 1.5833375420592E-16 -0.70710678118655 -1.2928932188135 +v 2.4492935982947E-16 -1 -2 +v 3.3152496545302E-16 -0.70710678118655 -2.7071067811865 +v -1.1480502970953 0 -2.7716385975339 +v -1.0359649148033 0.70710678118655 -2.5010405474608 +v -0.76536686473018 1 -1.8477590650226 +v -0.49476881465708 0.70710678118655 -1.1944775825844 +v -0.38268343236509 1.2246467991474E-16 -0.92387953251129 +v -0.49476881465708 -0.70710678118655 -1.1944775825844 +v -0.76536686473018 -1 -1.8477590650226 +v -1.0359649148033 -0.70710678118655 -2.5010405474608 +v -2.1213203435596 0 -2.1213203435596 +v -1.9142135623731 0.70710678118655 -1.9142135623731 +v -1.4142135623731 1 -1.4142135623731 +v -0.9142135623731 0.70710678118655 -0.9142135623731 +v -0.70710678118655 1.2246467991474E-16 -0.70710678118655 +v -0.91421356237309 -0.70710678118655 -0.9142135623731 +v -1.4142135623731 -1 -1.4142135623731 +v -1.9142135623731 -0.70710678118655 -1.9142135623731 +v -2.7716385975339 0 -1.1480502970953 +v -2.5010405474608 0.70710678118655 -1.0359649148033 +v -1.8477590650226 1 -0.76536686473018 +v -1.1944775825844 0.70710678118655 -0.49476881465708 +v -0.92387953251129 1.2246467991474E-16 -0.38268343236509 +v -1.1944775825844 -0.70710678118655 -0.49476881465708 +v -1.8477590650226 -1 -0.76536686473018 +v -2.5010405474608 -0.70710678118655 -1.0359649148033 +v -3 0 -5.5109105961631E-16 +v -2.7071067811865 0.70710678118655 -4.9728744817953E-16 +v -2 1 -3.6739403974421E-16 +v -1.2928932188135 0.70710678118655 -2.3750063130888E-16 +v -1 1.2246467991474E-16 -1.836970198721E-16 +v -1.2928932188135 -0.70710678118655 -2.3750063130888E-16 +v -2 -1 -3.6739403974421E-16 +v -2.7071067811865 -0.70710678118655 -4.9728744817953E-16 +v -2.7716385975339 0 1.1480502970953 +v -2.5010405474608 0.70710678118655 1.0359649148033 +v -1.8477590650226 1 0.76536686473018 +v -1.1944775825844 0.70710678118655 0.49476881465708 +v -0.92387953251129 1.2246467991474E-16 0.38268343236509 +v -1.1944775825844 -0.70710678118655 0.49476881465708 +v -1.8477590650226 -1 0.76536686473018 +v -2.5010405474608 -0.70710678118655 1.0359649148033 +v -2.1213203435596 0 2.1213203435596 +v -1.9142135623731 0.70710678118655 1.9142135623731 +v -1.4142135623731 1 1.4142135623731 +v -0.9142135623731 0.70710678118655 0.91421356237309 +v -0.70710678118655 1.2246467991474E-16 0.70710678118655 +v -0.9142135623731 -0.70710678118655 0.91421356237309 +v -1.4142135623731 -1 1.4142135623731 +v -1.9142135623731 -0.70710678118655 1.9142135623731 +v -1.1480502970953 0 2.7716385975339 +v -1.0359649148033 0.70710678118655 2.5010405474608 +v -0.76536686473018 1 1.8477590650226 +v -0.49476881465708 0.70710678118655 1.1944775825844 +v -0.38268343236509 1.2246467991474E-16 0.92387953251129 +v -0.49476881465708 -0.70710678118655 1.1944775825844 +v -0.76536686473018 -1 1.8477590650226 +v -1.0359649148033 -0.70710678118655 2.5010405474608 +vn 0 0 1 +vn 0 0.70710678118655 0.70710678118655 +vn 0 1 6.1232339957368E-17 +vn -0 0.70710678118655 -0.70710678118655 +vn -0 1.2246467991474E-16 -1 +vn -0 -0.70710678118655 -0.70710678118655 +vn -0 -1 -1.836970198721E-16 +vn 0 -0.70710678118655 0.70710678118655 +vn 0.38268343236509 0 0.92387953251129 +vn 0.2705980500731 0.70710678118655 0.65328148243819 +vn 2.3432602026631E-17 1 5.6571305614385E-17 +vn -0.2705980500731 0.70710678118655 -0.65328148243819 +vn -0.38268343236509 1.2246467991474E-16 -0.92387953251129 +vn -0.2705980500731 -0.70710678118655 -0.65328148243819 +vn -7.0297806079894E-17 -1 -1.6971391684316E-16 +vn 0.2705980500731 -0.70710678118655 0.65328148243819 +vn 0.70710678118655 0 0.70710678118655 +vn 0.5 0.70710678118655 0.5 +vn 4.3297802811775E-17 1 4.3297802811775E-17 +vn -0.5 0.70710678118655 -0.5 +vn -0.70710678118655 1.2246467991474E-16 -0.70710678118655 +vn -0.5 -0.70710678118655 -0.5 +vn -1.2989340843532E-16 -1 -1.2989340843532E-16 +vn 0.5 -0.70710678118655 0.5 +vn 0.92387953251129 0 0.38268343236509 +vn 0.65328148243819 0.70710678118655 0.2705980500731 +vn 5.6571305614385E-17 1 2.3432602026631E-17 +vn -0.65328148243819 0.70710678118655 -0.2705980500731 +vn -0.92387953251129 1.2246467991474E-16 -0.38268343236509 +vn -0.65328148243819 -0.70710678118655 -0.2705980500731 +vn -1.6971391684316E-16 -1 -7.0297806079894E-17 +vn 0.65328148243819 -0.70710678118655 0.2705980500731 +vn 1 0 6.1232339957368E-17 +vn 0.70710678118655 0.70710678118655 4.3297802811775E-17 +vn 6.1232339957368E-17 1 3.7493994566546E-33 +vn -0.70710678118655 0.70710678118655 -4.3297802811775E-17 +vn -1 1.2246467991474E-16 -6.1232339957368E-17 +vn -0.70710678118655 -0.70710678118655 -4.3297802811775E-17 +vn -1.836970198721E-16 -1 -1.1248198369964E-32 +vn 0.70710678118655 -0.70710678118655 4.3297802811775E-17 +vn 0.92387953251129 0 -0.38268343236509 +vn 0.65328148243819 0.70710678118655 -0.2705980500731 +vn 5.6571305614385E-17 1 -2.3432602026631E-17 +vn -0.65328148243819 0.70710678118655 0.2705980500731 +vn -0.92387953251129 1.2246467991474E-16 0.38268343236509 +vn -0.65328148243819 -0.70710678118655 0.2705980500731 +vn -1.6971391684316E-16 -1 7.0297806079894E-17 +vn 0.65328148243819 -0.70710678118655 -0.2705980500731 +vn 0.70710678118655 0 -0.70710678118655 +vn 0.5 0.70710678118655 -0.5 +vn 4.3297802811775E-17 1 -4.3297802811775E-17 +vn -0.5 0.70710678118655 0.5 +vn -0.70710678118655 1.2246467991474E-16 0.70710678118655 +vn -0.5 -0.70710678118655 0.5 +vn -1.2989340843532E-16 -1 1.2989340843532E-16 +vn 0.5 -0.70710678118655 -0.5 +vn 0.38268343236509 0 -0.92387953251129 +vn 0.2705980500731 0.70710678118655 -0.65328148243819 +vn 2.3432602026631E-17 1 -5.6571305614385E-17 +vn -0.2705980500731 0.70710678118655 0.65328148243819 +vn -0.38268343236509 1.2246467991474E-16 0.92387953251129 +vn -0.2705980500731 -0.70710678118655 0.65328148243819 +vn -7.0297806079894E-17 -1 1.6971391684316E-16 +vn 0.2705980500731 -0.70710678118655 -0.65328148243819 +vn 1.2246467991474E-16 0 -1 +vn 8.6595605623549E-17 0.70710678118655 -0.70710678118655 +vn 7.4987989133093E-33 1 -6.1232339957368E-17 +vn -8.6595605623549E-17 0.70710678118655 0.70710678118655 +vn -1.2246467991474E-16 1.2246467991474E-16 1 +vn -8.6595605623549E-17 -0.70710678118655 0.70710678118655 +vn -2.2496396739928E-32 -1 1.836970198721E-16 +vn 8.6595605623549E-17 -0.70710678118655 -0.70710678118655 +vn -0.38268343236509 0 -0.92387953251129 +vn -0.2705980500731 0.70710678118655 -0.65328148243819 +vn -2.3432602026631E-17 1 -5.6571305614385E-17 +vn 0.2705980500731 0.70710678118655 0.65328148243819 +vn 0.38268343236509 1.2246467991474E-16 0.92387953251129 +vn 0.2705980500731 -0.70710678118655 0.65328148243819 +vn 7.0297806079894E-17 -1 1.6971391684316E-16 +vn -0.2705980500731 -0.70710678118655 -0.65328148243819 +vn -0.70710678118655 0 -0.70710678118655 +vn -0.5 0.70710678118655 -0.5 +vn -4.3297802811775E-17 1 -4.3297802811775E-17 +vn 0.5 0.70710678118655 0.5 +vn 0.70710678118655 1.2246467991474E-16 0.70710678118655 +vn 0.5 -0.70710678118655 0.5 +vn 1.2989340843532E-16 -1 1.2989340843532E-16 +vn -0.5 -0.70710678118655 -0.5 +vn -0.92387953251129 0 -0.38268343236509 +vn -0.65328148243819 0.70710678118655 -0.2705980500731 +vn -5.6571305614385E-17 1 -2.3432602026632E-17 +vn 0.65328148243819 0.70710678118655 0.2705980500731 +vn 0.92387953251129 1.2246467991474E-16 0.38268343236509 +vn 0.65328148243819 -0.70710678118655 0.2705980500731 +vn 1.6971391684316E-16 -1 7.0297806079895E-17 +vn -0.65328148243819 -0.70710678118655 -0.2705980500731 +vn -1 0 -1.836970198721E-16 +vn -0.70710678118655 0.70710678118655 -1.2989340843532E-16 +vn -6.1232339957368E-17 1 -1.1248198369964E-32 +vn 0.70710678118655 0.70710678118655 1.2989340843532E-16 +vn 1 1.2246467991474E-16 1.836970198721E-16 +vn 0.70710678118655 -0.70710678118655 1.2989340843532E-16 +vn 1.836970198721E-16 -1 3.3744595109892E-32 +vn -0.70710678118655 -0.70710678118655 -1.2989340843532E-16 +vn -0.92387953251129 0 0.38268343236509 +vn -0.65328148243819 0.70710678118655 0.2705980500731 +vn -5.6571305614385E-17 1 2.3432602026632E-17 +vn 0.65328148243819 0.70710678118655 -0.2705980500731 +vn 0.92387953251129 1.2246467991474E-16 -0.38268343236509 +vn 0.65328148243819 -0.70710678118655 -0.2705980500731 +vn 1.6971391684316E-16 -1 -7.0297806079895E-17 +vn -0.65328148243819 -0.70710678118655 0.2705980500731 +vn -0.70710678118655 0 0.70710678118655 +vn -0.5 0.70710678118655 0.5 +vn -4.3297802811775E-17 1 4.3297802811775E-17 +vn 0.5 0.70710678118655 -0.5 +vn 0.70710678118655 1.2246467991474E-16 -0.70710678118655 +vn 0.5 -0.70710678118655 -0.5 +vn 1.2989340843532E-16 -1 -1.2989340843532E-16 +vn -0.5 -0.70710678118655 0.5 +vn -0.38268343236509 0 0.92387953251129 +vn -0.2705980500731 0.70710678118655 0.65328148243819 +vn -2.3432602026632E-17 1 5.6571305614385E-17 +vn 0.2705980500731 0.70710678118655 -0.65328148243819 +vn 0.38268343236509 1.2246467991474E-16 -0.92387953251129 +vn 0.2705980500731 -0.70710678118655 -0.65328148243819 +vn 7.0297806079895E-17 -1 -1.6971391684316E-16 +vn -0.2705980500731 -0.70710678118655 0.65328148243819 +f 1//1 9//9 10//10 2//2 +f 2//2 10//10 11//11 3//3 +f 3//3 11//11 12//12 4//4 +f 4//4 12//12 13//13 5//5 +f 5//5 13//13 14//14 6//6 +f 6//6 14//14 15//15 7//7 +f 7//7 15//15 16//16 8//8 +f 8//8 16//16 9//9 1//1 +f 9//9 17//17 18//18 10//10 +f 10//10 18//18 19//19 11//11 +f 11//11 19//19 20//20 12//12 +f 12//12 20//20 21//21 13//13 +f 13//13 21//21 22//22 14//14 +f 14//14 22//22 23//23 15//15 +f 15//15 23//23 24//24 16//16 +f 16//16 24//24 17//17 9//9 +f 17//17 25//25 26//26 18//18 +f 18//18 26//26 27//27 19//19 +f 19//19 27//27 28//28 20//20 +f 20//20 28//28 29//29 21//21 +f 21//21 29//29 30//30 22//22 +f 22//22 30//30 31//31 23//23 +f 23//23 31//31 32//32 24//24 +f 24//24 32//32 25//25 17//17 +f 25//25 33//33 34//34 26//26 +f 26//26 34//34 35//35 27//27 +f 27//27 35//35 36//36 28//28 +f 28//28 36//36 37//37 29//29 +f 29//29 37//37 38//38 30//30 +f 30//30 38//38 39//39 31//31 +f 31//31 39//39 40//40 32//32 +f 32//32 40//40 33//33 25//25 +f 33//33 41//41 42//42 34//34 +f 34//34 42//42 43//43 35//35 +f 35//35 43//43 44//44 36//36 +f 36//36 44//44 45//45 37//37 +f 37//37 45//45 46//46 38//38 +f 38//38 46//46 47//47 39//39 +f 39//39 47//47 48//48 40//40 +f 40//40 48//48 41//41 33//33 +f 41//41 49//49 50//50 42//42 +f 42//42 50//50 51//51 43//43 +f 43//43 51//51 52//52 44//44 +f 44//44 52//52 53//53 45//45 +f 45//45 53//53 54//54 46//46 +f 46//46 54//54 55//55 47//47 +f 47//47 55//55 56//56 48//48 +f 48//48 56//56 49//49 41//41 +f 49//49 57//57 58//58 50//50 +f 50//50 58//58 59//59 51//51 +f 51//51 59//59 60//60 52//52 +f 52//52 60//60 61//61 53//53 +f 53//53 61//61 62//62 54//54 +f 54//54 62//62 63//63 55//55 +f 55//55 63//63 64//64 56//56 +f 56//56 64//64 57//57 49//49 +f 57//57 65//65 66//66 58//58 +f 58//58 66//66 67//67 59//59 +f 59//59 67//67 68//68 60//60 +f 60//60 68//68 69//69 61//61 +f 61//61 69//69 70//70 62//62 +f 62//62 70//70 71//71 63//63 +f 63//63 71//71 72//72 64//64 +f 64//64 72//72 65//65 57//57 +f 65//65 73//73 74//74 66//66 +f 66//66 74//74 75//75 67//67 +f 67//67 75//75 76//76 68//68 +f 68//68 76//76 77//77 69//69 +f 69//69 77//77 78//78 70//70 +f 70//70 78//78 79//79 71//71 +f 71//71 79//79 80//80 72//72 +f 72//72 80//80 73//73 65//65 +f 73//73 81//81 82//82 74//74 +f 74//74 82//82 83//83 75//75 +f 75//75 83//83 84//84 76//76 +f 76//76 84//84 85//85 77//77 +f 77//77 85//85 86//86 78//78 +f 78//78 86//86 87//87 79//79 +f 79//79 87//87 88//88 80//80 +f 80//80 88//88 81//81 73//73 +f 81//81 89//89 90//90 82//82 +f 82//82 90//90 91//91 83//83 +f 83//83 91//91 92//92 84//84 +f 84//84 92//92 93//93 85//85 +f 85//85 93//93 94//94 86//86 +f 86//86 94//94 95//95 87//87 +f 87//87 95//95 96//96 88//88 +f 88//88 96//96 89//89 81//81 +f 89//89 97//97 98//98 90//90 +f 90//90 98//98 99//99 91//91 +f 91//91 99//99 100//100 92//92 +f 92//92 100//100 101//101 93//93 +f 93//93 101//101 102//102 94//94 +f 94//94 102//102 103//103 95//95 +f 95//95 103//103 104//104 96//96 +f 96//96 104//104 97//97 89//89 +f 97//97 105//105 106//106 98//98 +f 98//98 106//106 107//107 99//99 +f 99//99 107//107 108//108 100//100 +f 100//100 108//108 109//109 101//101 +f 101//101 109//109 110//110 102//102 +f 102//102 110//110 111//111 103//103 +f 103//103 111//111 112//112 104//104 +f 104//104 112//112 105//105 97//97 +f 105//105 113//113 114//114 106//106 +f 106//106 114//114 115//115 107//107 +f 107//107 115//115 116//116 108//108 +f 108//108 116//116 117//117 109//109 +f 109//109 117//117 118//118 110//110 +f 110//110 118//118 119//119 111//111 +f 111//111 119//119 120//120 112//112 +f 112//112 120//120 113//113 105//105 +f 113//113 121//121 122//122 114//114 +f 114//114 122//122 123//123 115//115 +f 115//115 123//123 124//124 116//116 +f 116//116 124//124 125//125 117//117 +f 117//117 125//125 126//126 118//118 +f 118//118 126//126 127//127 119//119 +f 119//119 127//127 128//128 120//120 +f 120//120 128//128 121//121 113//113 +f 121//121 1//1 2//2 122//122 +f 122//122 2//2 3//3 123//123 +f 123//123 3//3 4//4 124//124 +f 124//124 4//4 5//5 125//125 +f 125//125 5//5 6//6 126//126 +f 126//126 6//6 7//7 127//127 +f 127//127 7//7 8//8 128//128 +f 128//128 8//8 1//1 121//121 diff --git a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/CholmodSolverProxy.h b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/CholmodSolverProxy.h new file mode 100644 index 00000000000..701b709b0aa --- /dev/null +++ b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/CholmodSolverProxy.h @@ -0,0 +1,151 @@ +/****************************************************************************** +* 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 + +// This is a private header. Only include from .cpp files that also include +// . + +#include +#include + +namespace sofacholmod +{ + +/** + * Supernodal LL^T solver derived directly from CholmodBase (not from + * CholmodSupernodalLLT) so that the protected m_cholmodFactor and m_cholmod + * members remain accessible. CholmodSupernodalLLT makes them private via + * using-declarations, preventing access from further derived classes. + */ +template +class AccessibleCholmodLLT + : public Eigen::CholmodBase> +{ + using Base = Eigen::CholmodBase; + +public: + using typename Base::StorageIndex; + + AccessibleCholmodLLT() : Base() + { + this->m_cholmod.final_asis = 1; + this->m_cholmod.supernodal = CHOLMOD_SUPERNODAL; + } + + cholmod_factor* cholmodFactor() { return this->m_cholmodFactor; } + cholmod_common& cholmodCommon() { return this->m_cholmod; } +}; + +/** + * CHOLMOD-specific proxy that replaces the generic EigenSolverWrapper for + * CHOLMOD solvers. Exposes the raw CHOLMOD factor and common, enabling + * optimized partial solves (e.g. for compliance matrix computation). + */ +class CholmodSolverProxy : public sofa::component::linearsolver::direct::BaseEigenSolverProxy +{ + using CholmodEigenMatrix = Eigen::SparseMatrix; + + AccessibleCholmodLLT m_solver; + + // Persistent scratch reused by forwardSolveLP() across calls, so that + // cholmod_solve2 recycles these buffers instead of allocating/freeing a new + // dense at every constraint-solve step (matters for many small systems). + cholmod_dense* m_permutedB = nullptr; // P * B + cholmod_dense* m_Z = nullptr; // L^{-1} * P * B (the returned result) + cholmod_dense* m_Y = nullptr; // cholmod_solve2 workspace + cholmod_dense* m_E = nullptr; // cholmod_solve2 workspace + +public: + + ~CholmodSolverProxy() override + { + cholmod_common& c = m_solver.cholmodCommon(); + cholmod_free_dense(&m_permutedB, &c); + cholmod_free_dense(&m_Z, &c); + cholmod_free_dense(&m_Y, &c); + cholmod_free_dense(&m_E, &c); + } + + [[nodiscard]] Eigen::ComputationInfo info() const override + { + return m_solver.info(); + } + + void solve(const EigenVectorXdMap& /*b*/, EigenVectorXdMap& /*x*/) const override + { + msg_error("CholmodSolverProxy") << "CHOLMOD only supports double precision"; + } + + void solve(const EigenVectorXdMap& b, EigenVectorXdMap& x) const override + { + x = m_solver.solve(b); + } + + void analyzePattern(const BaseEigenSolverProxy::EigenSparseMatrixMap& /*a*/) override + { + msg_error("CholmodSolverProxy") << "CHOLMOD only supports double precision"; + } + + void analyzePattern(const BaseEigenSolverProxy::EigenSparseMatrixMap& a) override + { + m_solver.analyzePattern(a); + } + + void factorize(const BaseEigenSolverProxy::EigenSparseMatrixMap& /*a*/) override + { + msg_error("CholmodSolverProxy") << "CHOLMOD only supports double precision"; + } + + void factorize(const BaseEigenSolverProxy::EigenSparseMatrixMap& a) override + { + m_solver.factorize(a); + } + + // CHOLMOD-specific accessors for partial solves + cholmod_factor* cholmodFactor() { return m_solver.cholmodFactor(); } + cholmod_common& cholmodCommon() { return m_solver.cholmodCommon(); } + + /// Solve Z = L^{-1} * P * B for a dense right-hand side B, using + /// cholmod_solve2 so the result and workspace buffers are reused across + /// calls (no per-call allocation). The returned dense is owned by the proxy + /// and stays valid until the next call to this method or destruction. + /// Returns nullptr on failure. + cholmod_dense* forwardSolveLP(cholmod_dense* B) + { + cholmod_common& c = m_solver.cholmodCommon(); + cholmod_factor* f = m_solver.cholmodFactor(); + + if (!cholmod_solve2(CHOLMOD_P, f, B, nullptr, + &m_permutedB, nullptr, &m_Y, &m_E, &c)) + { + return nullptr; + } + if (!cholmod_solve2(CHOLMOD_L, f, m_permutedB, nullptr, + &m_Z, nullptr, &m_Y, &m_E, &c)) + { + return nullptr; + } + return m_Z; + } +}; + +} // namespace sofacholmod diff --git a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp index 487e325081d..f4c8813e1b4 100644 --- a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp +++ b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.cpp @@ -25,6 +25,7 @@ #include #include +#include #include #include @@ -44,6 +45,34 @@ namespace sofacholmod namespace { +/// Resolve a symbol from the BLAS backend already loaded in the process (the +/// same one CHOLMOD uses). Returns nullptr if not found. No link-time dependency +/// on BLAS is introduced. +void* resolveBlasSymbol(const char* name) +{ +#if defined(_WIN32) + // The BLAS DLL is already loaded (pulled in by CHOLMOD); look it up by name. + static const char* const moduleNames[] = { + "libopenblas", "openblas", "libblas", "mkl_rt" + }; + for (const char* moduleName : moduleNames) + { + if (HMODULE module = GetModuleHandleA(moduleName)) + { + if (auto* p = reinterpret_cast(GetProcAddress(module, name))) + { + return p; + } + } + } + return nullptr; +#else + // The BLAS symbols are already loaded in the process (via CHOLMOD), + // so search the global symbol table. + return dlsym(RTLD_DEFAULT, name); +#endif +} + /// Try to set the number of threads used by the BLAS backend at runtime. /// /// CHOLMOD's supernodal factorization delegates to dense BLAS3 kernels, so the @@ -64,55 +93,64 @@ bool trySetBlasNumThreads(int numThreads) "MKL_Set_Num_Threads", // Intel MKL }; -#if defined(_WIN32) - // The BLAS DLL is already loaded (pulled in by CHOLMOD); look it up by name. - static const char* const moduleNames[] = { - "libopenblas", "openblas", "libblas", "mkl_rt" - }; - for (const char* moduleName : moduleNames) - { - HMODULE module = GetModuleHandleA(moduleName); - if (!module) - { - continue; - } - for (const char* name : candidates) - { - if (auto fn = reinterpret_cast( - reinterpret_cast(GetProcAddress(module, name)))) - { - fn(numThreads); - return true; - } - } - } - return false; -#else - // The BLAS symbols are already loaded in the process (via CHOLMOD), - // so search the global symbol table. for (const char* name : candidates) { - if (auto fn = reinterpret_cast(dlsym(RTLD_DEFAULT, name))) + if (auto fn = reinterpret_cast(resolveBlasSymbol(name))) { fn(numThreads); return true; } } return false; -#endif +} + +/// Compute the lower triangle of W = fact * Z^T * Z, where Z is an n x m +/// column-major matrix with leading dimension ldZ (>= n). Uses the BLAS +/// symmetric rank-k update (dsyrk): it computes only the triangle we need (half +/// the flops of a full product) and is multi-threaded on OpenBLAS/MKL according +/// to the BLAS thread count (see the numThreads Data). Falls back to an Eigen +/// rank update if dsyrk is unavailable. +void computeLowerZtZ(const double* Z, Eigen::Index n, Eigen::Index m, Eigen::Index ldZ, + double fact, Eigen::MatrixXd& W) +{ + // CBLAS enum values (avoids depending on a cblas.h header being present). + constexpr int CblasColMajor = 102; + constexpr int CblasLower = 122; + constexpr int CblasTrans = 112; + + using DsyrkFn = void (*)(int, int, int, int, int, double, + const double*, int, double, double*, int); + static DsyrkFn dsyrk = reinterpret_cast(resolveBlasSymbol("cblas_dsyrk")); + + if (dsyrk) + { + // W := fact * Z^T * Z (lower triangle), with Z treated as A (n x m). + dsyrk(CblasColMajor, CblasLower, CblasTrans, + static_cast(m), static_cast(n), fact, + Z, static_cast(ldZ), 0.0, W.data(), static_cast(m)); + } + else + { + Eigen::Map> Zmap(Z, n, m, Eigen::OuterStride<>(ldZ)); + W.setZero(); + W.selfadjointView().rankUpdate(Zmap.transpose(), fact); + } } } // namespace template EigenCholmodSupernodalLLT::EigenCholmodSupernodalLLT() - : d_numThreads(initData(&d_numThreads, 0, + : d_numThreads(initData(&d_numThreads, 1, "numThreads", "Number of threads used by the BLAS backend of CHOLMOD's supernodal " - "factorization. A value <= 0 (default) keeps the BLAS default " - "(OPENBLAS_NUM_THREADS / OMP_NUM_THREADS). A small value (1-4) is often " - "faster on medium-sized systems by avoiding thread oversubscription. " - "Only effective with OpenBLAS or MKL.")) + "factorization. Default is 1: this is the fastest setting for the vast " + "majority of SOFA scenes, and it avoids catastrophic thread " + "oversubscription when several solvers are factorized in parallel (e.g. " + "with parallelODESolving). Increase it only for a single, very large " + "standalone system. A value <= 0 keeps the BLAS default " + "(OPENBLAS_NUM_THREADS / OMP_NUM_THREADS). Only effective with OpenBLAS " + "or MKL.")) { } @@ -130,6 +168,116 @@ void EigenCholmodSupernodalLLT::reinit() applyBlasNumThreads(); } +template +void EigenCholmodSupernodalLLT::invert(Matrix& A) +{ + if (const int nthreads = d_numThreads.getValue(); nthreads > 0) + { + if (auto* proxy = dynamic_cast(this->m_solver.get())) + { + proxy->cholmodCommon().nthreads_max = nthreads; + } + } + + Inherit1::invert(A); +} + + +// Adds the compliance block W = fact * J * A^{-1} * J^T into 'result'. +// +// CHOLMOD factorizes P*A*P^T = L*L^T, so A^{-1} = P^T * L^{-T} * L^{-1} * P and +// J * A^{-1} * J^T = (L^{-1} * P * J^T)^T * (L^{-1} * P * J^T) = Z^T * Z. +// The whole block is therefore one triangular forward-solve (with the m columns +// of J^T as right-hand sides) followed by one symmetric product Z^T*Z, avoiding +// a full solve per constraint row. +template +bool EigenCholmodSupernodalLLT::addJMInvJtLocal( + Matrix* M, ResMatrixType* result, const JMatrixType* J, SReal fact) +{ + SOFA_UNUSED(M); + + if (!this->isComponentStateValid()) + { + return true; + } + + // Access the raw CHOLMOD factor via our proxy. The generic EigenSolverWrapper + // does not expose it, so fall back to the base implementation if, for any + // reason, the registered solver is not our CholmodSolverProxy. + auto* proxy = dynamic_cast(this->m_solver.get()); + if (!proxy) + { + return Inherit1::addJMInvJtLocal(M, result, J, fact); + } + + if (J->rowSize() == 0) + { + return true; + } + + // Global row index of each (sparse) constraint row, used to scatter W back. + std::vector jLocal2Global; + jLocal2Global.reserve(J->rowSize()); + for (auto jit = J->begin(), jitend = J->end(); jit != jitend; ++jit) + { + jLocal2Global.push_back(jit->first); + } + + if (jLocal2Global.empty()) + { + return true; + } + + const auto m = static_cast(jLocal2Global.size()); // constraint rows + const auto n = static_cast(J->colSize()); // system size + + // 1. Expand the sparse J into a dense J^T (n x m, column-major as CHOLMOD expects). + // m_Jt is a reused buffer: setZero resizes only when the shape changes. + m_Jt.setZero(n, m); + { + Eigen::Index col = 0; + for (auto jit = J->begin(), jitend = J->end(); jit != jitend; ++jit, ++col) + { + for (auto it = jit->second.begin(), itend = jit->second.end(); it != itend; ++it) + { + m_Jt(it->first, col) = it->second; + } + } + } + + // 2. + 3. Z = L^{-1} * P * J^T, all m columns at once (permutation then supernodal + // BLAS3 forward-solve). The proxy reuses its workspace across calls, so Z_cd is + // owned by the proxy and must not be freed here. + cholmod_dense Jt_cd = Eigen::viewAsCholmod(m_Jt); + cholmod_dense* Z_cd = proxy->forwardSolveLP(&Jt_cd); + if (!Z_cd) + { + msg_error() << "CHOLMOD forward solve failed"; + return false; + } + + // 4. W = fact * Z^T * Z. W is symmetric, so only its lower triangle is computed. + Eigen::MatrixXd W(m, m); + computeLowerZtZ(static_cast(Z_cd->x), n, m, + static_cast(Z_cd->d), fact, W); + + // 5. Scatter the lower triangle (j >= i) into 'result', mirroring across the diagonal. + for (Eigen::Index i = 0; i < m; ++i) + { + for (Eigen::Index j = i; j < m; ++j) + { + const double value = W(j, i); + result->add(jLocal2Global[j], jLocal2Global[i], value); + if (i != j) + { + result->add(jLocal2Global[i], jLocal2Global[j], value); + } + } + } + + return true; +} + template void EigenCholmodSupernodalLLT::applyBlasNumThreads() { diff --git a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.h b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.h index a418e9e4cfd..4e6045a5e65 100644 --- a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.h +++ b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/EigenCholmodSupernodalLLT.h @@ -74,18 +74,28 @@ class EigenCholmodSupernodalLLT SOFA_CLASS(SOFA_TEMPLATE(EigenCholmodSupernodalLLT, TBlockType), SOFA_TEMPLATE2(sofa::component::linearsolver::direct::EigenDirectSparseSolver, TBlockType, MainCholmodSupernodalLLTFactory)); + using JMatrixType = typename Inherit1::JMatrixType; + using ResMatrixType = typename Inherit1::ResMatrixType; + void init() override; void reinit() override; + void invert(Matrix& A) override; + + bool addJMInvJtLocal(Matrix* M, ResMatrixType* result, + const JMatrixType* J, SReal fact) override; /// Number of threads the underlying BLAS (used by CHOLMOD's supernodal /// factorization) is allowed to use. - /// A value <= 0 (default) leaves the BLAS default untouched, i.e. controlled - /// by the OPENBLAS_NUM_THREADS / OMP_NUM_THREADS environment variables. - /// For the medium-sized systems typically solved here, a small value (1-4) - /// is often faster than the default (= number of cores) because it avoids - /// thread oversubscription. Only effective with OpenBLAS or MKL; ignored - /// with BLAS backends that expose no runtime thread-control API (e.g. Apple - /// Accelerate), where the environment variables must be used instead. + /// Default is 1, which is the fastest setting for the vast majority of SOFA + /// scenes: it avoids thread oversubscription, which is catastrophic when + /// several solvers are factorized in parallel (e.g. with parallelODESolving) + /// and still optimal for a single medium-sized system. Increase it only for + /// a single, very large standalone system. A value <= 0 leaves the BLAS + /// default untouched, i.e. controlled by the OPENBLAS_NUM_THREADS / + /// OMP_NUM_THREADS environment variables. Only effective with OpenBLAS or + /// MKL; ignored with BLAS backends that expose no runtime thread-control API + /// (e.g. Apple Accelerate), where the environment variables must be used + /// instead. sofa::core::objectmodel::Data d_numThreads; protected: @@ -94,6 +104,11 @@ class EigenCholmodSupernodalLLT /// Apply d_numThreads to the underlying BLAS backend if it exposes a runtime /// thread-count API (OpenBLAS/MKL). No-op if d_numThreads <= 0. void applyBlasNumThreads(); + + /// Reused dense J^T buffer for addJMInvJtLocal, to avoid reallocating at every + /// constraint-solve step. Not thread-safe, but each solver instance handles a + /// single object whose compliance is built sequentially. + Eigen::MatrixXd m_Jt; }; #ifndef SOFACHOLMOD_EIGENCHOLMODSUPERNODALLLT_CPP diff --git a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/init.cpp b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/init.cpp index 211c522908b..b8061a814bd 100644 --- a/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/init.cpp +++ b/applications/plugins/SofaCHOLMOD/src/SofaCHOLMOD/init.cpp @@ -25,6 +25,7 @@ #include #include +#include #include namespace sofacholmod @@ -36,9 +37,11 @@ template void MainCholmodSupernodalLLTFactory::registerSolver(const std::string& orderingMethodName) { std::lock_guard lock(s_mutex); - // CHOLMOD uses its own internal ordering; the OrderingMethodType is unused - getFactory().registerType< - Eigen::CholmodSupernodalLLT> >(orderingMethodName); + // CHOLMOD uses its own internal ordering; the OrderingMethodType is unused. + // Register our own proxy (not the generic EigenSolverWrapper) so that the + // solver can access the raw CHOLMOD factor for the optimized + // addJMInvJtLocal() compliance-matrix computation. + getFactory().registerProxyType(orderingMethodName); } From 7faa62a3efde915ab936c07660e9ef32af5a3827 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Mon, 13 Jul 2026 08:46:56 +0900 Subject: [PATCH 10/10] add tests and update comments about numthreads --- .../plugins/SofaCHOLMOD/CMakeLists.txt | 6 + applications/plugins/SofaCHOLMOD/README.md | 61 ++--- .../FEMBAR_EigenCholmodSupernodalLLT.scn | 4 + .../plugins/SofaCHOLMOD/tests/CMakeLists.txt | 15 ++ .../tests/EigenCholmodSupernodalLLT_test.cpp | 213 ++++++++++++++++++ 5 files changed, 273 insertions(+), 26 deletions(-) create mode 100644 applications/plugins/SofaCHOLMOD/tests/CMakeLists.txt create mode 100644 applications/plugins/SofaCHOLMOD/tests/EigenCholmodSupernodalLLT_test.cpp diff --git a/applications/plugins/SofaCHOLMOD/CMakeLists.txt b/applications/plugins/SofaCHOLMOD/CMakeLists.txt index 62f7466f734..4bd7c8cb3a9 100644 --- a/applications/plugins/SofaCHOLMOD/CMakeLists.txt +++ b/applications/plugins/SofaCHOLMOD/CMakeLists.txt @@ -52,3 +52,9 @@ sofa_create_package_with_targets( INCLUDE_SOURCE_DIR "src" INCLUDE_INSTALL_DIR "${PROJECT_NAME}" ) + +include(CMakeDependentOption) +cmake_dependent_option(SOFACHOLMOD_BUILD_TESTS "Compile the automatic tests" ON "SOFA_BUILD_TESTS OR NOT DEFINED SOFA_BUILD_TESTS" OFF) +if(SOFACHOLMOD_BUILD_TESTS) + add_subdirectory(tests) +endif() diff --git a/applications/plugins/SofaCHOLMOD/README.md b/applications/plugins/SofaCHOLMOD/README.md index bf5311a2313..a8ed646fe90 100644 --- a/applications/plugins/SofaCHOLMOD/README.md +++ b/applications/plugins/SofaCHOLMOD/README.md @@ -41,7 +41,7 @@ of the common linear-solver API), but its choice is **ignored** by this solver. | Data | Type | Default | Description | |--------------|-------|---------|-------------| -| `numThreads` | `int` | `1` | Number of threads the underlying BLAS backend may use for the factorization. The default of `1` is the fastest setting for the vast majority of SOFA scenes and avoids thread oversubscription (see [Threading](#threading)). Increase it only for a single, very large standalone system. A value `<= 0` leaves the BLAS default untouched (controlled by the `OPENBLAS_NUM_THREADS` / `OMP_NUM_THREADS` environment variables). Only effective with OpenBLAS or MKL. | +| `numThreads` | `int` | `1` | Number of threads the underlying BLAS backend may use for the factorization. The default of `1` is a safe choice that avoids thread oversubscription when several solvers run concurrently (see [Threading](#threading)). For a **single large standalone system**, a moderate value (roughly half the physical cores) is faster and should be tuned per machine. A value `<= 0` leaves the BLAS default untouched (controlled by the `OPENBLAS_NUM_THREADS` / `OMP_NUM_THREADS` environment variables). Only effective with OpenBLAS or MKL. | ## Usage @@ -109,39 +109,48 @@ CHOLMOD**, not by SOFA. Keep this in mind when benchmarking: ### Threading An optimized multithreaded BLAS (e.g. OpenBLAS-pthread) uses *all* cores by -default. For the medium-sized systems typical in SOFA this **oversubscribes** -threads: the launch/synchronization overhead makes the solver *slower* than with -a single thread, and the problem becomes catastrophic when several solvers are -factorized in parallel (e.g. multiple objects with `parallelODESolving="true"`), -where each solver would otherwise spawn one BLAS pool per core. - -For these reasons `numThreads` **defaults to `1`**, which is the fastest setting -in nearly all cases. Increase it only for a single, very large standalone -system. - -Illustrative measurements on a 24-core Linux box, OpenBLAS backend: - -Single medium system — `examples/FEMBAR_EigenCholmodSupernodalLLT.scn` -(10×10×50 grid, ~15k DOF): - -| Solver / configuration | ms / step | Speedup | -|-------------------------------------------------|-----------|---------| -| `SparseLDLSolver` | ~762 | 1.0x | -| `EigenCholmodSupernodalLLT`, reference BLAS | ~346 | 2.2x | -| `EigenCholmodSupernodalLLT`, OpenBLAS, all cores| ~92 | 8.3x | -| `EigenCholmodSupernodalLLT`, OpenBLAS, `numThreads="1"` | ~70 | 10.9x | +default. The right thread count depends entirely on the scene: + +- **A single large system** benefits from multithreaded BLAS, but only up to a + point: there is a sweet spot (roughly half the physical cores), beyond which + launch/synchronization overhead makes it *slower* again — using all cores is + usually the worst choice. +- **Many systems solved concurrently** (multiple objects, especially with + `parallelODESolving="true"`) already saturate the cores through object-level + parallelism. Adding BLAS threads on top **oversubscribes** them and is + catastrophic. + +`numThreads` therefore **defaults to `1`** — the safe choice that never +oversubscribes. For a single large system, raise it to the sweet spot (tune per +machine). The optimum is machine-dependent, so treat the numbers below as +illustrative (24-core Linux box, OpenBLAS backend). + +Single large system — `examples/FEMBAR_EigenCholmodSupernodalLLT.scn` +(10×10×50 grid, ~15k DOF). Effect of `numThreads`: + +| `numThreads` | ms / step | FPS | +|-------------------|-----------|------| +| 1 | ~67 | 14.9 | +| 4 | ~52 | 19.1 | +| **8** (sweet spot)| **~48** | 20.8 | +| 16 | ~64 | 15.6 | +| 24 (all cores) | ~96 | 10.4 | + +For reference, on the same scene `SparseLDLSolver` is ~762 ms/step, and CHOLMOD +on a non-optimized *reference* BLAS is ~346 ms/step — i.e. having an optimized +BLAS matters far more than the exact thread count. Many small systems in parallel — `examples/TorusFall.scn` (10 tori, -`parallelODESolving="true"`): +`parallelODESolving="true"`). Here more BLAS threads only hurt: | Solver / configuration | FPS | Speedup | |-------------------------------------------------|-------|---------| -| `EigenCholmodSupernodalLLT`, OpenBLAS, all cores| ~20 | 0.5x | +| `EigenCholmodSupernodalLLT`, all cores | ~20 | 0.5x | | `SparseLDLSolver` (`parallelInverseProduct`) | ~40 | 1.0x | | `EigenCholmodSupernodalLLT`, `numThreads="1"` | ~108 | 2.7x | -> Multithreaded BLAS only pays off for a single, very large system; whenever -> several solvers run concurrently, keep `numThreads="1"`. +> Rule of thumb: `numThreads="1"` whenever several solvers run concurrently; a +> moderate value (~half the cores) for a single large standalone system. ## Constraint solving (compliance matrix) diff --git a/applications/plugins/SofaCHOLMOD/examples/FEMBAR_EigenCholmodSupernodalLLT.scn b/applications/plugins/SofaCHOLMOD/examples/FEMBAR_EigenCholmodSupernodalLLT.scn index 28ca143592d..0929e3a34f0 100644 --- a/applications/plugins/SofaCHOLMOD/examples/FEMBAR_EigenCholmodSupernodalLLT.scn +++ b/applications/plugins/SofaCHOLMOD/examples/FEMBAR_EigenCholmodSupernodalLLT.scn @@ -5,6 +5,10 @@ + diff --git a/applications/plugins/SofaCHOLMOD/tests/CMakeLists.txt b/applications/plugins/SofaCHOLMOD/tests/CMakeLists.txt new file mode 100644 index 00000000000..e25025a8d8e --- /dev/null +++ b/applications/plugins/SofaCHOLMOD/tests/CMakeLists.txt @@ -0,0 +1,15 @@ +cmake_minimum_required(VERSION 3.22) + +project(SofaCHOLMOD_test) + +set(SOURCE_FILES + EigenCholmodSupernodalLLT_test.cpp +) + +add_executable(${PROJECT_NAME} ${SOURCE_FILES}) +target_link_libraries(${PROJECT_NAME} + Sofa.Testing + SofaCHOLMOD + Sofa.Component.LinearSolver.Ordering + Eigen3::Eigen +) diff --git a/applications/plugins/SofaCHOLMOD/tests/EigenCholmodSupernodalLLT_test.cpp b/applications/plugins/SofaCHOLMOD/tests/EigenCholmodSupernodalLLT_test.cpp new file mode 100644 index 00000000000..4fc1fd710c2 --- /dev/null +++ b/applications/plugins/SofaCHOLMOD/tests/EigenCholmodSupernodalLLT_test.cpp @@ -0,0 +1,213 @@ +/****************************************************************************** +* 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 +#include + +#include +#include +#include +#include +#include + +#include + +namespace sofacholmod +{ + +using sofa::core::objectmodel::New; +using NaturalOrdering = sofa::component::linearsolver::ordering::NaturalOrderingMethod; + +namespace +{ +/// Instantiate a CHOLMOD solver with a NaturalOrderingMethod and init() it. +/// sofacholmod::init() populates the internal solver factory (idempotent). +template +typename Solver::SPtr makeSolver(NaturalOrdering::SPtr& orderingOut) +{ + sofacholmod::init(); + auto solver = New(); + orderingOut = New(); + solver->l_orderingMethod.set(orderingOut.get()); + solver->init(); + return solver; +} +} + +/// The factory must produce a CholmodSolverProxy (not the generic wrapper), +/// otherwise the optimized addJMInvJtLocal path is silently bypassed and the +/// solver falls back to a slow full-solve-per-constraint-row implementation. +TEST(EigenCholmodSupernodalLLT, FactoryRegistersCholmodProxy) +{ + sofacholmod::init(); + + ASSERT_TRUE(MainCholmodSupernodalLLTFactory::hasSolver("Natural")); + + auto* proxy = MainCholmodSupernodalLLTFactory::getSolver("Natural"); + ASSERT_NE(proxy, nullptr); + EXPECT_NE(dynamic_cast(proxy), nullptr) + << "the CHOLMOD factory must produce a CholmodSolverProxy"; + delete proxy; +} + +/// Solve a SPD tridiagonal system and check the residual A*x == b. +TEST(EigenCholmodSupernodalLLT, SolveSPDScalar) +{ + using MatrixType = sofa::linearalgebra::CompressedRowSparseMatrix; + using VectorType = sofa::linearalgebra::FullVector; + using Solver = EigenCholmodSupernodalLLT; + + NaturalOrdering::SPtr ordering; + const Solver::SPtr solver = makeSolver(ordering); + + constexpr int n = 6; + MatrixType A(n, n); + for (int i = 0; i < n; ++i) + { + A.add(i, i, 2.0 + 0.1 * i); // SPD, diagonally dominant + if (i > 0) + { + A.add(i, i - 1, -1.0); + A.add(i - 1, i, -1.0); + } + } + A.compress(); + + VectorType b(n), x(n); + for (int i = 0; i < n; ++i) + { + b[i] = static_cast(i + 1); + } + + solver->invert(A); + solver->solve(A, x, b); + + for (int i = 0; i < n; ++i) + { + SReal Ax = 0; + for (int j = 0; j < n; ++j) + { + Ax += A(i, j) * x[j]; + } + EXPECT_NEAR(Ax, b[i], 1e-9) << "residual mismatch at row " << i; + } +} + +/// Same, using the Mat3x3 block template. +TEST(EigenCholmodSupernodalLLT, SolveSPDBlock3x3) +{ + using Block = sofa::type::Mat<3, 3, SReal>; + using MatrixType = sofa::linearalgebra::CompressedRowSparseMatrix; + using VectorType = sofa::linearalgebra::FullVector; + using Solver = EigenCholmodSupernodalLLT; + + NaturalOrdering::SPtr ordering; + const Solver::SPtr solver = makeSolver(ordering); + + // Two SPD diagonal blocks. + constexpr int n = 6; + MatrixType A(n, n); + Block B0; B0.identity(); + Block B1; B1.clear(); B1[0][0] = 2.0; B1[1][1] = 4.0; B1[2][2] = 8.0; + A.add(0, 0, B0); + A.add(3, 3, B1); + A.compress(); + + VectorType b(n), x(n); + b[0] = 1.0; b[1] = 2.0; b[2] = 3.0; + b[3] = 4.0; b[4] = 8.0; b[5] = 16.0; + + solver->invert(A); + solver->solve(A, x, b); + + const SReal expected[n] = {1.0, 2.0, 3.0, 2.0, 2.0, 2.0}; + for (int i = 0; i < n; ++i) + { + EXPECT_NEAR(x[i], expected[i], 1e-9) << "at index " << i; + } +} + +/// Verify the CHOLMOD-specific addJMInvJtLocal against a dense reference +/// W = fact * J * A^{-1} * J^T. This exercises the CholmodSolverProxy path +/// (forward-solve via cholmod_solve2 + the dsyrk symmetric product). +TEST(EigenCholmodSupernodalLLT, AddJMInvJtMatchesDenseReference) +{ + using MatrixType = sofa::linearalgebra::CompressedRowSparseMatrix; + using Solver = EigenCholmodSupernodalLLT; + + NaturalOrdering::SPtr ordering; + const Solver::SPtr solver = makeSolver(ordering); + + constexpr int n = 6; + MatrixType A(n, n); + Eigen::MatrixXd Adense = Eigen::MatrixXd::Zero(n, n); + const auto addA = [&](int i, int j, double v) { A.add(i, j, v); Adense(i, j) += v; }; + for (int i = 0; i < n; ++i) + { + addA(i, i, 2.0 + 0.1 * i); + if (i > 0) + { + addA(i, i - 1, -1.0); + addA(i - 1, i, -1.0); + } + } + A.compress(); + + solver->invert(A); + // A directly-instantiated solver is not necessarily left in the Valid state + // by init(); addJMInvJtLocal early-returns unless the component is valid. + solver->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); + + // Constraint Jacobian J (m x n), sparse. + constexpr int m = 3; + Solver::JMatrixType J; + J.resize(m, n); + Eigen::MatrixXd Jdense = Eigen::MatrixXd::Zero(m, n); + const auto setJ = [&](int i, int j, double v) { J.set(i, j, v); Jdense(i, j) = v; }; + setJ(0, 0, 1.0); setJ(0, 2, -0.5); + setJ(1, 1, 2.0); setJ(1, 4, 1.0); + setJ(2, 3, 1.0); setJ(2, 5, -1.0); + + const SReal fact = 0.75; + + sofa::linearalgebra::FullMatrix W; + W.resize(m, m); + + const bool ok = solver->addJMInvJtLocal(&A, &W, &J, fact); + ASSERT_TRUE(ok); + + const Eigen::MatrixXd Wref = fact * Jdense * Adense.inverse() * Jdense.transpose(); + + for (int i = 0; i < m; ++i) + { + for (int j = 0; j < m; ++j) + { + EXPECT_NEAR(W.element(i, j), Wref(i, j), 1e-8) + << "compliance mismatch at (" << i << "," << j << ")"; + } + } +} + +} // namespace sofacholmod