[Core][LinearSolver] Introduction of ordering method components (#4477)

* Introduction of ordering method components

* Give names to ordering methods

* Use OrderingMethod to select the solver

* Introduce EigenSolver factory

* cleanup

* Prefer metis when available

* deprecate traits

* documentation

* Fix by adding this

* Fix initialization of static variable

* Deprecate the Data 'ordering'

* Remove usage of ordering Data

* Initialize solver otherwise ordering method is not crated

Sofa/Component/LinearSolver/CMakeLists.txt

@@ -4,7 +4,11 @@ project(Sofa.Component.LinearSolver LANGUAGES CXX)
 set(SOFACOMPONENTLINEARSOLVER_SOURCE_DIR "src/sofa/component/linearsolver")
 
 sofa_add_subdirectory_modules(SOFACOMPONENTLINEARSOLVER_TARGETS
-    DIRECTORIES Iterative Direct Preconditioner
+    DIRECTORIES
+        Iterative
+        Ordering
+        Direct
+        Preconditioner
 )
 
 set(HEADER_FILES

Sofa/Component/LinearSolver/Direct/CMakeLists.txt

@@ -14,11 +14,13 @@ set(HEADER_FILES
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/CholeskySolver.inl
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenDirectSparseSolver.h
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenDirectSparseSolver.inl
+    ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenSolverFactory.h
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenSimplicialLDLT.h
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenSimplicialLLT.h
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenSparseLU.h
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenSparseQR.h
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/MatrixLinearSystem[BTDMatrix].h
+    ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/MetisOrderingMethod.h
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/PrecomputedLinearSolver.h
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/PrecomputedLinearSolver.inl
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/SVDLinearSolver.h
@@ -45,9 +47,11 @@ set(SOURCE_FILES
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/CholeskySolver.cpp
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenSimplicialLDLT.cpp
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenSimplicialLLT.cpp
+    ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenSolverFactory.cpp
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenSparseLU.cpp
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/EigenSparseQR.cpp
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/MatrixLinearSystem[BTDMatrix].cpp
+    ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/MetisOrderingMethod.cpp
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/PrecomputedLinearSolver.cpp
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/SVDLinearSolver.cpp
     ${SOFACOMPONENTLINEARSOLVERDIRECT_SOURCE_DIR}/SparseCommon.cpp
@@ -68,9 +72,10 @@ sofa_find_package(Threads REQUIRED)
 
 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)
 
 add_library(${PROJECT_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES} ${WRAPPER_FILES})
-target_link_libraries(${PROJECT_NAME} PUBLIC Sofa.Simulation.Core Sofa.Component.LinearSolver.Iterative)
+target_link_libraries(${PROJECT_NAME} PUBLIC Sofa.Simulation.Core Sofa.Component.LinearSolver.Iterative Sofa.Component.LinearSolver.Ordering)
 target_link_libraries(${PROJECT_NAME} PUBLIC metis)
 target_link_libraries(${PROJECT_NAME} PUBLIC Threads::Threads)
 

Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenDirectSparseSolver.h

@@ -22,8 +22,10 @@
 #pragma once
 #include <sofa/component/linearsolver/direct/config.h>
 #include <sofa/component/linearsolver/iterative/MatrixLinearSolver.h>
+#include <sofa/component/linearsolver/ordering/OrderingMethodAccessor.h>
 #include <variant>
 #include <Eigen/SparseCore>
+#include <sofa/component/linearsolver/direct/EigenSolverFactory.h>
 
 #include <sofa/helper/OptionsGroup.h>
 
@@ -33,24 +35,24 @@ namespace sofa::component::linearsolver::direct
 /**
  * Base class for all Eigen based direct sparse solvers
  */
-template<class TBlockType, class EigenSolver>
+template<class TBlockType, class TEigenSolverFactory>
 class EigenDirectSparseSolver
-    : public sofa::component::linearsolver::MatrixLinearSolver<
-        sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType>,
-        sofa::linearalgebra::FullVector<typename sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType>::Real> >
+    : public ordering::OrderingMethodAccessor<
+        sofa::component::linearsolver::MatrixLinearSolver<
+            sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType>,
+            sofa::linearalgebra::FullVector<typename sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType>::Real>
+        >
+    >
 {
 public:
-    typedef sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType> Matrix;
+    using Matrix = sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType>;
     using Real = typename Matrix::Real;
-    typedef sofa::linearalgebra::FullVector<Real> Vector;
+    using Vector = sofa::linearalgebra::FullVector<Real>;
 
-    SOFA_ABSTRACT_CLASS(SOFA_TEMPLATE2(EigenDirectSparseSolver, TBlockType, EigenSolver),
-        SOFA_TEMPLATE2(sofa::component::linearsolver::MatrixLinearSolver, Matrix, Vector));
+    using EigenSolverFactory = TEigenSolverFactory;
 
-    using NaturalOrderSolver = typename EigenSolver::NaturalOrderSolver;
-    using AMDOrderSolver     = typename EigenSolver::AMDOrderSolver;
-    using COLAMDOrderSolver  = typename EigenSolver::COLAMDOrderSolver;
-    using MetisOrderSolver   = typename EigenSolver::MetisOrderSolver;
+    SOFA_ABSTRACT_CLASS(SOFA_TEMPLATE2(EigenDirectSparseSolver, TBlockType, EigenSolverFactory),
+        SOFA_TEMPLATE(ordering::OrderingMethodAccessor, SOFA_TEMPLATE2(sofa::component::linearsolver::MatrixLinearSolver, Matrix, Vector)));
 
     ~EigenDirectSparseSolver() override = default;
 
@@ -66,12 +68,12 @@ class EigenDirectSparseSolver
 
 protected:
 
-    sofa::core::objectmodel::Data<sofa::helper::OptionsGroup> d_orderingMethod;
-    unsigned int m_selectedOrderingMethod { std::numeric_limits<unsigned int>::max() };
+    DeprecatedAndRemoved d_orderingMethod;
+    std::string m_selectedOrderingMethod;
 
-    std::variant<NaturalOrderSolver, AMDOrderSolver, COLAMDOrderSolver, MetisOrderSolver> m_solver;
+    std::unique_ptr<BaseEigenSolverProxy> m_solver;
 
-    Eigen::ComputationInfo getSolverInfo() const;
+    [[nodiscard]] Eigen::ComputationInfo getSolverInfo() const;
     void updateSolverOderingMethod();
 
     sofa::linearalgebra::CompressedRowSparseMatrix<Real> Mfiltered;
@@ -80,8 +82,6 @@ class EigenDirectSparseSolver
     typename sofa::linearalgebra::CompressedRowSparseMatrix<Real>::VecIndex MfilteredrowBegin;
     typename sofa::linearalgebra::CompressedRowSparseMatrix<Real>::VecIndex MfilteredcolsIndex;
 
-    EigenDirectSparseSolver();
-
     static constexpr unsigned int s_defaultOrderingMethod { 1 };
 };
 

Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenDirectSparseSolver.inl

@@ -21,6 +21,7 @@
 ******************************************************************************/
 #pragma once
 #include <sofa/component/linearsolver/direct/EigenDirectSparseSolver.h>
+#include <sofa/core/ComponentLibrary.h>
 
 #include <sofa/helper/ScopedAdvancedTimer.h>
 
@@ -50,10 +51,7 @@ void EigenDirectSparseSolver<TBlockType, EigenSolver>
     EigenVectorXdMap xMap(x.ptr(), x.size());
     EigenVectorXdMap bMap(b.ptr(), b.size());
 
-    std::visit([&bMap, &xMap](auto&& solver)
-    {
-        xMap = solver.solve(bMap);
-    }, m_solver);
+    m_solver->solve(bMap, xMap);
 }
 
 template <class TBlockType, class EigenSolver>
@@ -76,21 +74,15 @@ void EigenDirectSparseSolver<TBlockType, EigenSolver>
     if (analyzePattern)
     {
         SCOPED_TIMER_VARNAME(patternAnalysisTimer, "patternAnalysis");
-        std::visit([this](auto&& solver)
-        {
-            solver.analyzePattern(*m_map);
-        }, m_solver);
+        m_solver->analyzePattern(*m_map);
 
         MfilteredrowBegin = Mfiltered.rowBegin;
         MfilteredcolsIndex = Mfiltered.colsIndex;
     }
 
     {
         SCOPED_TIMER_VARNAME(factorizeTimer, "factorization");
-        std::visit([this](auto&& solver)
-        {
-            solver.factorize(*m_map);
-        }, m_solver);
+        m_solver->factorize(*m_map);
     }
 
     msg_error_when(getSolverInfo() == Eigen::ComputationInfo::InvalidInput) << "Solver cannot factorize: invalid input";
@@ -102,46 +94,50 @@ template <class TBlockType, class EigenSolver>
 Eigen::ComputationInfo EigenDirectSparseSolver<TBlockType, EigenSolver>
 ::getSolverInfo() const
 {
-    Eigen::ComputationInfo info;
-    std::visit([&info](auto&& solver)
-    {
-        info = solver.info();
-    }, m_solver);
-    return info;
+    return m_solver->info();
 }
 
 template <class TBlockType, class EigenSolver>
 void EigenDirectSparseSolver<TBlockType, EigenSolver>::updateSolverOderingMethod()
 {
-    if (m_selectedOrderingMethod != d_orderingMethod.getValue().getSelectedId())
+    if (this->l_orderingMethod)
     {
-        switch(d_orderingMethod.getValue().getSelectedId())
+        if (m_selectedOrderingMethod != this->l_orderingMethod->methodName())
         {
-        case 0:  m_solver.template emplace<std::variant_alternative_t<0, decltype(m_solver)> >(); break;
-        case 1:  m_solver.template emplace<std::variant_alternative_t<1, decltype(m_solver)> >(); break;
-        case 2:  m_solver.template emplace<std::variant_alternative_t<2, decltype(m_solver)> >(); break;
-        case 3:  m_solver.template emplace<std::variant_alternative_t<3, decltype(m_solver)> >(); break;
-        default: m_solver.template emplace<std::variant_alternative_t<s_defaultOrderingMethod, decltype(m_solver)> >(); break;
+            m_selectedOrderingMethod = this->l_orderingMethod->methodName();
+
+            if (EigenSolverFactory::template hasSolver<Real>(m_selectedOrderingMethod))
+            {
+                m_solver = std::unique_ptr<BaseEigenSolverProxy>(EigenSolverFactory::template getSolver<Real>(m_selectedOrderingMethod));
+            }
+            else
+            {
+                std::set<std::string> listAvailableOrderingMethods;
+                for (const auto& [orderingMethodName, _] : EigenSolverFactory::registeredSolvers())
+                {
+                    if (EigenSolverFactory::template hasSolver<Real>(orderingMethodName))
+                    {
+                        listAvailableOrderingMethods.insert(orderingMethodName);
+                    }
+                }
+
+                msg_error() << "This solver does not support the ordering method called '"
+                    << m_selectedOrderingMethod << "' found in the component "
+                    << this->l_orderingMethod->getPathName() << ". The list of available methods are: "
+                    << sofa::helper::join(listAvailableOrderingMethods, ",");
+                this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid);
+            }
+
+            MfilteredrowBegin.clear();
+            MfilteredcolsIndex.clear();
+            m_map.reset();
         }
-        m_selectedOrderingMethod = d_orderingMethod.getValue().getSelectedId();
-        if (m_selectedOrderingMethod >= std::variant_size_v<decltype(m_solver)>)
-            m_selectedOrderingMethod = s_defaultOrderingMethod;
-
-        MfilteredrowBegin.clear();
-        MfilteredcolsIndex.clear();
-        m_map.reset();
     }
-}
-
-template <class TBlockType, class EigenSolver>
-EigenDirectSparseSolver<TBlockType, EigenSolver>::EigenDirectSparseSolver()
-    : Inherit1()
-    , d_orderingMethod(initData(&d_orderingMethod, "ordering", "Ordering method"))
-{
-    sofa::helper::OptionsGroup d_orderingMethodOptions{"Natural", "AMD", "COLAMD", "Metis"};
-
-    d_orderingMethodOptions.setSelectedItem(s_defaultOrderingMethod);
-    d_orderingMethod.setValue(d_orderingMethodOptions);
+    else
+    {
+        msg_fatal() << "OrderingMethod missing.";
+        this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid);
+    }
 }
 
 }

Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSimplicialLDLT.h

