Splines linear problem 3x3 blocks

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/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/matrix.cpp

@@ -229,7 +229,7 @@ TEST(Matrix, 2x2Blocks)
     std::size_t const k = 10;
     std::unique_ptr<ddc::detail::SplinesLinearProblem<Kokkos::DefaultExecutionSpace>> top_left_block
             = std::make_unique<
-                    ddc::detail::SplinesLinearProblemDense<Kokkos::DefaultExecutionSpace>>(3);
+                    ddc::detail::SplinesLinearProblemDense<Kokkos::DefaultExecutionSpace>>(7);
     std::unique_ptr<ddc::detail::SplinesLinearProblem<Kokkos::DefaultExecutionSpace>> matrix
             = std::make_unique<ddc::detail::SplinesLinearProblem2x2Blocks<
                     Kokkos::DefaultExecutionSpace>>(N, std::move(top_left_block));
@@ -248,6 +248,31 @@ TEST(Matrix, 2x2Blocks)
     solve_and_validate(*matrix);
 }
 
+TEST(Matrix, 3x3Blocks)
+{
+    std::size_t const N = 10;
+    std::size_t const k = 10;
+    std::unique_ptr<ddc::detail::SplinesLinearProblem<Kokkos::DefaultExecutionSpace>> center_block
+            = std::make_unique<
+                    ddc::detail::SplinesLinearProblemDense<Kokkos::DefaultExecutionSpace>>(N - 5);
+    std::unique_ptr<ddc::detail::SplinesLinearProblem<Kokkos::DefaultExecutionSpace>> matrix
+            = std::make_unique<ddc::detail::SplinesLinearProblem3x3Blocks<
+                    Kokkos::DefaultExecutionSpace>>(N, 2, std::move(center_block));
+
+    // Build a non-symmetric full-rank matrix (without zero)
+    for (std::size_t i(0); i < N; ++i) {
+        std::cout << i;
+        matrix->set_element(i, i, 3. / 4 * ((N + 1) * i + 1));
+        for (std::size_t j(std::max(0, int(i) - int(k))); j < i; ++j) {
+            matrix->set_element(i, j, -(1. / 4) / k * (N * i + j + 1));
+        }
+        for (std::size_t j(i + 1); j < std::min(N, i + k + 1); ++j) {
+            matrix->set_element(i, j, -(1. / 4) / k * (N * i + j + 1));
+        }
+    }
+
+    solve_and_validate(*matrix);
+}
 
 class MatrixSizesFixture : public testing::TestWithParam<std::tuple<std::size_t, std::size_t>>
 {
@@ -339,6 +364,27 @@ TEST_P(MatrixSizesFixture, 2x2Blocks)
     solve_and_validate(*matrix);
 }
 
+TEST_P(MatrixSizesFixture, 3x3Blocks)
+{
+    auto const [N, k] = GetParam();
+    std::unique_ptr<ddc::detail::SplinesLinearProblem<Kokkos::DefaultExecutionSpace>> matrix
+            = ddc::detail::SplinesLinearProblemMaker::make_new_block_matrix_with_band_main_block<
+                    Kokkos::DefaultExecutionSpace>(N, k, k, false, 3, 2);
+
+    // Build a non-symmetric full-rank band matrix
+    for (std::size_t i(0); i < N; ++i) {
+        matrix->set_element(i, i, 3. / 4 * ((N + 1) * i + 1));
+        for (std::size_t j(std::max(0, int(i) - int(k))); j < i; ++j) {
+            matrix->set_element(i, j, -(1. / 4) / k * (N * i + j + 1));
+        }
+        for (std::size_t j(i + 1); j < std::min(N, i + k + 1); ++j) {
+            matrix->set_element(i, j, -(1. / 4) / k * (N * i + j + 1));
+        }
+    }
+
+    solve_and_validate(*matrix);
+}
+
 TEST_P(MatrixSizesFixture, PeriodicBand)
 {
     auto const [N, k] = GetParam();