WIP Polymorphic sparse scheme to solve following use case (#46)

[tkloczko]

Let us consider a Finite-Element or Finite-Volume code for which we have to solve a linear system using different libraries (Hypre, Petsc, Trilinos, Mumps, PastiX, ...). We would like to avoid conversions between different formats of sparse matrix. As the linear solver drives the choice of the sparse matrix format, we would like that the matrix the FE code has to fill has directly the format the linear solver requires. The choice of the solver is done at runtime, so it means that the underlying sparse scheme of the matrix should be instanciated at runtime. It implies that the sparse array or sparse tensor is a wrapper around a polymorphic sparse scheme. This latter is a bridge over an abstract scheme that uses type erasure to handle the different format implementations.

The following piece of code tries to illustrate which problem we try to solve. Let us say that we have wrapped several linear solvers into a minimal common API, namely MyLinearSolver. According to the choice of the linear solver, we want the Finite Volume solver to fill the sparse matrix A in the format that the linear solver requires. This choice is done after, for instance reading options in a file. It is thus at runtime that everything should work.

namespace xt {
template <T>
using polymorphic_sparse_array = xsparse_array<T, xsparse_polymorphic_scheme<T>>;
}

MyFVSolver fv_solver;
MyLinearSolver lin_solver;

// Init options and policy
fv_solver.init(options);
lin_solver.init(options);

// Init sparse matrix using the linear solver choice
xt::scheme_policy.setPolicy(lin_solver.format());
xt::polymorphic_sparse_array<double> A; 
// A uses now the scheme defined by the policy but the type remains the same
// enabling to use this unique type as argument for both MyFVSolver and MyLinearSolver API.

//  In the time loop

fv_solver.supplySparseMatrix(A); // No template except the value type of the matrix.
fv_solver.computeConvectiveTerm(); // fill the sparse matrix

lin_solver.setSparseMatrix(A); // No template except the value type of the matrix.
lin_solver.solve();

/// ...

[tkloczko]

In a recent past, we encountered the same difficulties, naming dealing with runtime options and template collections, using CGAL. A similar solution was found to address them.