@@ -23,7 +23,7 @@
 #include <sofa/component/linearsolver/direct/config.h>
 
 #include <sofa/component/linearsolver/direct/EigenDirectSparseSolver.h>
-#include <sofa/component/linearsolver/direct/SimplicialLDLTTraits.h>
+#include <sofa/component/linearsolver/direct/EigenSolverFactory.h>
 
 namespace sofa::component::linearsolver::direct
 {
@@ -37,15 +37,15 @@ template<class TBlockType>
 class EigenSimplicialLDLT
     : public EigenDirectSparseSolver<
         TBlockType,
-        SimplicialLDLTTraits<typename sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType>::Real>
+        MainSimplicialLDLTFactory
     >
 {
 public:
     typedef sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType> Matrix;
     using Real = typename Matrix::Real;
     typedef sofa::linearalgebra::FullVector<Real> Vector;
 
-    SOFA_CLASS(SOFA_TEMPLATE(EigenSimplicialLDLT, TBlockType), SOFA_TEMPLATE2(EigenDirectSparseSolver, TBlockType, SimplicialLDLTTraits<Real>));
+    SOFA_CLASS(SOFA_TEMPLATE(EigenSimplicialLDLT, TBlockType), SOFA_TEMPLATE2(EigenDirectSparseSolver, TBlockType, MainSimplicialLDLTFactory));
 
 };
 

Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSimplicialLLT.h

@@ -23,7 +23,6 @@
 #include <sofa/component/linearsolver/direct/config.h>
 
 #include <sofa/component/linearsolver/direct/EigenDirectSparseSolver.h>
-#include <sofa/component/linearsolver/direct/SimplicialLLTTraits.h>
 
 namespace sofa::component::linearsolver::direct
 {
@@ -37,15 +36,15 @@ template<class TBlockType>
 class EigenSimplicialLLT
     : public EigenDirectSparseSolver<
         TBlockType,
-        SimplicialLLTTraits<typename sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType>::Real>
+        MainSimplicialLLTFactory
     >
 {
 public:
     typedef sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType> Matrix;
     using Real = typename Matrix::Real;
     typedef sofa::linearalgebra::FullVector<Real> Vector;
 
-    SOFA_CLASS(SOFA_TEMPLATE(EigenSimplicialLLT, TBlockType), SOFA_TEMPLATE2(EigenDirectSparseSolver, TBlockType, SimplicialLLTTraits<Real>));
+    SOFA_CLASS(SOFA_TEMPLATE(EigenSimplicialLLT, TBlockType), SOFA_TEMPLATE2(EigenDirectSparseSolver, TBlockType, MainSimplicialLLTFactory));
 
 };
 

Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/EigenSolverFactory.cpp

@@ -0,0 +1,32 @@
+/******************************************************************************
+*                 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 <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#include <sofa/component/linearsolver/direct/EigenSolverFactory.h>
+
+namespace sofa::component::linearsolver::direct
+{
+
+MainSimplicialLDLTFactory::~MainSimplicialLDLTFactory() = default;
+MainSimplicialLLTFactory::~MainSimplicialLLTFactory() = default;
+MainQRFactory::~MainQRFactory() = default;
+MainLUFactory::~MainLUFactory() = default;
+
+}