diff --git a/include/ddc/kernels/splines.hpp b/include/ddc/kernels/splines.hpp index 67081eff4..27595fecd 100644 --- a/include/ddc/kernels/splines.hpp +++ b/include/ddc/kernels/splines.hpp @@ -26,6 +26,5 @@ #include "splines/splines_linear_problem_maker.hpp" #include "splines/splines_linear_problem_pds_band.hpp" #include "splines/splines_linear_problem_pds_tridiag.hpp" -#include "splines/splines_linear_problem_periodic_band.hpp" #include "splines/splines_linear_problem_sparse.hpp" #include "splines/view.hpp" diff --git a/include/ddc/kernels/splines/splines_linear_problem_maker.hpp b/include/ddc/kernels/splines/splines_linear_problem_maker.hpp index 9f0524c61..8fe8183d5 100644 --- a/include/ddc/kernels/splines/splines_linear_problem_maker.hpp +++ b/include/ddc/kernels/splines/splines_linear_problem_maker.hpp @@ -13,7 +13,6 @@ #include "splines_linear_problem_dense.hpp" #include "splines_linear_problem_pds_band.hpp" #include "splines_linear_problem_pds_tridiag.hpp" -#include "splines_linear_problem_periodic_band.hpp" #include "splines_linear_problem_sparse.hpp" namespace ddc::detail { @@ -136,44 +135,6 @@ class SplinesLinearProblemMaker return make_new_block_matrix_with_band_main_block(n, kl, ku, pds, bottom_size); } - /** - * @brief Construct a 2x2-blocks linear problem with band "main" block (the one called - * Q in SplinesLinearProblem2x2Blocks). - * - * @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. - * - * @return The SplinesLinearProblem instance. - */ - template - static std::unique_ptr> make_new_periodic_band_matrix( - int const n, - int const kl, - int const ku, - bool const pds) - { - int const bottom_size = std::max(kl, ku); - int const top_size = n - bottom_size; - std::unique_ptr> top_left_block; - if (bottom_size * n + bottom_size * (bottom_size + 1) + (2 * kl + 1 + ku) * top_size - >= n * n) { - return std::make_unique>(n); - } else if (pds && kl == ku && kl == 1) { - top_left_block = std::make_unique>(top_size); - } else if (kl == ku && pds) { - top_left_block = std::make_unique>(top_size, kl); - } else { - top_left_block - = std::make_unique>(top_size, kl, ku); - } - return std::make_unique< - SplinesLinearProblemPeriodicBand>(n, kl, ku, std::move(top_left_block)); - } - /** * @brief Construct a sparse matrix * diff --git a/include/ddc/kernels/splines/splines_linear_problem_periodic_band.hpp b/include/ddc/kernels/splines/splines_linear_problem_periodic_band.hpp deleted file mode 100644 index b92418813..000000000 --- a/include/ddc/kernels/splines/splines_linear_problem_periodic_band.hpp +++ /dev/null @@ -1,270 +0,0 @@ -// Copyright (C) The DDC development team, see COPYRIGHT.md file -// -// SPDX-License-Identifier: MIT - -#pragma once - -#include -#include -#include - -#include - -#include "splines_linear_problem.hpp" -#include "splines_linear_problem_dense.hpp" - -namespace ddc::detail { - -/** - * @brief A periodic-band linear problem dedicated to the computation of a spline approximation (taking in account boundary conditions), - * with all blocks except top-left one being stored in dense format. - * - * A = | Q | gamma | - * | lambda | delta | - * - * The storage format is dense row-major for top-left, top-right and bottom-left blocks, and determined by - * its type for the top-left block. - * - * This class implements a Schur complement method to perform a block-LU factorization and solve, - * calling top-left block and bottom-right block setup_solver() and solve() methods for internal operations. - * - * @tparam ExecSpace The Kokkos::ExecutionSpace on which operations related to the matrix are supposed to be performed. - */ -template -class SplinesLinearProblemPeriodicBand : public SplinesLinearProblem2x2Blocks -{ -public: - using typename SplinesLinearProblem2x2Blocks::MultiRHS; - using SplinesLinearProblem2x2Blocks::size; - -protected: - std::size_t m_kl; // no. of subdiagonals - std::size_t m_ku; // no. of superdiagonals - using SplinesLinearProblem2x2Blocks::m_top_left_block; - using SplinesLinearProblem2x2Blocks::m_top_right_block; - using SplinesLinearProblem2x2Blocks::m_bottom_left_block; - using SplinesLinearProblem2x2Blocks::m_bottom_right_block; - using SplinesLinearProblem2x2Blocks::gemv_minus1_1; - -public: - /** - * @brief SplinesLinearProblem2x2Blocks constructor. - * - * @param mat_size The size of one of the dimensions of the square matrix. - * @param q A pointer toward the top-left SplinesLinearProblem. - */ - explicit SplinesLinearProblemPeriodicBand( - std::size_t const mat_size, - std::size_t const kl, - std::size_t const ku, - std::unique_ptr> top_left_block) - : SplinesLinearProblem2x2Blocks( - mat_size, - std::move(top_left_block), - std::max(kl, ku), - std::max(kl, ku) + 1) - , m_kl(kl) - , m_ku(ku) - { - } - - double get_element(std::size_t const i, std::size_t const j) const override - { - assert(i < size()); - assert(j < size()); - - std::size_t const nq = m_top_left_block->size(); - std::size_t const ndelta = m_bottom_right_block->size(); - if (i >= nq && j < nq) { - std::ptrdiff_t d = j - i; - if (d > (std::ptrdiff_t)(size() / 2)) - d -= size(); - if (d < -(std::ptrdiff_t)(size() / 2)) - d += size(); - - if (d < -(std::ptrdiff_t)m_kl || d > (std::ptrdiff_t)m_ku) - return 0.0; - if (d > (std::ptrdiff_t)0) { - return m_bottom_left_block.h_view(i - nq, j); - } else { - return m_bottom_left_block.h_view(i - nq, j - nq + ndelta + 1); - } - } else { - return SplinesLinearProblem2x2Blocks::get_element(i, j); - } - } - - void set_element(std::size_t const i, std::size_t const j, double const aij) override - { - assert(i < size()); - assert(j < size()); - - std::size_t const nq = m_top_left_block->size(); - std::size_t const ndelta = m_bottom_right_block->size(); - if (i >= nq && j < nq) { - std::ptrdiff_t d = j - i; - if (d > (std::ptrdiff_t)(size() / 2)) - d -= size(); - if (d < -(std::ptrdiff_t)(size() / 2)) - d += size(); - - if (d < -(std::ptrdiff_t)m_kl || d > (std::ptrdiff_t)m_ku) { - assert(std::fabs(aij) < 1e-20); - return; - } - if (d > (std::ptrdiff_t)0) { - m_bottom_left_block.h_view(i - nq, j) = aij; - } else { - m_bottom_left_block.h_view(i - nq, j - nq + ndelta + 1) = aij; - } - } else { - SplinesLinearProblem2x2Blocks::set_element(i, j, aij); - } - } - -private: - // @brief Compute the Schur complement delta - lambda*Q^-1*gamma. - void compute_schur_complement() - { - Kokkos::parallel_for( - "compute_schur_complement", - Kokkos::MDRangePolicy>( - {0, 0}, - {m_bottom_right_block->size(), m_bottom_right_block->size()}), - [&](const int i, const int j) { - double val = 0.0; - // Upper diagonals in lambda, lower diagonals in gamma - for (int l = 0; l < i + 1; ++l) { - val += m_bottom_left_block.h_view(i, l) * m_top_right_block.h_view(l, j); - } - // Lower diagonals in lambda, upper diagonals in gamma - for (int l = i + 1; l < m_bottom_right_block->size() + 1; ++l) { - int const l_full - = m_top_left_block->size() - 1 - m_bottom_right_block->size() + l; - val += m_bottom_left_block.h_view(i, l) - * m_top_right_block.h_view(l_full, j); - } - m_bottom_right_block - ->set_element(i, j, m_bottom_right_block->get_element(i, j) - val); - }); - } - -public: - /** - * @brief Compute y <- y - LinOp*x or y <- y - LinOp^t*x. - * - * [SHOULD BE PRIVATE (GPU programming limitation)] - * - * @param x - * @param y - * @param LinOp - * @param transpose - */ - void per_gemv_minus1_1( - MultiRHS const x, - MultiRHS const y, - Kokkos::View const - LinOp, - bool const transpose = false) const - { - /* - assert(!transpose && LinOp.extent(0) == y.extent(0) - || transpose && LinOp.extent(1) == y.extent(0)); - assert(!transpose && LinOp.extent(1) == x.extent(0) - || transpose && LinOp.extent(0) == x.extent(0)); - */ - assert(x.extent(1) == y.extent(1)); - - std::size_t const nq = m_top_left_block->size(); - std::size_t const ndelta = m_bottom_right_block->size(); - Kokkos::parallel_for( - "per_gemv_minus1_1", - Kokkos::TeamPolicy< - ExecSpace>((y.extent(0) + transpose) * y.extent(1), Kokkos::AUTO), - KOKKOS_LAMBDA( - const typename Kokkos::TeamPolicy::member_type& teamMember) { - const int i = teamMember.league_rank() / y.extent(1); - const int j = teamMember.league_rank() % y.extent(1); - - if (!transpose) { - double LinOpTimesX = 0.; - Kokkos::parallel_reduce( - Kokkos::TeamThreadRange(teamMember, i + 1), - [&](const int l, double& LinOpTimesX_tmp) { - LinOpTimesX_tmp += LinOp(i, l) * x(l, j); - }, - LinOpTimesX); - teamMember.team_barrier(); - double LinOpTimesX2 = 0.; - Kokkos::parallel_reduce( - Kokkos::TeamThreadRange(teamMember, i + 1, ndelta), - [&](const int l, double& LinOpTimesX_tmp) { - int const l_full = nq - 1 - ndelta + l; - LinOpTimesX_tmp += LinOp(i, l) * x(l_full, j); - }, - LinOpTimesX2); - if (teamMember.team_rank() == 0) { - y(i, j) -= LinOpTimesX + LinOpTimesX2; - } - } else { - // Lower diagonals in lambda - for (int l = 0; l < i; ++l) { - if (teamMember.team_rank() == 0) { - y(nq - 1 - ndelta + i, j) -= LinOp(l, i) * x(l, j); - } - } - /// Upper diagonals in lambda - for (int l = i; l < ndelta; ++l) { - if (teamMember.team_rank() == 0) { - y(i, j) -= LinOp(l, i) * x(l, j); - } - } - } - }); - } - - /** - * @brief Solve the multiple right-hand sides linear problem Ax=b or its transposed version A^tx=b inplace. - * - * The solver method is the one known as Schur complement method. It can be summarized as follow, - * starting with the pre-computed elements of the matrix: - * - * | Q | Q^-1*gamma | - * | lambda | delta - lambda*Q^-1*gamma | - * - * For the non-transposed case: - * - Solve inplace Q * x'1 = b1 (using the solver internal to Q). - * - Compute inplace b'2 = b2 - lambda*x'1. - * - Solve inplace (delta - lambda*Q^-1*gamma) * x2 = b'2. - * - Compute inplace x1 = x'1 - (delta - lambda*Q^-1*gamma)*x2. - * - * @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 b, bool const transpose) const override - { - assert(b.extent(0) == size()); - - MultiRHS b1 = Kokkos:: - subview(b, - std::pair(0, m_top_left_block->size()), - Kokkos::ALL); - MultiRHS b2 = Kokkos:: - subview(b, - std::pair(m_top_left_block->size(), b.extent(0)), - Kokkos::ALL); - if (!transpose) { - m_top_left_block->solve(b1); - per_gemv_minus1_1(b1, b2, m_bottom_left_block.d_view); - m_bottom_right_block->solve(b2); - gemv_minus1_1(b2, b1, m_top_right_block.d_view); - } else { - gemv_minus1_1(b1, b2, m_top_right_block.d_view, true); - m_bottom_right_block->solve(b2, true); - per_gemv_minus1_1(b2, b1, m_bottom_left_block.d_view, true); - m_top_left_block->solve(b1, true); - } - } -}; - -} // namespace ddc::detail