Skip to content

Commit

Permalink
[LinearSolver.Direct] Introduce other ordering methods in SparseLDL t…
Browse files Browse the repository at this point in the history 
…hrough Eigen (#4370)

* [LinearSolver.Direct] Introduce other ordering methods in SparseLDL through Eigen

This gives more choice for the ordering method. It mimics the choices we have in Eigen-based linear system (e.g. EigenSimplicialLDLT).

AMD cannot be selected because it does not compile.

The data applyPermutation is deprecated

* deprecated unused functions

* Move definition of the deprecated data into the header file

Co-authored-by: Hugo <[email protected]>

* Fix compilation with AMD

* Fix unit test by going back to the initial calls to Metis

---------

Co-authored-by: Hugo <[email protected]>
  • Loading branch information
@alxbilger @hugtalbot
alxbilger and hugtalbot authored Jan 24, 2024
1 parent 308cfa5 commit 7e45178
Show file tree
Hide file tree
Showing 4 changed files with 84 additions and 27 deletions.
2 changes: 2 additions & 0 deletions Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseCommon.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,10 +37,12 @@ xadj[i+1]-xadj[i] is the number of neighbors of the i-th node
adj[xadj[i]] is the first neighbor of the i-th node
**/
SOFA_ATTRIBUTE_DEPRECATED__SPARSECOMMON()
SOFA_COMPONENT_LINEARSOLVER_DIRECT_API
void csrToAdj(int n, int * M_colptr, int * M_rowind, type::vector<int>& adj, type::vector<int>& xadj, type::vector<int>& t_adj, type::vector<int>& t_xadj, type::vector<int>& tran_countvec );

// compute the fill reducing permutation via METIS
SOFA_ATTRIBUTE_DEPRECATED__SPARSECOMMON()
SOFA_COMPONENT_LINEARSOLVER_DIRECT_API
void fillReducingPermutation(int nbColumns, int *columns, int* rowIndices,
int * perm,int * invperm);
Expand Down
100 changes: 74 additions & 26 deletions ...omponent/LinearSolver/Direct/src/sofa/component/linearsolver/direct/SparseLDLSolverImpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,12 +28,11 @@
#include <sofa/helper/OptionsGroup.h>
#include <sofa/linearalgebra/DiagonalSystemSolver.h>
#include <sofa/linearalgebra/TriangularSystemSolver.h>
#include <Eigen/SparseCore>
#include <Eigen/MetisSupport>
#include <Eigen/OrderingMethods>


extern "C" {
#include <metis.h>
}

namespace sofa::component::linearsolver::direct
{

Expand Down Expand Up @@ -166,14 +165,19 @@ public :
protected :

Data<bool> d_precomputeSymbolicDecomposition; ///< If true the solver will reuse the precomputed symbolic decomposition. Otherwise it will recompute it at each step.
Data<bool> d_applyPermutation; ///< If true the solver will apply a fill-reducing permutation to the matrix of the system.
core::objectmodel::lifecycle::DeprecatedData d_applyPermutation{this, "v24.06", "v24.12", "applyPermutation", "Use the Data 'ordering'"};
Data<sofa::helper::OptionsGroup> d_orderingMethod; ///< Ordering method
Data<int> d_L_nnz; ///< Number of non-zero values in the lower triangular matrix of the factorization. The lower, the faster the system is solved.

static constexpr const char* s_defaultOrderingMethod { "Metis" };

SparseLDLSolverImpl() : Inherit()
, d_precomputeSymbolicDecomposition(initData(&d_precomputeSymbolicDecomposition, true ,"precomputeSymbolicDecomposition", "If true the solver will reuse the precomputed symbolic decomposition. Otherwise it will recompute it at each step."))
, d_applyPermutation(initData(&d_applyPermutation, true ,"applyPermutation", "If true the solver will apply a fill-reducing permutation to the matrix of the system."))
, d_orderingMethod(initData(&d_orderingMethod, sofa::helper::OptionsGroup{"Natural", "AMD", "COLAMD", "Metis"}, "ordering", "Ordering method"))
, d_L_nnz(initData(&d_L_nnz, 0, "L_nnz", "Number of non-zero values in the lower triangular matrix of the factorization. The lower, the faster the system is solved.", true, true))
{}
{
sofa::helper::getWriteAccessor(d_orderingMethod)->setSelectedItem(s_defaultOrderingMethod);
}

template<class VecInt,class VecReal>
void solve_cpu(Real * x,const Real * b,SparseLDLImplInvertData<VecInt,VecReal> * data)
Expand Down Expand Up @@ -227,29 +231,73 @@ protected :
}
}

void LDL_ordering(int n,int * M_colptr,int * M_rowind,int * perm,int * invperm)
void naturalOrder(int n, int* perm, int* invperm)
{
if( d_applyPermutation.getValue() )
for (int j = 0; j < n; ++j)
{
// METIS
csrToAdj( n, M_colptr, M_rowind, adj, xadj, t_adj, t_xadj, tran_countvec );

/*
int numflag = 0, options = 0;
The new API of metis requires pointers on numflag and "options" which are "structure" to parametrize the factorization
We give NULL and NULL to use the default option (see doc of metis for details) !
If you have the error "SparseLDLSolver failure to factorize, D(k,k) is zero" that probably means that you use the previsou version of metis.
In this case you have to download and install the last version from : www.cs.umn.edu/~metis‎
*/
METIS_NodeND(&n , xadj.data(), adj.data(), nullptr, nullptr, perm,invperm);
perm[j] = j;
invperm[j] = j;
}
else
}

