Splines LAPACK backend

include/ddc/kernels/splines.hpp

@@ -20,6 +20,7 @@
 #include "splines/spline_evaluator_2d.hpp"
 #include "splines/splines_linear_problem.hpp"
 #include "splines/splines_linear_problem_2x2_blocks.hpp"
+#include "splines/splines_linear_problem_3x3_blocks.hpp"
 #include "splines/splines_linear_problem_band.hpp"
 #include "splines/splines_linear_problem_dense.hpp"
 #include "splines/splines_linear_problem_maker.hpp"

include/ddc/kernels/splines/spline_builder.hpp

@@ -22,7 +22,8 @@ namespace ddc {
  * An enum determining the backend solver of a SplineBuilder or SplineBuilder2d.
  */
 enum class SplineSolver {
-    GINKGO ///< Enum member to identify the Ginkgo-based solver (iterative method)
+    GINKGO, ///< Enum member to identify the Ginkgo-based solver (iterative method)
+    LAPACK ///< Enum member to identify the LAPACK-based solver (direct method)
 };
 
 /**
@@ -603,10 +604,36 @@ void SplineBuilder<
         return;
 	*/
 
-    matrix = ddc::detail::SplinesLinearProblemMaker::make_new_sparse<ExecSpace>(
-            ddc::discrete_space<BSplines>().nbasis(),
-            cols_per_chunk,
-            preconditionner_max_block_size);
+    if constexpr (Solver == ddc::SplineSolver::LAPACK) {
+        int upper_band_width;
+        if (bsplines_type::is_uniform()) {
+            upper_band_width = bsplines_type::degree() / 2;
+        } else {
+            upper_band_width = bsplines_type::degree() - 1;
+        }
+        if constexpr (bsplines_type::is_periodic()) {
+            matrix = ddc::detail::SplinesLinearProblemMaker::make_new_periodic_band_matrix<
+                    ExecSpace>(
+                    ddc::discrete_space<BSplines>().nbasis(),
+                    upper_band_width,
+                    upper_band_width,
+                    bsplines_type::is_uniform());
+        } else {
+            matrix = ddc::detail::SplinesLinearProblemMaker::
+                    make_new_block_matrix_with_band_main_block<ExecSpace>(
+                            ddc::discrete_space<BSplines>().nbasis(),
+                            upper_band_width,
+                            upper_band_width,
+                            bsplines_type::is_uniform(),
+                            lower_block_size,
+                            upper_block_size);
+        }
+    } else if constexpr (Solver == ddc::SplineSolver::GINKGO) {
+        matrix = ddc::detail::SplinesLinearProblemMaker::make_new_sparse<ExecSpace>(
+                ddc::discrete_space<BSplines>().nbasis(),
+                cols_per_chunk,
+                preconditionner_max_block_size);
+    }
 
     build_matrix_system();
 

include/ddc/kernels/splines/splines_linear_problem_2x2_blocks.hpp

@@ -50,7 +50,7 @@ class SplinesLinearProblem2x2Blocks : public SplinesLinearProblem<ExecSpace>
      * @brief SplinesLinearProblem2x2Blocks constructor.
      *
      * @param mat_size The size of one of the dimensions of the square matrix.
-     * @param top_left_block A pointer toward the top-left SplinesLinearProblem. `setup_solver` must not have been called on `q`.
+     * @param top_left_block A pointer toward the top-left SplinesLinearProblem. `setup_solver` must not have been called on it.
      */
     explicit SplinesLinearProblem2x2Blocks(
             std::size_t const mat_size,

include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp

@@ -0,0 +1,185 @@
+// Copyright (C) The DDC development team, see COPYRIGHT.md file
+//
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include <cassert>
+#include <memory>
+#include <string>
+
+#include "splines_linear_problem.hpp"
+#include "splines_linear_problem_2x2_blocks.hpp"
+
+namespace ddc::detail {
+
+/**
+ * @brief A 3x3-blocks linear problem dedicated to the computation of a spline approximation,
+ * with all blocks except center one being stored in dense format.
+ *
+ * A = | a | b | c |
+ *     | d | Q | e |
+ *     | f | g | h |
+ *
+ * The storage format is dense for all blocks except center one, whose storage format is determined by its type.
+ *
+ * The matrix itself and blocks a, Q and h are square (which fully determines the dimensions of the others).
+ *
+ * This class implements row & columns interchanges of the matrix and of multiple right-hand sides to restructure the
+ * 3x3-blocks linear problem into a 2x2-blocks linear problem, relying then on the operations implemented in SplinesLinearProblem2x2Blocks.
+ *
+ * @tparam ExecSpace The Kokkos::ExecutionSpace on which operations related to the matrix are supposed to be performed.
+ */
+template <class ExecSpace>
+class SplinesLinearProblem3x3Blocks : public SplinesLinearProblem2x2Blocks<ExecSpace>
+{
+public:
+    using typename SplinesLinearProblem2x2Blocks<ExecSpace>::MultiRHS;
+    using SplinesLinearProblem2x2Blocks<ExecSpace>::size;
+    using SplinesLinearProblem2x2Blocks<ExecSpace>::solve;
+    using SplinesLinearProblem2x2Blocks<ExecSpace>::m_top_left_block;
+
+protected:
+    std::size_t m_top_size;
+
+public:
+    /**
+     * @brief SplinesLinearProblem3x3Blocks constructor.
+     *
+     * @param mat_size The size of one of the dimensions of the square matrix.
+     * @param top_size The size of one of the dimensions of the top-left square block.
+     * @param center_block A pointer toward the center SplinesLinearProblem. `setup_solver` must not have been called on it.
+     */
+    explicit SplinesLinearProblem3x3Blocks(
+            std::size_t const mat_size,
+            std::size_t const top_size,
+            std::unique_ptr<SplinesLinearProblem<ExecSpace>> center_block)
+        : SplinesLinearProblem2x2Blocks<ExecSpace>(mat_size, std::move(center_block))
+        , m_top_size(top_size)
+    {
+    }
+
+private:
+    /// @brief Adjust indices, governs the row & columns interchanges to restructure the 3x3-blocks matrix into a 2x2-blocks matrix.
+    void adjust_indices(std::size_t& i, std::size_t& j) const
+    {
+        std::size_t const nq = m_top_left_block->size(); // size of the center block
+
+        if (i < m_top_size) {
+            i += nq;
+        } else if (i < m_top_size + nq) {
+            i -= m_top_size;
+        }
+
+        if (j < m_top_size) {
+            j += nq;
+        } else if (j < m_top_size + nq) {
+            j -= m_top_size;
+        }
+    }
+
+public:
+    double get_element(std::size_t i, std::size_t j) const override
+    {
+        adjust_indices(i, j);
+        return SplinesLinearProblem2x2Blocks<ExecSpace>::get_element(i, j);
+    }
+
+    void set_element(std::size_t i, std::size_t j, double const aij) override
+    {
+        adjust_indices(i, j);
+        return SplinesLinearProblem2x2Blocks<ExecSpace>::set_element(i, j, aij);
+    }
+
+private:
+    /**
+     * @brief Perform row interchanges on multiple right-hand sides to get a 2-blocks structure (matching the requirements
+     * of the SplinesLinearProblem2x2Blocks solver).
+     *
+     * |  b_top   |    | b_center |
+     * | b_center | -> |  b_top   | -- Considered as a
+     * | b_bottom |    | b_bottom | -- single bottom block
+     *
+     * @param b The multiple right-hand sides.
+     */
+    void interchange_rows_from_3_to_2_blocks_rhs(MultiRHS const b) const
+    {
+        std::size_t const nq = m_top_left_block->size(); // size of the center block
+
+        // prevent Kokkos::deep_copy(b_top_dst, b_top) to be a deep_copy between overlapping allocations
+        assert(nq >= m_top_size);
+
+        MultiRHS const b_top = Kokkos::
+                subview(b, std::pair<std::size_t, std::size_t> {0, m_top_size}, Kokkos::ALL);
+        MultiRHS const b_center = Kokkos::
+                subview(b,
+                        std::pair<std::size_t, std::size_t> {m_top_size, m_top_size + nq},
+                        Kokkos::ALL);
+
+        MultiRHS const b_center_dst
+                = Kokkos::subview(b, std::pair<std::size_t, std::size_t> {0, nq}, Kokkos::ALL);
+        MultiRHS const b_top_dst = Kokkos::
+                subview(b, std::pair<std::size_t, std::size_t> {nq, nq + m_top_size}, Kokkos::ALL);
+
+        MultiRHS const buffer = Kokkos::create_mirror(ExecSpace(), b_center);
+
+        Kokkos::deep_copy(buffer, b_center);
+        Kokkos::deep_copy(b_top_dst, b_top);
+        Kokkos::deep_copy(b_center_dst, buffer);
+    }
+
+    /**
+     * @brief Perform row interchanges on multiple right-hand sides to restore its 3-blocks structure.
+     *
+     * | b_center |    |  b_top   |
+     * |  b_top   | -> | b_center |
+     * | b_bottom |    | b_bottom |
+     *
+     * @param b The multiple right-hand sides.
+     */
+    void interchange_rows_from_2_to_3_blocks_rhs(MultiRHS const b) const
+    {
+        std::size_t const nq = m_top_left_block->size(); // size of the center block
+
+        // prevent Kokkos::deep_copy(b_top, b_top_src) to be a deep_copy between overlapping allocations
+        assert(nq >= m_top_size);
+
+        MultiRHS const b_center_src
+                = Kokkos::subview(b, std::pair<std::size_t, std::size_t> {0, nq}, Kokkos::ALL);
+        MultiRHS const b_top_src = Kokkos::
+                subview(b, std::pair<std::size_t, std::size_t> {nq, nq + m_top_size}, Kokkos::ALL);
+
+        MultiRHS const b_top = Kokkos::
+                subview(b, std::pair<std::size_t, std::size_t> {0, m_top_size}, Kokkos::ALL);
+        MultiRHS const b_center = Kokkos::
+                subview(b,
+                        std::pair<std::size_t, std::size_t> {m_top_size, m_top_size + nq},
+                        Kokkos::ALL);
+
+        MultiRHS const buffer = Kokkos::create_mirror(ExecSpace(), b_center);
+
+        Kokkos::deep_copy(buffer, b_center_src);
+        Kokkos::deep_copy(b_top, b_top_src);
+        Kokkos::deep_copy(b_center, buffer);
+    }
+
+public:
+    /**
+     * @brief Solve the multiple right-hand sides linear problem Ax=b or its transposed version A^tx=b inplace.
+     *
+     * Perform row interchanges on multiple right-hand sides to obtain a 2x2-blocks linear problem and call the SplinesLinearProblem2x2Blocks solver.
+     *
+     * @param[in, out] b A 2D Kokkos::View storing the multiple right-hand sides of the problem and receiving the corresponding solution.
+     * @param transpose Choose between the direct or transposed version of the linear problem.
+     */
+    void solve(MultiRHS const b, bool const transpose) const override
+    {
+        assert(b.extent(0) == size());
+
+        interchange_rows_from_3_to_2_blocks_rhs(b);
+        SplinesLinearProblem2x2Blocks<ExecSpace>::solve(b, transpose);
+        interchange_rows_from_2_to_3_blocks_rhs(b);
+    }
+};
+
+} // namespace ddc::detail

include/ddc/kernels/splines/splines_linear_problem_maker.hpp

@@ -8,6 +8,7 @@
 #include <optional>
 
 #include "splines_linear_problem_2x2_blocks.hpp"
+#include "splines_linear_problem_3x3_blocks.hpp"
 #include "splines_linear_problem_band.hpp"
 #include "splines_linear_problem_dense.hpp"
 #include "splines_linear_problem_pds_band.hpp"
@@ -65,15 +66,16 @@ class SplinesLinearProblemMaker
     }
 
     /**
-     * @brief Construct a 2x2-blocks linear problem with band "main" block (the one called
-     * Q in SplinesLinearProblem2x2Blocks).
+     * @brief Construct a 2x2-blocks or 3x3-blocks linear problem with band "main" block (the one called
+     * Q in SplinesLinearProblem2x2Blocks and SplinesLinearProblem3x3Blocks).
      *
      * @tparam the Kokkos::ExecutionSpace on which matrix-related operation will be performed.
      * @param n The size of one of the dimensions of the whole square matrix.
      * @param kl The number of subdiagonals in the band block.
      * @param ku The number of superdiagonals in the band block.
      * @param pds A boolean indicating if the band block is positive-definite symetric or not.
-     * @param bottom_right_size The size of one of the dimensions of the bottom-right block.
+     * @param bottom_right_size The size of one of the dimensions of the bottom-right square block.
+     * @param top_left_size The size of one of the dimensions of the top-left square block.
      *
      * @return The SplinesLinearProblem instance.
      */
@@ -84,13 +86,18 @@ class SplinesLinearProblemMaker
             int const kl,
             int const ku,
             bool const pds,
-            int const bottom_size)
+            int const bottom_right_size,
+            int const top_left_size = 0)
     {
-        int const top_size = n - bottom_size;
-        std::unique_ptr<SplinesLinearProblem<ExecSpace>> top_left_block
-                = make_new_band<ExecSpace>(top_size, kl, ku, pds);
+        int const main_size = n - top_left_size - bottom_right_size;
+        std::unique_ptr<SplinesLinearProblem<ExecSpace>> main_block
+                = make_new_band<ExecSpace>(main_size, kl, ku, pds);
+        if (top_left_size == 0) {
+            return std::make_unique<
+                    SplinesLinearProblem2x2Blocks<ExecSpace>>(n, std::move(main_block));
+        }
         return std::make_unique<
-                SplinesLinearProblem2x2Blocks<ExecSpace>>(n, std::move(top_left_block));
+                SplinesLinearProblem3x3Blocks<ExecSpace>>(n, top_left_size, std::move(main_block));
     }
 
     /**

tests/splines/CMakeLists.txt

@@ -101,18 +101,38 @@ foreach(BC "BC_PERIODIC" "BC_GREVILLE")
   endforeach()
 endforeach()
 
+# LAPACK
 foreach(BC "BC_PERIODIC" "BC_GREVILLE" "BC_HERMITE")
   foreach(DEGREE_X RANGE "${SPLINES_TEST_DEGREE_MIN}" "${SPLINES_TEST_DEGREE_MAX}")
     foreach(BSPLINES_TYPE "BSPLINES_TYPE_UNIFORM" "BSPLINES_TYPE_NON_UNIFORM")
-      set(test_name "splines_tests_BATCHED_DEGREE_X_${DEGREE_X}_${BSPLINES_TYPE}_${BC}")
+      set(test_name "splines_tests_BATCHED_DEGREE_X_${DEGREE_X}_${BSPLINES_TYPE}_${BC}_LAPACK")
       add_executable("${test_name}" ../main.cpp batched_spline_builder.cpp)
       target_compile_features("${test_name}" PUBLIC cxx_std_17)
       target_link_libraries("${test_name}"
         PUBLIC
         GTest::gtest
         DDC::DDC
       )
-      target_compile_definitions("${test_name}" PUBLIC -DDEGREE_X=${DEGREE_X} -D${BSPLINES_TYPE} -D${BC})
+    target_compile_definitions("${test_name}" PUBLIC -DDEGREE_X=${DEGREE_X} -D${BSPLINES_TYPE} -D${BC} -DSOLVER_LAPACK)
+      # add_test("${test_name}" "${test_name}")
+      gtest_discover_tests("${test_name}" DISCOVERY_MODE PRE_TEST)
+    endforeach()
+  endforeach()
+endforeach()
+
+# GINKGO
+foreach(BC "BC_PERIODIC")
+  foreach(DEGREE_X RANGE "${SPLINES_TEST_DEGREE_MIN}" "${SPLINES_TEST_DEGREE_MAX}")
+    foreach(BSPLINES_TYPE "BSPLINES_TYPE_UNIFORM")
+      set(test_name "splines_tests_BATCHED_DEGREE_X_${DEGREE_X}_${BSPLINES_TYPE}_${BC}_GINKGO")
+      add_executable("${test_name}" ../main.cpp batched_spline_builder.cpp)
+      target_compile_features("${test_name}" PUBLIC cxx_std_17)
+      target_link_libraries("${test_name}"
+        PUBLIC
+        GTest::gtest
+        DDC::DDC
+      )
+    target_compile_definitions("${test_name}" PUBLIC -DDEGREE_X=${DEGREE_X} -D${BSPLINES_TYPE} -D${BC} -DSOLVER_GINKGO)
       # add_test("${test_name}" "${test_name}")
       gtest_discover_tests("${test_name}" DISCOVERY_MODE PRE_TEST)
     endforeach()

tests/splines/batched_spline_builder.cpp

@@ -231,7 +231,11 @@ static void BatchedSplineTest()
             IDim<I, I>,
             s_bcl,
             s_bcr,
+#if defined(SOLVER_LAPACK)
+            ddc::SplineSolver::LAPACK,
+#elif defined(SOLVER_GINKGO)
             ddc::SplineSolver::GINKGO,
+#endif
             IDim<X, I>...>
             spline_builder(dom_vals);
 
@@ -416,18 +420,20 @@ static void BatchedSplineTest()
                         1.0e-14 * max_norm_int));
 }
 
-#if defined(BC_PERIODIC) && defined(BSPLINES_TYPE_UNIFORM)
-#define SUFFIX(name) name##Periodic##Uniform
-#elif defined(BC_PERIODIC) && defined(BSPLINES_TYPE_NON_UNIFORM)
-#define SUFFIX(name) name##Periodic##NonUniform
-#elif defined(BC_GREVILLE) && defined(BSPLINES_TYPE_UNIFORM)
-#define SUFFIX(name) name##Greville##Uniform
-#elif defined(BC_GREVILLE) && defined(BSPLINES_TYPE_NON_UNIFORM)
-#define SUFFIX(name) name##Greville##NonUniform
-#elif defined(BC_HERMITE) && defined(BSPLINES_TYPE_UNIFORM)
-#define SUFFIX(name) name##Hermite##Uniform
-#elif defined(BC_HERMITE) && defined(BSPLINES_TYPE_NON_UNIFORM)
-#define SUFFIX(name) name##Hermite##NonUniform
+#if defined(BC_PERIODIC) && defined(BSPLINES_TYPE_UNIFORM) && defined(SOLVER_LAPACK)
+#define SUFFIX(name) name##Lapack##Periodic##Uniform
+#elif defined(BC_PERIODIC) && defined(BSPLINES_TYPE_NON_UNIFORM) && defined(SOLVER_LAPACK)
+#define SUFFIX(name) name##Lapack##Periodic##NonUniform
+#elif defined(BC_GREVILLE) && defined(BSPLINES_TYPE_UNIFORM) && defined(SOLVER_LAPACK)
+#define SUFFIX(name) name##Lapack##Greville##Uniform
+#elif defined(BC_GREVILLE) && defined(BSPLINES_TYPE_NON_UNIFORM) && defined(SOLVER_LAPACK)
+#define SUFFIX(name) name##Lapack##Greville##NonUniform
+#elif defined(BC_HERMITE) && defined(BSPLINES_TYPE_UNIFORM) && defined(SOLVER_LAPACK)
+#define SUFFIX(name) name##Lapack##Hermite##Uniform
+#elif defined(BC_HERMITE) && defined(BSPLINES_TYPE_NON_UNIFORM) && defined(SOLVER_LAPACK)
+#define SUFFIX(name) name##Lapack##Hermite##NonUniform
+#elif defined(BC_PERIODIC) && defined(BSPLINES_TYPE_UNIFORM) && defined(SOLVER_GINKGO)
+#define SUFFIX(name) name##Ginkgo##Periodic##Uniform
 #endif
 
 TEST(SUFFIX(BatchedSplineHost), 1DX)