SplinesLinearSolver new API

benchmarks/splines.cpp

@@ -185,15 +185,15 @@ static void characteristics_advection(benchmark::State& state)
 
 #ifdef KOKKOS_ENABLE_CUDA
 std::string chip = "gpu";
-int cols_per_chunk_ref = 65535;
+std::size_t cols_per_chunk_ref = 65535;
 unsigned int preconditionner_max_block_size_ref = 1u;
 #elif defined(KOKKOS_ENABLE_OPENMP)
 std::string chip = "cpu";
-int cols_per_chunk_ref = 8192;
+std::size_t cols_per_chunk_ref = 8192;
 unsigned int preconditionner_max_block_size_ref = 32u;
 #elif defined(KOKKOS_ENABLE_SERIAL)
 std::string chip = "cpu";
-int cols_per_chunk_ref = 8192;
+std::size_t cols_per_chunk_ref = 8192;
 unsigned int preconditionner_max_block_size_ref = 32u;
 #endif
 

include/ddc/kernels/splines.hpp

@@ -11,14 +11,14 @@
 #include "splines/greville_interpolation_points.hpp"
 #include "splines/knots_as_interpolation_points.hpp"
 #include "splines/math_tools.hpp"
-#include "splines/matrix.hpp"
-#include "splines/matrix_maker.hpp"
-#include "splines/matrix_sparse.hpp"
 #include "splines/null_extrapolation_rule.hpp"
 #include "splines/periodic_extrapolation_rule.hpp"
 #include "splines/spline_boundary_conditions.hpp"
 #include "splines/spline_builder.hpp"
 #include "splines/spline_builder_2d.hpp"
 #include "splines/spline_evaluator.hpp"
 #include "splines/spline_evaluator_2d.hpp"
+#include "splines/splines_linear_solver.hpp"
+#include "splines/splines_linear_solver_maker.hpp"
+#include "splines/splines_linear_solver_sparse.hpp"
 #include "splines/view.hpp"

include/ddc/kernels/splines/spline_builder.hpp

@@ -8,6 +8,7 @@
 #include "ddc/kokkos_allocator.hpp"
 
 #include "deriv.hpp"
+#include "splines_linear_solver_maker.hpp"
 
 namespace ddc {
 enum class SplineSolver {
@@ -136,15 +137,15 @@ class SplineBuilder
     double m_dx; // average cell size for normalization of derivatives
 
     // interpolator specific
-    std::unique_ptr<ddc::detail::Matrix> matrix;
+    std::unique_ptr<ddc::detail::SplinesLinearProblem<exec_space>> matrix;
 
     /// Calculate offset so that the matrix is diagonally dominant
     int compute_offset(interpolation_domain_type const& interpolation_domain);
 
 public:
     explicit SplineBuilder(
             batched_interpolation_domain_type const& batched_interpolation_domain,
-            std::optional<int> cols_per_chunk = std::nullopt,
+            std::optional<std::size_t> cols_per_chunk = std::nullopt,
             std::optional<unsigned int> preconditionner_max_block_size = std::nullopt)
         : m_batched_interpolation_domain(batched_interpolation_domain)
         , m_offset(compute_offset(interpolation_domain()))
@@ -264,7 +265,7 @@ class SplineBuilder
      *
      * @return A reference to the interpolation matrix.
      */
-    const ddc::detail::Matrix& get_interpolation_matrix() const noexcept
+    const ddc::detail::SplinesLinearProblem<exec_space>& get_interpolation_matrix() const noexcept
     {
         return *matrix;
     }
@@ -311,7 +312,7 @@ class SplineBuilder
     void allocate_matrix(
             int lower_block_size,
             int upper_block_size,
-            std::optional<int> cols_per_chunk = std::nullopt,
+            std::optional<std::size_t> cols_per_chunk = std::nullopt,
             std::optional<unsigned int> preconditionner_max_block_size = std::nullopt);
 
     void build_matrix_system();
@@ -476,7 +477,7 @@ void SplineBuilder<
         allocate_matrix(
                 [[maybe_unused]] int lower_block_size,
                 [[maybe_unused]] int upper_block_size,
-                std::optional<int> cols_per_chunk,
+                std::optional<std::size_t> cols_per_chunk,
                 std::optional<unsigned int> preconditionner_max_block_size)
 {
     // Special case: linear spline
@@ -487,14 +488,14 @@ void SplineBuilder<
         return;
 	*/
 
-    matrix = ddc::detail::MatrixMaker::make_new_sparse<ExecSpace>(
+    matrix = ddc::detail::SplinesLinearProblemMaker::make_new_sparse<ExecSpace>(
             ddc::discrete_space<BSplines>().nbasis(),
             cols_per_chunk,
             preconditionner_max_block_size);
 
     build_matrix_system();
 
-    matrix->factorize();
+    matrix->setup_solver();
 }
 
 template <
@@ -724,7 +725,7 @@ operator()(
             ddc::discrete_space<bsplines_type>().nbasis(),
             batch_domain().size());
     // Compute spline coef
-    matrix->solve_batch_inplace(bcoef_section);
+    matrix->solve(bcoef_section);
     // Transpose back spline_tr in spline
     ddc::parallel_for_each(
             exec_space(),

include/ddc/kernels/splines/spline_builder_2d.hpp

@@ -155,7 +155,7 @@ class SplineBuilder2D
      */
     explicit SplineBuilder2D(
             batched_interpolation_domain_type const& batched_interpolation_domain,
-            std::optional<int> cols_per_chunk = std::nullopt,
+            std::optional<std::size_t> cols_per_chunk = std::nullopt,
             std::optional<unsigned int> preconditionner_max_block_size = std::nullopt)
         : m_spline_builder1(
                 batched_interpolation_domain,