void LDL_ordering(int n, int nnz, int * M_colptr,int * M_rowind, Real* M_values, int * perm,int * invperm)
{
const auto orderingMethod = sofa::helper::getReadAccessor(d_orderingMethod)->getSelectedItem();
if (orderingMethod != "Natural")
{
for(int j=0; j<n ;++j)
const auto computeOrdering = [&](auto ordering)
{
using PermutationType = typename decltype(ordering)::PermutationType;
PermutationType permutation;
using EigenSparseMatrix = Eigen::SparseMatrix<Real, Eigen::ColMajor>;
using EigenSparseMatrixMap = Eigen::Map<const EigenSparseMatrix>;
const auto map = EigenSparseMatrixMap( n, n, nnz, M_colptr, M_rowind, M_values);
ordering.template operator()<const EigenSparseMatrix>(map, permutation);

const PermutationType inv = permutation.inverse();

for (int j = 0; j < n; ++j)
{
perm[j] = j ;
invperm[j] = j ;
perm[j] = permutation.indices()(j);
invperm[j] = inv.indices()(j) ;
}
};

if (orderingMethod == "Metis")
{
// This could be the way to call Metis using the Eigen API, but
// it triggers failing unit tests compared to the other method
// computeOrdering(Eigen::MetisOrdering<int>());

type::vector<int> xadj,adj,t_xadj,t_adj;
csrToAdj( n, M_colptr, M_rowind, adj, xadj, t_adj, t_xadj, tran_countvec );

/*
int numflag = 0, options = 0;
The new API of metis requires pointers on numflag and "options" which are "structure" to parametrize the factorization
We give NULL and NULL to use the default option (see doc of metis for details) !
If you have the error "SparseLDLSolver failure to factorize, D(k,k) is zero" that probably means that you use the previsou version of metis.
In this case you have to download and install the last version from : www.cs.umn.edu/~metis‎
*/
METIS_NodeND(&n , xadj.data(), adj.data(), nullptr, nullptr, perm,invperm);
}
else if (orderingMethod == "COLAMD")
{
computeOrdering(Eigen::COLAMDOrdering<int>());
}
else if (orderingMethod == "AMD")
{
computeOrdering(Eigen::AMDOrdering<int>());
}
else
{
msg_error() << "Ordering method '" << orderingMethod << "' not supported: fallback to natural order";
naturalOrder(n, perm, invperm);
}
}
else
{
naturalOrder(n, perm, invperm);
}
}

Expand Down Expand Up @@ -307,7 +355,7 @@ protected :
memcpy(data->P_rowind.data(),M_rowind,data->P_nnz * sizeof(int));

//ordering function
LDL_ordering( data->n , M_colptr , M_rowind , data->perm.data(), data->invperm.data() );
LDL_ordering( data->n , data->P_nnz, M_colptr , M_rowind , M_values, data->perm.data(), data->invperm.data() );

data->Parent.clear();
data->Parent.resize(data->n);
Expand Down Expand Up @@ -380,7 +428,7 @@ protected :

type::vector<Real> Tmp;
protected : //the following variables are used during the factorization they cannot be used in the main thread !
type::vector<int> xadj,adj,t_xadj,t_adj;

type::vector<Real> Y;
type::vector<int> Lnz,Flag,Pattern;
type::vector<int> tran_countvec;
Expand Down
7 changes: 7 additions & 0 deletions Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/config.h.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,3 +65,10 @@ namespace sofa::component::linearsolver::direct
#define SOFA_ATTRIBUTE_DEPRECATED__SOLVER_DIRECT_VERBOSEDATA() \
SOFA_ATTRIBUTE_DEPRECATED("v23.12", "v24.06", "This Data is no longer used")
#endif // SOFA_BUILD_SOFA_COMPONENT_LINEARSOLVER_DIRECT

#ifdef SOFA_BUILD_SOFA_COMPONENT_LINEARSOLVER_DIRECT
#define SOFA_ATTRIBUTE_DEPRECATED__SPARSECOMMON()
#else
#define SOFA_ATTRIBUTE_DEPRECATED__SPARSECOMMON() \
SOFA_ATTRIBUTE_DEPRECATED("v24.06", "v24.12", "This function is no longer used")
#endif // SOFA_BUILD_SOFA_COMPONENT_LINEARSOLVER_DIRECT
2 changes: 1 addition & 1 deletion examples/Component/LinearSolver/Direct/FEMBAR_SparseLDLSolver.scn
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

<include href="../FEMBAR-common.xml"/>

<SparseLDLSolver applyPermutation="false" template="CompressedRowSparseMatrixMat3x3"/>
<SparseLDLSolver ordering="Natural" template="CompressedRowSparseMatrixMat3x3"/>
<HexahedronFEMForceField name="FEM" youngModulus="4000" poissonRatio="0.3" method="large" />

</Node>

0 comments on commit 7e45178

Please sign in to comment.