From 703c2caf30f452fb64955d34d54e065fc779036b Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 14 Apr 2022 10:46:29 +0200 Subject: [PATCH] Add all built-in solvers and several ordering methods (#1) * Introduce LLT, LU and QR solvers * Choice of the ordering method * Add metis support * Optimal plugins * Link with Metis * find package metis * File encoding * Another try * Find metis for v21.12 * Fix CMakeLists.txt * unused variable * Remove find_package of metis * Missing template keyword * find_package SofaFramework * Make sure the linear solver module is found * Replace hard-coded value * Condition on finding the metis include file * Fix missing guard * Remove warning directive --- CMakeLists.txt | 42 +++++- scenes/Cylinder_EigenSimplicialLDLT.scn | 17 ++- scenes/FEMBAR_EigenConjugateGradient.scn | 17 ++- scenes/FEMBAR_EigenSimplicialLDLT.scn | 17 ++- scenes/FEMBAR_EigenSimplicialLLT.scn | 33 +++++ scenes/FEMBAR_EigenSparseLU.scn | 33 +++++ .../EigenConjugateGradient.h | 2 +- .../EigenDirectSparseSolver.h | 14 ++ .../EigenDirectSparseSolver[CRS].h | 83 +++++++++++ .../EigenDirectSparseSolver[CRS].inl | 140 ++++++++++++++++++ src/EigenLinearSolvers/EigenSimplicialLDLT.h | 2 +- .../EigenSimplicialLDLT[CRS].cpp | 9 +- .../EigenSimplicialLDLT[CRS].h | 35 +---- .../EigenSimplicialLDLT[CRS].inl | 55 ------- src/EigenLinearSolvers/EigenSimplicialLLT.h | 19 +++ .../EigenSimplicialLLT[CRS].cpp | 19 +++ .../EigenSimplicialLLT[CRS].h | 40 +++++ src/EigenLinearSolvers/EigenSparseLU.h | 19 +++ src/EigenLinearSolvers/EigenSparseLU[CRS].cpp | 19 +++ src/EigenLinearSolvers/EigenSparseLU[CRS].h | 40 +++++ src/EigenLinearSolvers/EigenSparseQR.h | 19 +++ src/EigenLinearSolvers/EigenSparseQR[CRS].cpp | 19 +++ src/EigenLinearSolvers/EigenSparseQR[CRS].h | 40 +++++ src/EigenLinearSolvers/FindMetis.h | 14 ++ src/EigenLinearSolvers/SimplicialLDLTTraits.h | 21 +++ src/EigenLinearSolvers/SimplicialLLTTraits.h | 21 +++ src/EigenLinearSolvers/SparseLUTraits.h | 21 +++ src/EigenLinearSolvers/SparseQRTraits.h | 21 +++ 28 files changed, 710 insertions(+), 121 deletions(-) create mode 100644 scenes/FEMBAR_EigenSimplicialLLT.scn create mode 100644 scenes/FEMBAR_EigenSparseLU.scn create mode 100644 src/EigenLinearSolvers/EigenDirectSparseSolver.h create mode 100644 src/EigenLinearSolvers/EigenDirectSparseSolver[CRS].h create mode 100644 src/EigenLinearSolvers/EigenDirectSparseSolver[CRS].inl delete mode 100644 src/EigenLinearSolvers/EigenSimplicialLDLT[CRS].inl create mode 100644 src/EigenLinearSolvers/EigenSimplicialLLT.h create mode 100644 src/EigenLinearSolvers/EigenSimplicialLLT[CRS].cpp create mode 100644 src/EigenLinearSolvers/EigenSimplicialLLT[CRS].h create mode 100644 src/EigenLinearSolvers/EigenSparseLU.h create mode 100644 src/EigenLinearSolvers/EigenSparseLU[CRS].cpp create mode 100644 src/EigenLinearSolvers/EigenSparseLU[CRS].h create mode 100644 src/EigenLinearSolvers/EigenSparseQR.h create mode 100644 src/EigenLinearSolvers/EigenSparseQR[CRS].cpp create mode 100644 src/EigenLinearSolvers/EigenSparseQR[CRS].h create mode 100644 src/EigenLinearSolvers/FindMetis.h create mode 100644 src/EigenLinearSolvers/SimplicialLDLTTraits.h create mode 100644 src/EigenLinearSolvers/SimplicialLLTTraits.h create mode 100644 src/EigenLinearSolvers/SparseLUTraits.h create mode 100644 src/EigenLinearSolvers/SparseQRTraits.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 5c76092..243c1bd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,29 +1,55 @@ cmake_minimum_required(VERSION 3.12) project(EigenLinearSolvers VERSION 1.0 LANGUAGES CXX) -find_package(Sofa.Component.LinearSolver QUIET) -if (NOT Sofa.Component.LinearSolver_FOUND) - find_package(SofaBaseLinearSolver QUIET) +find_package(SofaFramework REQUIRED) +find_package(Sofa.Component.LinearSolver.Direct QUIET) +if (NOT Sofa.Component.LinearSolver.Direct_FOUND) + find_package(SofaSparseSolver QUIET) +endif() + +if (NOT Sofa.Component.LinearSolver.Direct_FOUND AND NOT SofaSparseSolver_FOUND) + message(FATAL "Cannot find the linear solver module in SOFA") endif() # List all files set(EIGENLINEARSOLVERS_SRC_DIR src/${PROJECT_NAME}) set(HEADER_FILES ${EIGENLINEARSOLVERS_SRC_DIR}/config.h.in + ${EIGENLINEARSOLVERS_SRC_DIR}/FindMetis.h + + ${EIGENLINEARSOLVERS_SRC_DIR}/SimplicialLDLTTraits.h + ${EIGENLINEARSOLVERS_SRC_DIR}/SimplicialLLTTraits.h + ${EIGENLINEARSOLVERS_SRC_DIR}/SparseLUTraits.h + ${EIGENLINEARSOLVERS_SRC_DIR}/SparseQRTraits.h ${EIGENLINEARSOLVERS_SRC_DIR}/EigenConjugateGradient.h ${EIGENLINEARSOLVERS_SRC_DIR}/EigenConjugateGradient[CRS].h ${EIGENLINEARSOLVERS_SRC_DIR}/EigenConjugateGradient[CRS].inl + ${EIGENLINEARSOLVERS_SRC_DIR}/EigenDirectSparseSolver.h + ${EIGENLINEARSOLVERS_SRC_DIR}/EigenDirectSparseSolver[CRS].h + ${EIGENLINEARSOLVERS_SRC_DIR}/EigenDirectSparseSolver[CRS].inl + ${EIGENLINEARSOLVERS_SRC_DIR}/EigenSimplicialLDLT.h ${EIGENLINEARSOLVERS_SRC_DIR}/EigenSimplicialLDLT[CRS].h - ${EIGENLINEARSOLVERS_SRC_DIR}/EigenSimplicialLDLT[CRS].inl + + ${EIGENLINEARSOLVERS_SRC_DIR}/EigenSimplicialLLT.h + ${EIGENLINEARSOLVERS_SRC_DIR}/EigenSimplicialLLT[CRS].h + + ${EIGENLINEARSOLVERS_SRC_DIR}/EigenSparseLU.h + ${EIGENLINEARSOLVERS_SRC_DIR}/EigenSparseLU[CRS].h + + ${EIGENLINEARSOLVERS_SRC_DIR}/EigenSparseQR.h + ${EIGENLINEARSOLVERS_SRC_DIR}/EigenSparseQR[CRS].h ) set(SOURCE_FILES ${EIGENLINEARSOLVERS_SRC_DIR}/init.cpp ${EIGENLINEARSOLVERS_SRC_DIR}/EigenConjugateGradient[CRS].cpp ${EIGENLINEARSOLVERS_SRC_DIR}/EigenSimplicialLDLT[CRS].cpp + ${EIGENLINEARSOLVERS_SRC_DIR}/EigenSimplicialLLT[CRS].cpp + ${EIGENLINEARSOLVERS_SRC_DIR}/EigenSparseLU[CRS].cpp + ${EIGENLINEARSOLVERS_SRC_DIR}/EigenSparseQR[CRS].cpp ) set(README_FILES README.md @@ -33,10 +59,10 @@ set(README_FILES add_library(${PROJECT_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES} ${README_FILES}) # Link the plugin library to its dependency(ies). -if (Sofa.Component.LinearSolver_FOUND) - target_link_libraries(${PROJECT_NAME} Sofa.Component.LinearSolver) -elseif(SofaBaseLinearSolver_FOUND) - target_link_libraries(${PROJECT_NAME} SofaBaseLinearSolver) +if (Sofa.Component.LinearSolver.Direct_FOUND) + target_link_libraries(${PROJECT_NAME} Sofa.Component.LinearSolver.Direct) +elseif(SofaSparseSolver_FOUND) + target_link_libraries(${PROJECT_NAME} SofaSparseSolver) endif() # Create package Config, Version & Target files. diff --git a/scenes/Cylinder_EigenSimplicialLDLT.scn b/scenes/Cylinder_EigenSimplicialLDLT.scn index 16c1e3b..bbbf46f 100644 --- a/scenes/Cylinder_EigenSimplicialLDLT.scn +++ b/scenes/Cylinder_EigenSimplicialLDLT.scn @@ -1,17 +1,18 @@ - - - - - - - + + + + + + + + - + diff --git a/scenes/FEMBAR_EigenConjugateGradient.scn b/scenes/FEMBAR_EigenConjugateGradient.scn index 4cdf5ee..891e068 100644 --- a/scenes/FEMBAR_EigenConjugateGradient.scn +++ b/scenes/FEMBAR_EigenConjugateGradient.scn @@ -1,13 +1,16 @@ - - - + + + + + + + + + - - - @@ -16,7 +19,7 @@ - + diff --git a/scenes/FEMBAR_EigenSimplicialLDLT.scn b/scenes/FEMBAR_EigenSimplicialLDLT.scn index 859856e..d7c4bf3 100644 --- a/scenes/FEMBAR_EigenSimplicialLDLT.scn +++ b/scenes/FEMBAR_EigenSimplicialLDLT.scn @@ -1,13 +1,16 @@ - - - + + + + + + + + + - - - @@ -16,7 +19,7 @@ - + diff --git a/scenes/FEMBAR_EigenSimplicialLLT.scn b/scenes/FEMBAR_EigenSimplicialLLT.scn new file mode 100644 index 0000000..2527acf --- /dev/null +++ b/scenes/FEMBAR_EigenSimplicialLLT.scn @@ -0,0 +1,33 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/scenes/FEMBAR_EigenSparseLU.scn b/scenes/FEMBAR_EigenSparseLU.scn new file mode 100644 index 0000000..72c0960 --- /dev/null +++ b/scenes/FEMBAR_EigenSparseLU.scn @@ -0,0 +1,33 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/EigenLinearSolvers/EigenConjugateGradient.h b/src/EigenLinearSolvers/EigenConjugateGradient.h index de113df..900c6ab 100644 --- a/src/EigenLinearSolvers/EigenConjugateGradient.h +++ b/src/EigenLinearSolvers/EigenConjugateGradient.h @@ -5,7 +5,7 @@ namespace EigenLinearSolvers { /** - * EigenConjugateGradient is a SOFA component allowing to solve a of linear equations. This is an essential component + * EigenConjugateGradient is a SOFA component allowing to solve a set of linear equations. This is an essential component * in a SOFA simulation. * * The class is empty because it is partially specialized in other files. diff --git a/src/EigenLinearSolvers/EigenDirectSparseSolver.h b/src/EigenLinearSolvers/EigenDirectSparseSolver.h new file mode 100644 index 0000000..b001187 --- /dev/null +++ b/src/EigenLinearSolvers/EigenDirectSparseSolver.h @@ -0,0 +1,14 @@ +#pragma once +#include + +namespace EigenLinearSolvers +{ + +template +class EigenDirectSparseSolver +{ +public: + virtual ~EigenDirectSparseSolver() = default; +}; + +} diff --git a/src/EigenLinearSolvers/EigenDirectSparseSolver[CRS].h b/src/EigenLinearSolvers/EigenDirectSparseSolver[CRS].h new file mode 100644 index 0000000..a8eb077 --- /dev/null +++ b/src/EigenLinearSolvers/EigenDirectSparseSolver[CRS].h @@ -0,0 +1,83 @@ +#pragma once +#include + +#include + +#if __has_include() +#include +#else +#include +#endif + +#include +#include + +#include +#include + +namespace EigenLinearSolvers +{ + +/** + * Partial template specialization of EigenDirectSparseSolver for a matrix of type CompressedRowSparseMatrix + */ +template +class EigenDirectSparseSolver< + sofa::linearalgebra::CompressedRowSparseMatrix, + sofa::linearalgebra::FullVector::Real>, + EigenSolver > + : public sofa::component::linearsolver::MatrixLinearSolver< + sofa::linearalgebra::CompressedRowSparseMatrix, + sofa::linearalgebra::FullVector::Real> > +{ +public: + typedef sofa::linearalgebra::CompressedRowSparseMatrix Matrix; + using Real = typename Matrix::Real; + typedef sofa::linearalgebra::FullVector Vector; + + SOFA_ABSTRACT_CLASS(SOFA_TEMPLATE3(EigenDirectSparseSolver, Matrix, Vector, EigenSolver), + SOFA_TEMPLATE2(sofa::component::linearsolver::MatrixLinearSolver, Matrix, Vector)); + + using NaturalOrderSolver = typename EigenSolver::NaturalOrderSolver; + using AMDOrderSolver = typename EigenSolver::AMDOrderSolver; + using COLAMDOrderSolver = typename EigenSolver::COLAMDOrderSolver; +#if EIGENLINEARSOLVERS_HAS_METIS_INCLUDE == 1 + using MetisOrderSolver = typename EigenSolver::MetisOrderSolver; +#endif + + ~EigenDirectSparseSolver() override = default; + + void init() override; + void reinit() override; + + using EigenSparseMatrix = Eigen::SparseMatrix; + using EigenSparseMatrixMap = Eigen::Map; + using EigenVectorXdMap = Eigen::Map >; + + void solve (Matrix& A, Vector& x, Vector& b) override; + void invert(Matrix& A) override; + +protected: + + sofa::core::objectmodel::Data d_orderingMethod; + unsigned int m_selectedOrderingMethod { std::numeric_limits::max() }; + + std::variant m_solver; + + Eigen::ComputationInfo getSolverInfo() const; + void updateSolverOderingMethod(); + + sofa::linearalgebra::CompressedRowSparseMatrix Mfiltered; + std::unique_ptr m_map; + + typename sofa::linearalgebra::CompressedRowSparseMatrix::VecIndex MfilteredrowBegin; + typename sofa::linearalgebra::CompressedRowSparseMatrix::VecIndex MfilteredcolsIndex; + + EigenDirectSparseSolver(); +}; + +} diff --git a/src/EigenLinearSolvers/EigenDirectSparseSolver[CRS].inl b/src/EigenLinearSolvers/EigenDirectSparseSolver[CRS].inl new file mode 100644 index 0000000..29cfd5c --- /dev/null +++ b/src/EigenLinearSolvers/EigenDirectSparseSolver[CRS].inl @@ -0,0 +1,140 @@ +#pragma once +#include + +#include + +namespace EigenLinearSolvers +{ +template +void EigenDirectSparseSolver< + sofa::linearalgebra::CompressedRowSparseMatrix, + sofa::linearalgebra::FullVector::Real>, EigenSolver> + ::init() +{ + updateSolverOderingMethod(); +} + +template +void EigenDirectSparseSolver< + sofa::linearalgebra::CompressedRowSparseMatrix, + sofa::linearalgebra::FullVector::Real>, EigenSolver> + ::reinit() +{ + updateSolverOderingMethod(); +} + +template +void EigenDirectSparseSolver, sofa::linearalgebra::FullVector::Real>, EigenSolver > + ::solve(Matrix& A, Vector& x, Vector& b) +{ + SOFA_UNUSED(A); + + EigenVectorXdMap xMap(x.ptr(), x.size()); + EigenVectorXdMap bMap(b.ptr(), b.size()); + + std::visit([&bMap, &xMap](auto&& solver) + { + xMap = solver.solve(bMap); + }, m_solver); +} + +template +void EigenDirectSparseSolver, sofa::linearalgebra::FullVector::Real>, EigenSolver > + ::invert(Matrix& A) +{ + { + sofa::helper::ScopedAdvancedTimer copyTimer("copyMatrixData"); + Mfiltered.copyNonZeros(A); + Mfiltered.compress(); + } + + if (!m_map) + { + m_map = std::make_unique(Mfiltered.rows(), Mfiltered.cols(), Mfiltered.getColsValue().size(), + (typename EigenSparseMatrixMap::StorageIndex*)Mfiltered.rowBegin.data(), (typename EigenSparseMatrixMap::StorageIndex*)Mfiltered.colsIndex.data(), Mfiltered.colsValue.data()); + } + + const bool analyzePattern = (MfilteredrowBegin != Mfiltered.rowBegin) || (MfilteredcolsIndex != Mfiltered.colsIndex); + + if (analyzePattern) + { + sofa::helper::ScopedAdvancedTimer patternAnalysisTimer("patternAnalysis"); + std::visit([this](auto&& solver) + { + solver.analyzePattern(*m_map); + }, m_solver); + + MfilteredrowBegin = Mfiltered.rowBegin; + MfilteredcolsIndex = Mfiltered.colsIndex; + } + + { + sofa::helper::ScopedAdvancedTimer factorizeTimer("factorization"); + std::visit([this](auto&& solver) + { + solver.factorize(*m_map); + }, m_solver); + } + + msg_error_when(getSolverInfo() == Eigen::ComputationInfo::InvalidInput) << "Solver cannot factorize: invalid input"; + msg_error_when(getSolverInfo() == Eigen::ComputationInfo::NoConvergence) << "Solver cannot factorize: no convergence"; + msg_error_when(getSolverInfo() == Eigen::ComputationInfo::NumericalIssue) << "Solver cannot factorize: numerical issue"; +} + +template +Eigen::ComputationInfo EigenDirectSparseSolver< + sofa::linearalgebra::CompressedRowSparseMatrix, + sofa::linearalgebra::FullVector::Real>, EigenSolver> +::getSolverInfo() const +{ + Eigen::ComputationInfo info; + std::visit([&info](auto&& solver) + { + info = solver.info(); + }, m_solver); + return info; +} + +template +void EigenDirectSparseSolver, sofa::linearalgebra::FullVector +::Real>, EigenSolver>::updateSolverOderingMethod() +{ + if (m_selectedOrderingMethod != d_orderingMethod.getValue().getSelectedId()) + { + switch(d_orderingMethod.getValue().getSelectedId()) + { + case 0: m_solver.template emplace >(); break; + case 1: m_solver.template emplace >(); break; + case 2: m_solver.template emplace >(); break; +#if EIGENLINEARSOLVERS_HAS_METIS_INCLUDE == 1 + case 3: m_solver.template emplace >(); break; +#endif + default: m_solver.template emplace >(); break; + } + m_selectedOrderingMethod = d_orderingMethod.getValue().getSelectedId(); + if (m_selectedOrderingMethod >= std::variant_size_v) + m_selectedOrderingMethod = 1; + + MfilteredrowBegin.clear(); + MfilteredcolsIndex.clear(); + m_map.reset(); + } +} + +template +EigenDirectSparseSolver, sofa::linearalgebra::FullVector< +typename sofa::linearalgebra::CompressedRowSparseMatrix::Real>, EigenSolver>::EigenDirectSparseSolver() + : Inherit1() + , d_orderingMethod(initData(&d_orderingMethod, "ordering", "Ordering method")) +{ +#if EIGENLINEARSOLVERS_HAS_METIS_INCLUDE == 1 + sofa::helper::OptionsGroup d_orderingMethodOptions(4,"Natural", "AMD", "COLAMD", "Metis"); +#else + sofa::helper::OptionsGroup d_orderingMethodOptions(3,"Natural", "AMD", "COLAMD"); +#endif + + d_orderingMethodOptions.setSelectedItem(1); // default None + d_orderingMethod.setValue(d_orderingMethodOptions); +} + +} diff --git a/src/EigenLinearSolvers/EigenSimplicialLDLT.h b/src/EigenLinearSolvers/EigenSimplicialLDLT.h index e0ecb41..9f26a50 100644 --- a/src/EigenLinearSolvers/EigenSimplicialLDLT.h +++ b/src/EigenLinearSolvers/EigenSimplicialLDLT.h @@ -4,7 +4,7 @@ namespace EigenLinearSolvers { /** - * EigenSimplicialLDLT is a SOFA component allowing to solve a of linear equations. This is an essential component + * EigenSimplicialLDLT is a SOFA component allowing to solve a set of linear equations. This is an essential component * in a SOFA simulation. * * The class is empty because it is partially specialized in other files. diff --git a/src/EigenLinearSolvers/EigenSimplicialLDLT[CRS].cpp b/src/EigenLinearSolvers/EigenSimplicialLDLT[CRS].cpp index e00c20e..4afa776 100644 --- a/src/EigenLinearSolvers/EigenSimplicialLDLT[CRS].cpp +++ b/src/EigenLinearSolvers/EigenSimplicialLDLT[CRS].cpp @@ -1,11 +1,6 @@ #define EIGENLINEARSOLVERS_EIGENEIGENSIMPLICIALLDLT_CRS_CPP -#include - -#if __has_include() -#include -#else -#include -#endif +#include +#include #include diff --git a/src/EigenLinearSolvers/EigenSimplicialLDLT[CRS].h b/src/EigenLinearSolvers/EigenSimplicialLDLT[CRS].h index 1188dac..7c14362 100644 --- a/src/EigenLinearSolvers/EigenSimplicialLDLT[CRS].h +++ b/src/EigenLinearSolvers/EigenSimplicialLDLT[CRS].h @@ -2,15 +2,8 @@ #include #include - -#if __has_include() -#include -#else -#include -#endif - -#include -#include +#include +#include namespace EigenLinearSolvers { @@ -22,31 +15,19 @@ template class EigenSimplicialLDLT< sofa::linearalgebra::CompressedRowSparseMatrix, sofa::linearalgebra::FullVector::Real> > - : public sofa::component::linearsolver::MatrixLinearSolver< - sofa::linearalgebra::CompressedRowSparseMatrix, sofa::linearalgebra::FullVector::Real> > + : public EigenDirectSparseSolver< + sofa::linearalgebra::CompressedRowSparseMatrix, + sofa::linearalgebra::FullVector::Real>, + SimplicialLDLTTraits::Real> + > { public: typedef sofa::linearalgebra::CompressedRowSparseMatrix Matrix; using Real = typename Matrix::Real; typedef sofa::linearalgebra::FullVector Vector; - SOFA_CLASS(SOFA_TEMPLATE2(EigenSimplicialLDLT, Matrix, Vector),SOFA_TEMPLATE2(sofa::component::linearsolver::MatrixLinearSolver, Matrix, Vector)); - - using EigenSparseMatrix = Eigen::SparseMatrix; - using EigenSparseMatrixMap = Eigen::Map; - using EigenVectorXdMap = Eigen::Map >; - - void solve (Matrix& A, Vector& x, Vector& b) override; - void invert(Matrix& A) override; - -protected: - Eigen::SimplicialLDLT > m_solver; - - sofa::linearalgebra::CompressedRowSparseMatrix Mfiltered; - std::unique_ptr m_map; + SOFA_CLASS(SOFA_TEMPLATE2(EigenSimplicialLDLT, Matrix, Vector),SOFA_TEMPLATE3(EigenDirectSparseSolver, Matrix, Vector, SimplicialLDLTTraits)); - typename sofa::linearalgebra::CompressedRowSparseMatrix::VecIndex MfilteredrowBegin; - typename sofa::linearalgebra::CompressedRowSparseMatrix::VecIndex MfilteredcolsIndex; }; #ifndef EIGENLINEARSOLVERS_EIGENEIGENSIMPLICIALLDLT_CRS_CPP diff --git a/src/EigenLinearSolvers/EigenSimplicialLDLT[CRS].inl b/src/EigenLinearSolvers/EigenSimplicialLDLT[CRS].inl deleted file mode 100644 index a94a348..0000000 --- a/src/EigenLinearSolvers/EigenSimplicialLDLT[CRS].inl +++ /dev/null @@ -1,55 +0,0 @@ -#pragma once -#include - -#include - -namespace EigenLinearSolvers -{ - -template -void EigenSimplicialLDLT, sofa::linearalgebra::FullVector::Real> > - ::solve(Matrix& A, Vector& x, Vector& b) -{ - EigenVectorXdMap xMap(x.ptr(), x.size()); - EigenVectorXdMap bMap(b.ptr(), b.size()); - xMap = m_solver.solve(bMap); -} - -template -void EigenSimplicialLDLT, sofa::linearalgebra::FullVector::Real> > - ::invert(Matrix& A) -{ - { - sofa::helper::ScopedAdvancedTimer copyTimer("copyMatrixData"); - Mfiltered.copyNonZeros(A); - Mfiltered.compress(); - } - - if (!m_map) - { - m_map = std::make_unique(Mfiltered.rows(), Mfiltered.cols(), Mfiltered.getColsValue().size(), - (typename EigenSparseMatrixMap::StorageIndex*)Mfiltered.rowBegin.data(), (typename EigenSparseMatrixMap::StorageIndex*)Mfiltered.colsIndex.data(), Mfiltered.colsValue.data()); - } - - const bool analyzePattern = (MfilteredrowBegin != Mfiltered.rowBegin) || (MfilteredcolsIndex != Mfiltered.colsIndex); - - if (analyzePattern) - { - sofa::helper::ScopedAdvancedTimer patternAnalysisTimer("patternAnalysis"); - m_solver.analyzePattern(*m_map); - - MfilteredrowBegin = Mfiltered.rowBegin; - MfilteredcolsIndex = Mfiltered.colsIndex; - } - - { - sofa::helper::ScopedAdvancedTimer factorizeTimer("factorization"); - m_solver.factorize(*m_map); - } - - msg_error_when(m_solver.info() == Eigen::ComputationInfo::InvalidInput) << "Solver cannot factorize: invalid input"; - msg_error_when(m_solver.info() == Eigen::ComputationInfo::NoConvergence) << "Solver cannot factorize: no convergence"; - msg_error_when(m_solver.info() == Eigen::ComputationInfo::NumericalIssue) << "Solver cannot factorize: numerical issue"; -} - -} diff --git a/src/EigenLinearSolvers/EigenSimplicialLLT.h b/src/EigenLinearSolvers/EigenSimplicialLLT.h new file mode 100644 index 0000000..82333e7 --- /dev/null +++ b/src/EigenLinearSolvers/EigenSimplicialLLT.h @@ -0,0 +1,19 @@ +#pragma once + +namespace EigenLinearSolvers +{ + +/** + * EigenSimplicialLLT is a SOFA component allowing to solve a set of linear equations. This is an essential component + * in a SOFA simulation. + * + * The class is empty because it is partially specialized in other files. + * + * In all specializations, a linear solver from the Eigen library is used to solve the linear system. + * In this class, the solver is Eigen::SimplicialLLT (https://eigen.tuxfamily.org/dox/classEigen_1_1SimplicialLLT.html) + */ +template +class EigenSimplicialLLT +{}; + +} diff --git a/src/EigenLinearSolvers/EigenSimplicialLLT[CRS].cpp b/src/EigenLinearSolvers/EigenSimplicialLLT[CRS].cpp new file mode 100644 index 0000000..a75e2f0 --- /dev/null +++ b/src/EigenLinearSolvers/EigenSimplicialLLT[CRS].cpp @@ -0,0 +1,19 @@ +#define EIGENLINEARSOLVERS_EIGENEIGENSIMPLICIALLLT_CRS_CPP +#include +#include + +#include + +namespace EigenLinearSolvers +{ + template class SOFA_EIGENLINEARSOLVERS_API EigenSimplicialLLT< sofa::linearalgebra::CompressedRowSparseMatrix, sofa::linearalgebra::FullVector >; + template class SOFA_EIGENLINEARSOLVERS_API EigenSimplicialLLT< + sofa::linearalgebra::CompressedRowSparseMatrix< sofa::type::Mat<3,3,SReal> >, + sofa::linearalgebra::FullVector >; + + int EigenSimplicialLLTCRSClass = sofa::core::RegisterObject("Direct Linear Solver using a Sparse LDL^T factorization.") + .add< EigenSimplicialLLT< sofa::linearalgebra::CompressedRowSparseMatrix, sofa::linearalgebra::FullVector > >() + .add< EigenSimplicialLLT< sofa::linearalgebra::CompressedRowSparseMatrix< sofa::type::Mat<3,3,SReal> >, sofa::linearalgebra::FullVector > >() + ; + +} diff --git a/src/EigenLinearSolvers/EigenSimplicialLLT[CRS].h b/src/EigenLinearSolvers/EigenSimplicialLLT[CRS].h new file mode 100644 index 0000000..ff09619 --- /dev/null +++ b/src/EigenLinearSolvers/EigenSimplicialLLT[CRS].h @@ -0,0 +1,40 @@ +#pragma once +#include + +#include +#include +#include + +namespace EigenLinearSolvers +{ + +/** + * Partial template specialization of EigenSimplicialLLT for a matrix of type CompressedRowSparseMatrix + */ +template +class EigenSimplicialLLT< + sofa::linearalgebra::CompressedRowSparseMatrix, + sofa::linearalgebra::FullVector::Real> > + : public EigenDirectSparseSolver< + sofa::linearalgebra::CompressedRowSparseMatrix, + sofa::linearalgebra::FullVector::Real>, + SimplicialLLTTraits::Real> + > +{ +public: + typedef sofa::linearalgebra::CompressedRowSparseMatrix Matrix; + using Real = typename Matrix::Real; + typedef sofa::linearalgebra::FullVector Vector; + + SOFA_CLASS(SOFA_TEMPLATE2(EigenSimplicialLLT, Matrix, Vector), SOFA_TEMPLATE3(EigenDirectSparseSolver, Matrix, Vector, SimplicialLLTTraits)); + +}; + +#ifndef EIGENLINEARSOLVERS_EIGENEIGENSIMPLICIALLLT_CRS_CPP + extern template class SOFA_EIGENLINEARSOLVERS_API EigenSimplicialLLT< sofa::linearalgebra::CompressedRowSparseMatrix, sofa::linearalgebra::FullVector >; + extern template class SOFA_EIGENLINEARSOLVERS_API EigenSimplicialLLT< + sofa::linearalgebra::CompressedRowSparseMatrix< sofa::type::Mat<3,3,SReal> >, + sofa::linearalgebra::FullVector >; +#endif + +} diff --git a/src/EigenLinearSolvers/EigenSparseLU.h b/src/EigenLinearSolvers/EigenSparseLU.h new file mode 100644 index 0000000..d72cee3 --- /dev/null +++ b/src/EigenLinearSolvers/EigenSparseLU.h @@ -0,0 +1,19 @@ +#pragma once + +namespace EigenLinearSolvers +{ + +/** + * EigenSparseLU is a SOFA component allowing to solve a set of linear equations. This is an essential component + * in a SOFA simulation. + * + * The class is empty because it is partially specialized in other files. + * + * In all specializations, a linear solver from the Eigen library is used to solve the linear system. + * In this class, the solver is Eigen::SparseLU (https://eigen.tuxfamily.org/dox/classEigen_1_1SparseLU.html) + */ +template +class EigenSparseLU +{}; + +} diff --git a/src/EigenLinearSolvers/EigenSparseLU[CRS].cpp b/src/EigenLinearSolvers/EigenSparseLU[CRS].cpp new file mode 100644 index 0000000..d1a0909 --- /dev/null +++ b/src/EigenLinearSolvers/EigenSparseLU[CRS].cpp @@ -0,0 +1,19 @@ +#define EIGENLINEARSOLVERS_EIGENSPARSELU_CRS_CPP +#include +#include + +#include + +namespace EigenLinearSolvers +{ + template class SOFA_EIGENLINEARSOLVERS_API EigenSparseLU< sofa::linearalgebra::CompressedRowSparseMatrix, sofa::linearalgebra::FullVector >; + template class SOFA_EIGENLINEARSOLVERS_API EigenSparseLU< + sofa::linearalgebra::CompressedRowSparseMatrix< sofa::type::Mat<3,3,SReal> >, + sofa::linearalgebra::FullVector >; + + int EigenSparseLUCRSClass = sofa::core::RegisterObject("Direct Linear Solver using a Sparse LDL^T factorization.") + .add< EigenSparseLU< sofa::linearalgebra::CompressedRowSparseMatrix, sofa::linearalgebra::FullVector > >() + .add< EigenSparseLU< sofa::linearalgebra::CompressedRowSparseMatrix< sofa::type::Mat<3,3,SReal> >, sofa::linearalgebra::FullVector > >() + ; + +} diff --git a/src/EigenLinearSolvers/EigenSparseLU[CRS].h b/src/EigenLinearSolvers/EigenSparseLU[CRS].h new file mode 100644 index 0000000..e696bbe --- /dev/null +++ b/src/EigenLinearSolvers/EigenSparseLU[CRS].h @@ -0,0 +1,40 @@ +#pragma once +#include + +#include +#include +#include + +namespace EigenLinearSolvers +{ + +/** + * Partial template specialization of EigenSparseLU for a matrix of type CompressedRowSparseMatrix + */ +template +class EigenSparseLU< + sofa::linearalgebra::CompressedRowSparseMatrix, + sofa::linearalgebra::FullVector::Real> > + : public EigenDirectSparseSolver< + sofa::linearalgebra::CompressedRowSparseMatrix, + sofa::linearalgebra::FullVector::Real>, + SparseLUTraits::Real> + > +{ +public: + typedef sofa::linearalgebra::CompressedRowSparseMatrix Matrix; + using Real = typename Matrix::Real; + typedef sofa::linearalgebra::FullVector Vector; + + SOFA_CLASS(SOFA_TEMPLATE2(EigenSparseLU, Matrix, Vector), SOFA_TEMPLATE3(EigenDirectSparseSolver, Matrix, Vector, SparseLUTraits)); + +}; + +#ifndef EIGENLINEARSOLVERS_EIGENSPARSELU_CRS_CPP + extern template class SOFA_EIGENLINEARSOLVERS_API EigenSparseLU< sofa::linearalgebra::CompressedRowSparseMatrix, sofa::linearalgebra::FullVector >; + extern template class SOFA_EIGENLINEARSOLVERS_API EigenSparseLU< + sofa::linearalgebra::CompressedRowSparseMatrix< sofa::type::Mat<3,3,SReal> >, + sofa::linearalgebra::FullVector >; +#endif + +} diff --git a/src/EigenLinearSolvers/EigenSparseQR.h b/src/EigenLinearSolvers/EigenSparseQR.h new file mode 100644 index 0000000..398ce1a --- /dev/null +++ b/src/EigenLinearSolvers/EigenSparseQR.h @@ -0,0 +1,19 @@ +#pragma once + +namespace EigenLinearSolvers +{ + +/** + * EigenSparseQR is a SOFA component allowing to solve a set of linear equations. This is an essential component + * in a SOFA simulation. + * + * The class is empty because it is partially specialized in other files. + * + * In all specializations, a linear solver from the Eigen library is used to solve the linear system. + * In this class, the solver is Eigen::SparseLU (https://eigen.tuxfamily.org/dox/classEigen_1_1SparseQR.html) + */ +template +class EigenSparseQR +{}; + +} diff --git a/src/EigenLinearSolvers/EigenSparseQR[CRS].cpp b/src/EigenLinearSolvers/EigenSparseQR[CRS].cpp new file mode 100644 index 0000000..1fece7a --- /dev/null +++ b/src/EigenLinearSolvers/EigenSparseQR[CRS].cpp @@ -0,0 +1,19 @@ +#define EIGENLINEARSOLVERS_EIGENSPARSELU_CRS_CPP +#include +#include + +#include + +namespace EigenLinearSolvers +{ + template class SOFA_EIGENLINEARSOLVERS_API EigenSparseQR< sofa::linearalgebra::CompressedRowSparseMatrix, sofa::linearalgebra::FullVector >; + template class SOFA_EIGENLINEARSOLVERS_API EigenSparseQR< + sofa::linearalgebra::CompressedRowSparseMatrix< sofa::type::Mat<3,3,SReal> >, + sofa::linearalgebra::FullVector >; + + int EigenSparseQRCRSClass = sofa::core::RegisterObject("Direct Linear Solver using a Sparse LDL^T factorization.") + .add< EigenSparseQR< sofa::linearalgebra::CompressedRowSparseMatrix, sofa::linearalgebra::FullVector > >() + .add< EigenSparseQR< sofa::linearalgebra::CompressedRowSparseMatrix< sofa::type::Mat<3,3,SReal> >, sofa::linearalgebra::FullVector > >() + ; + +} diff --git a/src/EigenLinearSolvers/EigenSparseQR[CRS].h b/src/EigenLinearSolvers/EigenSparseQR[CRS].h new file mode 100644 index 0000000..2d69ffd --- /dev/null +++ b/src/EigenLinearSolvers/EigenSparseQR[CRS].h @@ -0,0 +1,40 @@ +#pragma once +#include + +#include +#include +#include + +namespace EigenLinearSolvers +{ + +/** + * Partial template specialization of EigenSparseQR for a matrix of type CompressedRowSparseMatrix + */ +template +class EigenSparseQR< + sofa::linearalgebra::CompressedRowSparseMatrix, + sofa::linearalgebra::FullVector::Real> > + : public EigenDirectSparseSolver< + sofa::linearalgebra::CompressedRowSparseMatrix, + sofa::linearalgebra::FullVector::Real>, + SparseQRTraits::Real> + > +{ +public: + typedef sofa::linearalgebra::CompressedRowSparseMatrix Matrix; + using Real = typename Matrix::Real; + typedef sofa::linearalgebra::FullVector Vector; + + SOFA_CLASS(SOFA_TEMPLATE2(EigenSparseQR, Matrix, Vector), SOFA_TEMPLATE3(EigenDirectSparseSolver, Matrix, Vector, SparseQRTraits)); + +}; + +#ifndef EIGENLINEARSOLVERS_EIGENSPARSELU_CRS_CPP + extern template class SOFA_EIGENLINEARSOLVERS_API EigenSparseQR< sofa::linearalgebra::CompressedRowSparseMatrix, sofa::linearalgebra::FullVector >; + extern template class SOFA_EIGENLINEARSOLVERS_API EigenSparseQR< + sofa::linearalgebra::CompressedRowSparseMatrix< sofa::type::Mat<3,3,SReal> >, + sofa::linearalgebra::FullVector >; +#endif + +} diff --git a/src/EigenLinearSolvers/FindMetis.h b/src/EigenLinearSolvers/FindMetis.h new file mode 100644 index 0000000..eeda452 --- /dev/null +++ b/src/EigenLinearSolvers/FindMetis.h @@ -0,0 +1,14 @@ +#pragma once + +/** + * It seems that some binary versions of SOFA do not ship the metis.h include file. + * See https://github.com/sofa-framework/sofa/issues/2866 + * Therefore, metis cannot be supported. + */ +#if __has_include("metis.h") +#define EIGENLINEARSOLVERS_HAS_METIS_INCLUDE 1 +#include +#else +#define EIGENLINEARSOLVERS_HAS_METIS_INCLUDE 0 +#endif + diff --git a/src/EigenLinearSolvers/SimplicialLDLTTraits.h b/src/EigenLinearSolvers/SimplicialLDLTTraits.h new file mode 100644 index 0000000..a92ff4d --- /dev/null +++ b/src/EigenLinearSolvers/SimplicialLDLTTraits.h @@ -0,0 +1,21 @@ +#pragma once +#include + +#include +#include +#include + +namespace EigenLinearSolvers +{ + template + struct SimplicialLDLTTraits + { + using EigenSolver = Eigen::SparseMatrix; + using AMDOrderSolver = Eigen::SimplicialLDLT, Eigen::Lower, Eigen::AMDOrdering >; + using COLAMDOrderSolver = Eigen::SimplicialLDLT, Eigen::Lower, Eigen::COLAMDOrdering >; + using NaturalOrderSolver = Eigen::SimplicialLDLT, Eigen::Lower, Eigen::NaturalOrdering >; +#if EIGENLINEARSOLVERS_HAS_METIS_INCLUDE == 1 + using MetisOrderSolver = Eigen::SimplicialLDLT, Eigen::Lower, Eigen::MetisOrdering >; +#endif + }; +} diff --git a/src/EigenLinearSolvers/SimplicialLLTTraits.h b/src/EigenLinearSolvers/SimplicialLLTTraits.h new file mode 100644 index 0000000..5f3a91c --- /dev/null +++ b/src/EigenLinearSolvers/SimplicialLLTTraits.h @@ -0,0 +1,21 @@ +#pragma once +#include + +#include +#include +#include + +namespace EigenLinearSolvers +{ + template + struct SimplicialLLTTraits + { + using EigenSolver = Eigen::SparseMatrix; + using AMDOrderSolver = Eigen::SimplicialLLT, Eigen::Lower, Eigen::AMDOrdering >; + using COLAMDOrderSolver = Eigen::SimplicialLLT, Eigen::Lower, Eigen::COLAMDOrdering >; + using NaturalOrderSolver = Eigen::SimplicialLLT, Eigen::Lower, Eigen::NaturalOrdering >; +#if EIGENLINEARSOLVERS_HAS_METIS_INCLUDE == 1 + using MetisOrderSolver = Eigen::SimplicialLLT, Eigen::Lower, Eigen::MetisOrdering >; +#endif + }; +} diff --git a/src/EigenLinearSolvers/SparseLUTraits.h b/src/EigenLinearSolvers/SparseLUTraits.h new file mode 100644 index 0000000..54e6dc9 --- /dev/null +++ b/src/EigenLinearSolvers/SparseLUTraits.h @@ -0,0 +1,21 @@ +#pragma once +#include + +#include +#include +#include + +namespace EigenLinearSolvers +{ + template + struct SparseLUTraits + { + using EigenSolver = Eigen::SparseMatrix; + using AMDOrderSolver = Eigen::SparseLU, Eigen::AMDOrdering >; + using COLAMDOrderSolver = Eigen::SparseLU, Eigen::COLAMDOrdering >; + using NaturalOrderSolver = Eigen::SparseLU, Eigen::NaturalOrdering >; +#if EIGENLINEARSOLVERS_HAS_METIS_INCLUDE == 1 + using MetisOrderSolver = Eigen::SparseLU, Eigen::MetisOrdering >; +#endif + }; +} diff --git a/src/EigenLinearSolvers/SparseQRTraits.h b/src/EigenLinearSolvers/SparseQRTraits.h new file mode 100644 index 0000000..6974ad6 --- /dev/null +++ b/src/EigenLinearSolvers/SparseQRTraits.h @@ -0,0 +1,21 @@ +#pragma once +#include + +#include +#include +#include + +namespace EigenLinearSolvers +{ + template + struct SparseQRTraits + { + using EigenSolver = Eigen::SparseMatrix; + using AMDOrderSolver = Eigen::SparseQR, Eigen::AMDOrdering >; + using COLAMDOrderSolver = Eigen::SparseQR, Eigen::COLAMDOrdering >; + using NaturalOrderSolver = Eigen::SparseQR, Eigen::NaturalOrdering >; +#if EIGENLINEARSOLVERS_HAS_METIS_INCLUDE == 1 + using MetisOrderSolver = Eigen::SparseQR, Eigen::MetisOrdering >; +#endif + }; +}