spline_builder.hpp

// Copyright (C) The DDC development team, see COPYRIGHT.md file
//
// SPDX-License-Identifier: MIT
 
#pragma once
 
#include <array>
#include <cassert>
#include <cstddef>
#include <memory>
#include <optional>
#include <stdexcept>
#include <string>
#include <tuple>
#include <utility>
 
#include <ddc/ddc.hpp>
 
#include <Kokkos_Core.hpp>
 
#include "deriv.hpp"
#include "integrals.hpp"
#include "math_tools.hpp"
#include "spline_builder_closures.hpp"
#include "splines_linear_problem.hpp"
#include "splines_linear_problem_maker.hpp"
#include "view.hpp"
 
namespace ddc {
 
/**
 * @brief An enum determining the backend solver of a SplineBuilder or SplineBuilder2d.
 *
 * 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)
    LAPACK ///< Enum member to identify the LAPACK-based solver (direct method)
};
 
/**
 * @brief A class for creating a spline approximation of a function.
 *
 * A class which contains an operator () which can be used to build a spline approximation
 * of a function. A spline approximation is represented by coefficients stored in a Chunk
 * of B-splines. The spline is constructed such that it respects the closure relations
 * SBCLower and SBCUpper, and it interpolates the function at the points on the interpolation_discrete_dimension
 * associated with interpolation_discrete_dimension_type.
 * @tparam ExecSpace The Kokkos execution space on which the spline approximation is performed.
 * @tparam MemorySpace The Kokkos memory space on which the data (interpolation function and splines coefficients) is stored.
 * @tparam BSplines The discrete dimension representing the B-splines.
 * @tparam InterpolationDDim The discrete dimension on which interpolation points are defined.
 * @tparam SBCLower The lower closure relation.
 * @tparam SBCUpper The upper closure relation.
 * @tparam Solver The SplineSolver giving the backend used to perform the spline approximation.
 */
template <
        class ExecSpace,
        class MemorySpace,
        class BSplines,
        class InterpolationDDim,
        ddc::SplineBuilderClosure SBCLower,
        ddc::SplineBuilderClosure SBCUpper,
        SplineSolver Solver>
class SplineBuilder
{
    static_assert(
            (BSplines::is_periodic() && (SBCLower == ddc::SplineBuilderClosure::PERIODIC)
             && (SBCUpper == ddc::SplineBuilderClosure::PERIODIC))
            || (!BSplines::is_periodic() && (SBCLower != ddc::SplineBuilderClosure::PERIODIC)
                && (SBCUpper != ddc::SplineBuilderClosure::PERIODIC)));
 
public:
    /// @brief The type of the Kokkos execution space used by this class.
    using exec_space = ExecSpace;
 
    /// @brief The type of the Kokkos memory space used by this class.
    using memory_space = MemorySpace;
 
    /// @brief The type of the interpolation continuous dimension (continuous dimension of interest) used by this class.
    using continuous_dimension_type = InterpolationDDim::continuous_dimension_type;
 
    /// @brief The type of the interpolation discrete dimension (discrete dimension of interest) used by this class.
    using interpolation_discrete_dimension_type = InterpolationDDim;
 
    /// @brief The discrete dimension representing the B-splines.
    using bsplines_type = BSplines;
 
    /// @brief The type of the Deriv dimension at the boundaries.
    using deriv_type = ddc::Deriv<continuous_dimension_type>;
 
    /// @brief The type of the domain for the 1D interpolation mesh used by this class.
    using interpolation_domain_type = ddc::DiscreteDomain<interpolation_discrete_dimension_type>;
 
    /**
     * @brief The type of the whole domain representing interpolation points.
     *
     * @tparam The batched discrete domain on which the interpolation points are defined.
     */
    template <concepts::discrete_domain BatchedInterpolationDDom>
    using batched_interpolation_domain_type = BatchedInterpolationDDom;
 
    /**
     * @brief The type of the batch domain (obtained by removing the dimension of interest
     * from the whole domain).
     *
     * @tparam The batched discrete domain on which the interpolation points are defined.
     *
     * Example: For batched_interpolation_domain_type = DiscreteDomain<X,Y,Z> and a dimension of interest Y,
     * this is DiscreteDomain<X,Z>
     */
    template <concepts::discrete_domain BatchedInterpolationDDom>
    using batch_domain_type = ddc::
            remove_dims_of_t<BatchedInterpolationDDom, interpolation_discrete_dimension_type>;
 
    /**
     * @brief The type of the whole spline domain (cartesian product of 1D spline domain
     * and batch domain) preserving the underlying memory layout (order of dimensions).
     *
     * @tparam The batched discrete domain on which the interpolation points are defined.
     *
     * Example: For batched_interpolation_domain_type = DiscreteDomain<X,Y,Z> and a dimension of interest Y
     * (associated to a B-splines tag BSplinesY), this is DiscreteDomain<X,BSplinesY,Z>.
     */
    template <concepts::discrete_domain BatchedInterpolationDDom>
    using batched_spline_domain_type = ddc::replace_dim_of_t<
            BatchedInterpolationDDom,
            interpolation_discrete_dimension_type,
            bsplines_type>;
 
private:
    /**
     * @brief The type of the whole spline domain (cartesian product of the 1D spline domain
     * and the batch domain) with 1D spline dimension being the leading dimension.
     *
     * @tparam The batched discrete domain on which the interpolation points are defined.
     *
     * Example: For batched_interpolation_domain_type = DiscreteDomain<X,Y,Z> and a dimension of interest Y
     * (associated to a B-splines tag BSplinesY), this is DiscreteDomain<BSplinesY,X,Z>.
     */
    template <concepts::discrete_domain BatchedInterpolationDDom>
    using batched_spline_tr_domain_type
            = ddc::detail::convert_type_seq_to_discrete_domain_t<ddc::type_seq_merge_t<
                    ddc::detail::TypeSeq<bsplines_type>,
                    ddc::type_seq_remove_t<
                            ddc::to_type_seq_t<BatchedInterpolationDDom>,
                            ddc::detail::TypeSeq<interpolation_discrete_dimension_type>>>>;
 
public:
    /**
     * @brief The type of the whole Deriv domain (cartesian product of 1D Deriv domain
     * and batch domain) preserving the underlying memory layout (order of dimensions).
     *
     * @tparam The batched discrete domain on which the interpolation points are defined.
     *
     * Example: For batched_interpolation_domain_type = DiscreteDomain<X,Y,Z> and a dimension of interest Y,
     * this is DiscreteDomain<X,Deriv<Y>,Z>
     */
    template <concepts::discrete_domain BatchedInterpolationDDom>
    using batched_derivs_domain_type = ddc::replace_dim_of_t<
            BatchedInterpolationDDom,
            interpolation_discrete_dimension_type,
            deriv_type>;
 
    /// @brief Indicates if the degree of the splines is odd or even.
    static constexpr bool s_odd = BSplines::degree() % 2;
 
    /// @brief The number of equations defining the closure relation at the lower bound.
    static constexpr int s_nbe_xmin = n_boundary_equations(SBCLower, BSplines::degree());
 
    /// @brief The number of equations defining the closure relation at the upper bound.
    static constexpr int s_nbe_xmax = n_boundary_equations(SBCUpper, BSplines::degree());
 
    /// @brief The number of input values defining the closure relation at the lower bound.
    static constexpr int s_nbv_xmin = SBCLower == SplineBuilderClosure::HOMOGENEOUS_HERMITE
                                              ? 0
                                              : n_boundary_equations(SBCLower, BSplines::degree());
 
    /// @brief The number of input values defining the closure relation at the upper bound.
    static constexpr int s_nbv_xmax = SBCUpper == SplineBuilderClosure::HOMOGENEOUS_HERMITE
                                              ? 0
                                              : n_boundary_equations(SBCUpper, BSplines::degree());
 
#if DDC_BUILD_DEPRECATED_CODE()
    /// @brief The boundary condition implemented at the lower bound.
    /// @deprecated Use s_sbc_xmin instead
    [[deprecated("Use s_sbc_xmin instead")]] static constexpr ddc::SplineBuilderClosure s_bc_xmin
            = SBCLower;
 
    /// @brief The boundary condition implemented at the upper bound.
    /// @deprecated Use s_sbc_xmax instead
    [[deprecated("Use s_sbc_xmax instead")]] static constexpr ddc::SplineBuilderClosure s_bc_xmax
            = SBCUpper;
#endif
 
    /// @brief The closure relation implemented at the lower bound.
    static constexpr ddc::SplineBuilderClosure s_sbc_xmin = SBCLower;
 
    /// @brief The closure relation implemented at the upper bound.
    static constexpr ddc::SplineBuilderClosure s_sbc_xmax = SBCUpper;
 
    /// @brief The SplineSolver giving the backend used to perform the spline approximation.
    static constexpr SplineSolver s_spline_solver = Solver;
 
private:
    interpolation_domain_type m_interpolation_domain;
 
    int m_offset = 0;
 
    double m_dx; // average cell size for normalization of derivatives
 
    // interpolator specific
    std::unique_ptr<ddc::detail::SplinesLinearProblem<exec_space>> m_matrix;
 
    std::string m_label;
 
    /// Calculate offset so that the matrix is diagonally dominant
    void compute_offset(interpolation_domain_type const& interpolation_domain, int& offset);
 
public:
    /**
     * @brief Build a SplineBuilder acting on interpolation_domain.
     *
     * @param label A label used to tag parallel regions and memory allocations for profiling.
     *
     * @param interpolation_domain The domain on which the interpolation points are defined.
     *
     * @param cols_per_chunk A parameter used by the slicer (internal to the solver) to define the size
     * of a chunk of right-hand sides of the linear problem to be computed in parallel (chunks are treated
     * by the linear solver one-after-the-other).
     * This value is optional. If no value is provided then the default value is chosen by the requested solver.
     *
     * @param preconditioner_max_block_size A parameter used by the slicer (internal to the solver) to
     * define the size of a block used by the Block-Jacobi preconditioner.
     * This value is optional. If no value is provided then the default value is chosen by the requested solver.
     *
     * @see MatrixSparse
     */
    explicit SplineBuilder(
            std::string label,
            interpolation_domain_type const& interpolation_domain,
            std::optional<std::size_t> cols_per_chunk = std::nullopt,
            std::optional<unsigned int> preconditioner_max_block_size = std::nullopt)
        : m_interpolation_domain(interpolation_domain)
        , m_dx((ddc::discrete_space<BSplines>().rmax() - ddc::discrete_space<BSplines>().rmin())
               / ddc::discrete_space<BSplines>().ncells())
        , m_label(std::move(label))
    {
        static_assert(
                ((SBCLower == SplineBuilderClosure::PERIODIC)
                 == (SBCUpper == SplineBuilderClosure::PERIODIC)),
                "Incompatible closure relations");
        check_valid_grid();
 
        compute_offset(this->interpolation_domain(), m_offset);
 
        // Calculate block sizes
        int lower_block_size;
        int upper_block_size;
        if constexpr (bsplines_type::is_uniform()) {
            upper_block_size = compute_block_sizes_uniform(SBCLower, s_nbe_xmin);
            lower_block_size = compute_block_sizes_uniform(SBCUpper, s_nbe_xmax);
        } else {
            upper_block_size = compute_block_sizes_non_uniform(SBCLower, s_nbe_xmin);
            lower_block_size = compute_block_sizes_non_uniform(SBCUpper, s_nbe_xmax);
        }
        allocate_matrix(
                lower_block_size,
                upper_block_size,
                cols_per_chunk,
                preconditioner_max_block_size);
    }
 
    /**
     * @brief Build a SplineBuilder acting on interpolation_domain.
     *
     * @param interpolation_domain The domain on which the interpolation points are defined.
     *
     * @param cols_per_chunk A parameter used by the slicer (internal to the solver) to define the size
     * of a chunk of right-hand sides of the linear problem to be computed in parallel (chunks are treated
     * by the linear solver one-after-the-other).
     * This value is optional. If no value is provided then the default value is chosen by the requested solver.
     *
     * @param preconditioner_max_block_size A parameter used by the slicer (internal to the solver) to
     * define the size of a block used by the Block-Jacobi preconditioner.
     * This value is optional. If no value is provided then the default value is chosen by the requested solver.
     *
     * @see MatrixSparse
     */
    explicit SplineBuilder(
            interpolation_domain_type const& interpolation_domain,
            std::optional<std::size_t> cols_per_chunk = std::nullopt,
            std::optional<unsigned int> preconditioner_max_block_size = std::nullopt)
        : SplineBuilder(
                  "no-label",
                  interpolation_domain,
                  cols_per_chunk,
                  preconditioner_max_block_size)
    {
    }
 
    /**
     * @brief Build a SplineBuilder acting on the interpolation domain contained by batched_interpolation_domain.
     *
     * @param label A label used to tag parallel regions and memory allocations for profiling.
     *
     * @param batched_interpolation_domain The whole domain on which the interpolation points are defined.
     *
     * @param cols_per_chunk A parameter used by the slicer (internal to the solver) to define the size
     * of a chunk of right-hand sides of the linear problem to be computed in parallel (chunks are treated
     * by the linear solver one-after-the-other).
     * This value is optional. If no value is provided then the default value is chosen by the requested solver.
     *
     * @param preconditioner_max_block_size A parameter used by the slicer (internal to the solver) to
     * define the size of a block used by the Block-Jacobi preconditioner.
     * This value is optional. If no value is provided then the default value is chosen by the requested solver.
     *
     * @see MatrixSparse
     */
    template <concepts::discrete_domain BatchedInterpolationDDom>
    explicit SplineBuilder(
            std::string label,
            BatchedInterpolationDDom const& batched_interpolation_domain,
            std::optional<std::size_t> cols_per_chunk = std::nullopt,
            std::optional<unsigned int> preconditioner_max_block_size = std::nullopt)
        : SplineBuilder(
                  std::move(label),
                  interpolation_domain_type(batched_interpolation_domain),
                  cols_per_chunk,
                  preconditioner_max_block_size)
    {
    }
 
    /**
     * @brief Build a SplineBuilder acting on the interpolation domain contained by batched_interpolation_domain.
     *
     * @param batched_interpolation_domain The whole domain on which the interpolation points are defined.
     *
     * @param cols_per_chunk A parameter used by the slicer (internal to the solver) to define the size
     * of a chunk of right-hand sides of the linear problem to be computed in parallel (chunks are treated
     * by the linear solver one-after-the-other).
     * This value is optional. If no value is provided then the default value is chosen by the requested solver.
     *
     * @param preconditioner_max_block_size A parameter used by the slicer (internal to the solver) to
     * define the size of a block used by the Block-Jacobi preconditioner.
     * This value is optional. If no value is provided then the default value is chosen by the requested solver.
     *
     * @see MatrixSparse
     */
    template <concepts::discrete_domain BatchedInterpolationDDom>
    explicit SplineBuilder(
            BatchedInterpolationDDom const& batched_interpolation_domain,
            std::optional<std::size_t> cols_per_chunk = std::nullopt,
            std::optional<unsigned int> preconditioner_max_block_size = std::nullopt)
        : SplineBuilder(
                  "no-label",
                  interpolation_domain_type(batched_interpolation_domain),
                  cols_per_chunk,
                  preconditioner_max_block_size)
    {
    }
 
 
    /// @brief Copy-constructor is deleted.
    SplineBuilder(SplineBuilder const& x) = delete;
 
    /** @brief Move-constructs.
     *
     * @param x An rvalue to another SplineBuilder.
     */
    SplineBuilder(SplineBuilder&& x) = default;
 
    /// @brief Destructs.
    ~SplineBuilder() = default;
 
    /// @brief Copy-assignment is deleted.
    SplineBuilder& operator=(SplineBuilder const& x) = delete;
 
    /** @brief Move-assigns.
     *
     * @param x An rvalue to another SplineBuilder.
     * @return A reference to this object.
     */
    SplineBuilder& operator=(SplineBuilder&& x) = default;
 
    /**
     * @brief Get the domain for the 1D interpolation mesh used by this class.
     *
     * This is 1D because it is defined along the dimension of interest.
     *
     * @return The 1D domain for the interpolation mesh.
     */
    interpolation_domain_type interpolation_domain() const noexcept
    {
        return m_interpolation_domain;
    }
 
    /**
     * @brief Get the whole domain representing interpolation points.
     *
     * Values of the function must be provided on this domain in order
     * to build a spline representation of the function (cartesian product of 1D interpolation_domain and batch_domain).
     *
     * @param batched_interpolation_domain The whole domain on which the interpolation points are defined.
     *
     * @return The domain for the interpolation mesh.
     */
    template <concepts::discrete_domain BatchedInterpolationDDom>
    BatchedInterpolationDDom batched_interpolation_domain(
            BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept
    {
        assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain));
        return batched_interpolation_domain;
    }
 
    /**
     * @brief Get the batch domain.
     *
     * Obtained by removing the dimension of interest from the whole interpolation domain.
     *
     * @param batched_interpolation_domain The whole domain on which the interpolation points are defined.
     *
     * @return The batch domain.
     */
    template <class BatchedInterpolationDDom>
    batch_domain_type<BatchedInterpolationDDom> batch_domain(
            BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept
    {
        assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain));
        return ddc::remove_dims_of(batched_interpolation_domain, interpolation_domain());
    }
 
    /**
     * @brief Get the 1D domain on which spline coefficients are defined.
     *
     * The 1D spline domain corresponding to the dimension of interest.
     *
     * @return The 1D domain for the spline coefficients.
     */
    ddc::DiscreteDomain<bsplines_type> spline_domain() const noexcept
    {
        return ddc::discrete_space<bsplines_type>().full_domain();
    }
 
    /**
     * @brief Get the whole domain on which spline coefficients are defined.
     *
     * Spline approximations (spline-transformed functions) are computed on this domain.
     *
     * @param batched_interpolation_domain The whole domain on which the interpolation points are defined.
     *
     * @return The domain for the spline coefficients.
     */
    template <class BatchedInterpolationDDom>
    batched_spline_domain_type<BatchedInterpolationDDom> batched_spline_domain(
            BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept
    {
        assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain));
        return ddc::replace_dim_of<
                interpolation_discrete_dimension_type,
                bsplines_type>(batched_interpolation_domain, spline_domain());
    }
 
private:
    /**
     * @brief Get the whole domain on which spline coefficients are defined, with the dimension of interest being the leading dimension.
     *
     * This is used internally due to solver limitation and because it may be beneficial to computation performance. For LAPACK backend and non-periodic closure relation, we are using SplinesLinearSolver3x3Blocks which requires upper_block_size additional rows for internal operations.
     *
     * @param batched_interpolation_domain The whole domain on which the interpolation points are defined.
     *
     * @return The (transposed) domain for the spline coefficients.
     */
    template <class BatchedInterpolationDDom>
    batched_spline_tr_domain_type<BatchedInterpolationDDom> batched_spline_tr_domain(
            BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept
    {
        assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain));
        return batched_spline_tr_domain_type<BatchedInterpolationDDom>(
                ddc::replace_dim_of<bsplines_type, bsplines_type>(
                        batched_spline_domain(batched_interpolation_domain),
                        ddc::DiscreteDomain<bsplines_type>(
                                ddc::DiscreteElement<bsplines_type>(0),
                                ddc::DiscreteVector<bsplines_type>(
                                        m_matrix->required_number_of_rhs_rows()))));
    }
 
public:
    /**
     * @brief Get the whole domain on which derivatives on lower boundary are defined.
     *
     * This is only used with SplineBuilderClosure::HERMITE closure relations.
     *
     * @param batched_interpolation_domain The whole domain on which the interpolation points are defined.
     *
     * @return The domain for the Derivs values.
     */
    template <class BatchedInterpolationDDom>
    batched_derivs_domain_type<BatchedInterpolationDDom> batched_derivs_xmin_domain(
            BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept
    {
        assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain));
        return ddc::replace_dim_of<interpolation_discrete_dimension_type, deriv_type>(
                batched_interpolation_domain,
                ddc::DiscreteDomain<deriv_type>(
                        ddc::DiscreteElement<deriv_type>(1),
                        ddc::DiscreteVector<deriv_type>(s_nbv_xmin)));
    }
 
    /**
     * @brief Get the whole domain on which derivatives on upper boundary are defined.
     *
     * This is only used with SplineBuilderClosure::HERMITE closure relations.
     *
     * @param batched_interpolation_domain The whole domain on which the interpolation points are defined.
     *
     * @return The domain for the Derivs values.
     */
    template <class BatchedInterpolationDDom>
    batched_derivs_domain_type<BatchedInterpolationDDom> batched_derivs_xmax_domain(
            BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept
    {
        assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain));
        return ddc::replace_dim_of<interpolation_discrete_dimension_type, deriv_type>(
                batched_interpolation_domain,
                ddc::DiscreteDomain<deriv_type>(
                        ddc::DiscreteElement<deriv_type>(1),
                        ddc::DiscreteVector<deriv_type>(s_nbv_xmax)));
    }
 
    /**
     * @brief Compute a spline approximation of a function.
     *
     * Use the values of a function (defined on
     * SplineBuilder::batched_interpolation_domain) and the derivatives of the
     * function at the boundaries (in the case of SplineBuilderClosure::HERMITE only, defined
     * on SplineBuilder::batched_derivs_xmin_domain and SplineBuilder::batched_derivs_xmax_domain)
     * to calculate a spline approximation of this function.
     *
     * The spline approximation is stored as a ChunkSpan of coefficients
     * associated with B-splines.
     *
     * @param[out] spline The coefficients of the spline computed by this SplineBuilder.
     * @param[in] vals The values of the function on the interpolation mesh.
     * @param[in] derivs_xmin The values of the derivatives at the lower boundary
     * (used only with SplineBuilderClosure::HERMITE lower closure relation).
     * @param[in] derivs_xmax The values of the derivatives at the upper boundary
     * (used only with SplineBuilderClosure::HERMITE upper closure relation).
     */
    template <class Layout, class BatchedInterpolationDDom>
    void operator()(
            ddc::ChunkSpan<
                    double,
                    batched_spline_domain_type<BatchedInterpolationDDom>,
                    Layout,
                    memory_space> spline,
            ddc::ChunkSpan<double const, BatchedInterpolationDDom, Layout, memory_space> vals,
            std::optional<ddc::ChunkSpan<
                    double const,
                    batched_derivs_domain_type<BatchedInterpolationDDom>,
                    Layout,
                    memory_space>> derivs_xmin
            = std::nullopt,
            std::optional<ddc::ChunkSpan<
                    double const,
                    batched_derivs_domain_type<BatchedInterpolationDDom>,
                    Layout,
                    memory_space>> derivs_xmax
            = std::nullopt) const;
 
    /**
     * @brief Compute the quadrature coefficients associated to the b-splines used by this SplineBuilder.
     *
     * Those coefficients can be used to perform integration way faster than SplineEvaluator::integrate().
     *
     * This function solves matrix equation A^t*Q=integral_bsplines. In case of HERMITE closure relations,
     * integral_bsplines contains the integral coefficients at the boundaries, and Q thus has to
     * be split in three parts (quadrature coefficients for the derivatives at lower boundary,
     * for the values inside the domain and for the derivatives at upper boundary).
     *
     * A discrete function f can then be integrated using sum_j Q_j*f_j for j in interpolation_domain.
     * If closure relation is HERMITE, sum_j Qderiv_j*(d^j f/dx^j) for j in derivs_domain
     * must be added at the boundary.
     *
     * Please refer to section 2.8.1 of Emily's Bourne phd (https://theses.fr/2022AIXM0412) for more information and to
     * the (Non)PeriodicSplineBuilderTest for example usage to compute integrals.
     *
     * @tparam OutMemorySpace The Kokkos::MemorySpace on which the quadrature coefficients are be returned
     * (but they are computed on ExecSpace then copied).
     *
     * @return A tuple containing the three Chunks containing the quadrature coefficients (if HERMITE
     * is not used, first and third are empty).
     */
    template <class OutMemorySpace = MemorySpace>
    std::tuple<
            ddc::Chunk<
                    double,
                    ddc::DiscreteDomain<
                            ddc::Deriv<typename InterpolationDDim::continuous_dimension_type>>,
                    ddc::KokkosAllocator<double, OutMemorySpace>>,
            ddc::Chunk<
                    double,
                    ddc::DiscreteDomain<InterpolationDDim>,
                    ddc::KokkosAllocator<double, OutMemorySpace>>,
            ddc::Chunk<
                    double,
                    ddc::DiscreteDomain<
                            ddc::Deriv<typename InterpolationDDim::continuous_dimension_type>>,
                    ddc::KokkosAllocator<double, OutMemorySpace>>>
    quadrature_coefficients() const;
 
private:
    static int compute_block_sizes_uniform(ddc::SplineBuilderClosure bound_cond, int nbc);
 
    static int compute_block_sizes_non_uniform(ddc::SplineBuilderClosure bound_cond, int nbc);
 
    void allocate_matrix(
            int lower_block_size,
            int upper_block_size,
            std::optional<std::size_t> cols_per_chunk = std::nullopt,
            std::optional<unsigned int> preconditioner_max_block_size = std::nullopt);
 
    void build_matrix_system();
 
    void check_valid_grid();
 
    template <class KnotElement>
    static void check_n_points_in_cell(int n_points_in_cell, KnotElement current_cell_end_idx);
};
 
template <
        class ExecSpace,
        class MemorySpace,
        class BSplines,
        class InterpolationDDim,
        ddc::SplineBuilderClosure SBCLower,
        ddc::SplineBuilderClosure SBCUpper,
        SplineSolver Solver>
void SplineBuilder<
        ExecSpace,
        MemorySpace,
        BSplines,
        InterpolationDDim,
        SBCLower,
        SBCUpper,
        Solver>::compute_offset(interpolation_domain_type const& interpolation_domain, int& offset)
{
    if constexpr (bsplines_type::is_periodic()) {
        // Calculate offset so that the matrix is diagonally dominant
        std::array<double, bsplines_type::degree() + 1> values_ptr;
        Kokkos::mdspan<double, Kokkos::extents<std::size_t, bsplines_type::degree() + 1>> const
                values(values_ptr.data());
        ddc::DiscreteElement<interpolation_discrete_dimension_type> start(
                interpolation_domain.front());
        auto jmin = ddc::discrete_space<BSplines>()
                            .eval_basis(values, ddc::coordinate(start + BSplines::degree()));
        if constexpr (bsplines_type::degree() % 2 == 0) {
            offset = jmin.uid() - start.uid() + bsplines_type::degree() / 2 - BSplines::degree();
        } else {
            int const mid = bsplines_type::degree() / 2;
            offset = jmin.uid() - start.uid()
                     + (DDC_MDSPAN_ACCESS_OP(values, mid) > DDC_MDSPAN_ACCESS_OP(values, mid + 1)
                                ? mid
                                : mid + 1)
                     - BSplines::degree();
        }
    } else {
        offset = 0;
    }
}
 
template <
        class ExecSpace,
        class MemorySpace,
        class BSplines,
        class InterpolationDDim,
        ddc::SplineBuilderClosure SBCLower,
        ddc::SplineBuilderClosure SBCUpper,
        SplineSolver Solver>
int SplineBuilder<ExecSpace, MemorySpace, BSplines, InterpolationDDim, SBCLower, SBCUpper, Solver>::
        compute_block_sizes_uniform(ddc::SplineBuilderClosure const bound_cond, int const nbc)
{
    if (bound_cond == ddc::SplineBuilderClosure::PERIODIC) {
        return static_cast<int>(bsplines_type::degree()) / 2;
    }
 
    if (bound_cond == ddc::SplineBuilderClosure::HERMITE
        || bound_cond == ddc::SplineBuilderClosure::HOMOGENEOUS_HERMITE) {
        return nbc;
    }
 
    if (bound_cond == ddc::SplineBuilderClosure::GREVILLE) {
        return static_cast<int>(bsplines_type::degree()) - 1;
    }
 
    throw std::runtime_error("ddc::SplineBuilderClosure not handled");
}
 
template <
        class ExecSpace,
        class MemorySpace,
        class BSplines,
        class InterpolationDDim,
        ddc::SplineBuilderClosure SBCLower,
        ddc::SplineBuilderClosure SBCUpper,
        SplineSolver Solver>
int SplineBuilder<ExecSpace, MemorySpace, BSplines, InterpolationDDim, SBCLower, SBCUpper, Solver>::
        compute_block_sizes_non_uniform(ddc::SplineBuilderClosure const bound_cond, int const nbc)
{
    if (bound_cond == ddc::SplineBuilderClosure::PERIODIC
        || bound_cond == ddc::SplineBuilderClosure::GREVILLE) {
        return static_cast<int>(bsplines_type::degree()) - 1;
    }
 
    if (bound_cond == ddc::SplineBuilderClosure::HERMITE
        || bound_cond == ddc::SplineBuilderClosure::HOMOGENEOUS_HERMITE) {
        return nbc + 1;
    }
 
    throw std::runtime_error("ddc::SplineBuilderClosure not handled");
}
 
template <
        class ExecSpace,
        class MemorySpace,
        class BSplines,
        class InterpolationDDim,
        ddc::SplineBuilderClosure SBCLower,
        ddc::SplineBuilderClosure SBCUpper,
        SplineSolver Solver>
void SplineBuilder<
        ExecSpace,
        MemorySpace,
        BSplines,
        InterpolationDDim,
        SBCLower,
        SBCUpper,
        Solver>::
        allocate_matrix(
                [[maybe_unused]] int lower_block_size,
                [[maybe_unused]] int upper_block_size,
                std::optional<std::size_t> cols_per_chunk,
                std::optional<unsigned int> preconditioner_max_block_size)
{
    // Special case: linear spline
    // No need for matrix assembly
    // (disabled)
    // if constexpr (bsplines_type::degree() == 1) {
    //     return;
    // }
 
    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()) {
            m_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 {
            m_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) {
        m_matrix = ddc::detail::SplinesLinearProblemMaker::make_new_sparse<ExecSpace>(
                ddc::discrete_space<BSplines>().nbasis(),
                cols_per_chunk,
                preconditioner_max_block_size);
    }
 
    build_matrix_system();
 
    m_matrix->setup_solver();
}
 
template <
        class ExecSpace,
        class MemorySpace,
        class BSplines,
        class InterpolationDDim,
        ddc::SplineBuilderClosure SBCLower,
        ddc::SplineBuilderClosure SBCUpper,
        SplineSolver Solver>
void SplineBuilder<
        ExecSpace,
        MemorySpace,
        BSplines,
        InterpolationDDim,
        SBCLower,
        SBCUpper,
        Solver>::build_matrix_system()
{
    // Hermite closure relations at xmin, if any
    if constexpr (
            SBCLower == ddc::SplineBuilderClosure::HERMITE
            || SBCLower == ddc::SplineBuilderClosure::HOMOGENEOUS_HERMITE) {
        std::array<double, (bsplines_type::degree() / 2 + 1) * (bsplines_type::degree() + 1)>
                derivs_ptr;
        ddc::DSpan2D const
                derivs(derivs_ptr.data(),
                       bsplines_type::degree() + 1,
                       bsplines_type::degree() / 2 + 1);
        ddc::discrete_space<BSplines>().eval_basis_and_n_derivs(
                derivs,
                ddc::discrete_space<BSplines>().rmin(),
                s_nbe_xmin);
 
        // In order to improve the condition number of the matrix, we normalize
        // all derivatives by multiplying the i-th derivative by dx^i
        for (std::size_t i = 0; i < bsplines_type::degree() + 1; ++i) {
            for (std::size_t j = 1; j < bsplines_type::degree() / 2 + 1; ++j) {
                DDC_MDSPAN_ACCESS_OP(derivs, i, j) *= ddc::detail::ipow(m_dx, j);
            }
        }
 
        if constexpr (s_nbe_xmin > 0) {
            // iterate only to deg as last bspline is 0
            for (std::size_t i = 0; i < s_nbe_xmin; ++i) {
                for (std::size_t j = 0; j < bsplines_type::degree(); ++j) {
                    m_matrix->set_element(i, j, DDC_MDSPAN_ACCESS_OP(derivs, j, i + s_odd));
                }
            }
        }
    }
 
    // Interpolation points
    std::array<double, bsplines_type::degree() + 1> values_ptr;
    Kokkos::mdspan<double, Kokkos::extents<std::size_t, bsplines_type::degree() + 1>> const values(
            values_ptr.data());
 
    int start = interpolation_domain().front().uid();
    ddc::host_for_each(interpolation_domain(), [&](auto ix) {
        auto jmin = ddc::discrete_space<BSplines>().eval_basis(
                values,
                ddc::coordinate(ddc::DiscreteElement<interpolation_discrete_dimension_type>(ix)));
        for (std::size_t s = 0; s < bsplines_type::degree() + 1; ++s) {
            int const j = ddc::detail::
                    modulo(int(jmin.uid() - m_offset + s),
                           static_cast<int>(ddc::discrete_space<BSplines>().nbasis()));
            m_matrix->set_element(
                    ix.uid() - start + s_nbe_xmin,
                    j,
                    DDC_MDSPAN_ACCESS_OP(values, s));
        }
    });
 
    // Hermite closure relations at xmax, if any
    if constexpr (
            SBCUpper == ddc::SplineBuilderClosure::HERMITE
            || SBCUpper == ddc::SplineBuilderClosure::HOMOGENEOUS_HERMITE) {
        std::array<double, (bsplines_type::degree() / 2 + 1) * (bsplines_type::degree() + 1)>
                derivs_ptr;
        Kokkos::mdspan<
                double,
                Kokkos::extents<
                        std::size_t,
                        bsplines_type::degree() + 1,
                        bsplines_type::degree() / 2 + 1>> const derivs(derivs_ptr.data());
 
        ddc::discrete_space<BSplines>().eval_basis_and_n_derivs(
                derivs,
                ddc::discrete_space<BSplines>().rmax(),
                s_nbe_xmax);
 
        // In order to improve the condition number of the matrix, we normalize
        // all derivatives by multiplying the i-th derivative by dx^i
        for (std::size_t i = 0; i < bsplines_type::degree() + 1; ++i) {
            for (std::size_t j = 1; j < bsplines_type::degree() / 2 + 1; ++j) {
                DDC_MDSPAN_ACCESS_OP(derivs, i, j) *= ddc::detail::ipow(m_dx, j);
            }
        }
 
        if constexpr (s_nbe_xmax > 0) {
            int const i0 = ddc::discrete_space<BSplines>().nbasis() - s_nbe_xmax;
            int const j0 = ddc::discrete_space<BSplines>().nbasis() - bsplines_type::degree();
            for (std::size_t j = 0; j < bsplines_type::degree(); ++j) {
                for (std::size_t i = 0; i < s_nbe_xmax; ++i) {
                    m_matrix->set_element(
                            i0 + i,
                            j0 + j,
                            DDC_MDSPAN_ACCESS_OP(derivs, j + 1, i + s_odd));
                }
            }
        }
    }
}
 
template <
        class ExecSpace,
        class MemorySpace,
        class BSplines,
        class InterpolationDDim,
        ddc::SplineBuilderClosure SBCLower,
        ddc::SplineBuilderClosure SBCUpper,
        SplineSolver Solver>
template <class Layout, class BatchedInterpolationDDom>
void SplineBuilder<
        ExecSpace,
        MemorySpace,
        BSplines,
        InterpolationDDim,
        SBCLower,
        SBCUpper,
        Solver>::
operator()(
        ddc::ChunkSpan<
                double,
                batched_spline_domain_type<BatchedInterpolationDDom>,
                Layout,
                memory_space> spline,
        ddc::ChunkSpan<double const, BatchedInterpolationDDom, Layout, memory_space> vals,
        std::optional<ddc::ChunkSpan<
                double const,
                batched_derivs_domain_type<BatchedInterpolationDDom>,
                Layout,
                memory_space>> const derivs_xmin,
        std::optional<ddc::ChunkSpan<
                double const,
                batched_derivs_domain_type<BatchedInterpolationDDom>,
                Layout,
                memory_space>> const derivs_xmax) const
{
    auto const batched_interpolation_domain = vals.domain();
 
    assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain));
    assert(batch_domain_type<BatchedInterpolationDDom>(batched_interpolation_domain)
           == batch_domain_type<BatchedInterpolationDDom>(spline.domain()));
 
    if (batch_domain(batched_interpolation_domain).empty()) {
        return;
    }
 
    assert(vals.template extent<interpolation_discrete_dimension_type>()
           == ddc::discrete_space<bsplines_type>().nbasis() - s_nbe_xmin - s_nbe_xmax);
 
    if constexpr (SBCLower == SplineBuilderClosure::HERMITE) {
        assert(ddc::DiscreteElement<deriv_type>(derivs_xmin->domain().front()).uid() == s_odd);
        assert(derivs_xmin.has_value() || s_nbe_xmin == 0);
    } else {
        assert(!derivs_xmin.has_value() || derivs_xmin->template extent<deriv_type>() == 0);
    }
    if constexpr (SBCUpper == SplineBuilderClosure::HERMITE) {
        assert(ddc::DiscreteElement<deriv_type>(derivs_xmax->domain().front()).uid() == s_odd);
        assert(derivs_xmax.has_value() || s_nbe_xmax == 0);
    } else {
        assert(!derivs_xmax.has_value() || derivs_xmax->template extent<deriv_type>() == 0);
    }
 
    // Hermite closure relations at xmin, if any
    // NOTE: For consistency with the linear system, the i-th derivative
    //       provided by the user must be multiplied by dx^i
    if constexpr (SBCLower == SplineBuilderClosure::HERMITE) {
        assert(derivs_xmin->template extent<deriv_type>() == s_nbe_xmin);
        auto derivs_xmin_values = *derivs_xmin;
        auto const dx_proxy = m_dx;
        auto const odd_proxy = s_odd;
        ddc::parallel_for_each(
                "ddc_splines_hermite_compute_lower_coefficients",
                exec_space(),
                batch_domain(batched_interpolation_domain),
                KOKKOS_LAMBDA(
                        batch_domain_type<BatchedInterpolationDDom>::discrete_element_type j) {
                    for (int i = 0; i < s_nbe_xmin; ++i) {
                        spline(ddc::DiscreteElement<bsplines_type>(i), j)
                                = derivs_xmin_values(
                                          ddc::DiscreteElement<deriv_type>(i + odd_proxy),
                                          j)
                                  * ddc::detail::ipow(dx_proxy, i + odd_proxy);
                    }
                });
    } else if constexpr (SBCLower == SplineBuilderClosure::HOMOGENEOUS_HERMITE) {
        ddc::DiscreteDomain<bsplines_type> const dx_splines(
                ddc::DiscreteElement<bsplines_type>(0),
                ddc::DiscreteVector<bsplines_type>(s_nbe_xmin));
        batched_spline_domain_type<BatchedInterpolationDDom> const
                dx_spline_domain(dx_splines, batch_domain(batched_interpolation_domain));
        ddc::parallel_fill(exec_space(), spline[dx_spline_domain], 0.0);
    }
 
    // Fill spline with vals (to work in spline afterward and preserve vals)
    ddc::parallel_fill(
            exec_space(),
            spline[ddc::DiscreteDomain<bsplines_type>(
                    ddc::DiscreteElement<bsplines_type>(s_nbe_xmin),
                    ddc::DiscreteVector<bsplines_type>(m_offset))],
            0.);
    // NOTE: We rely on Kokkos::deep_copy because ddc::parallel_deepcopy do not support
    //       different domain-typed Chunks.
    Kokkos::deep_copy(
            exec_space(),
            spline[ddc::DiscreteDomain<bsplines_type>(
                           ddc::DiscreteElement<bsplines_type>(s_nbe_xmin + m_offset),
                           ddc::DiscreteVector<bsplines_type>(static_cast<std::size_t>(
                                   vals.domain()
                                           .template extent<
                                                   interpolation_discrete_dimension_type>())))]
                    .allocation_kokkos_view(),
            vals.allocation_kokkos_view());
 
 
 
    // Hermite closure relations at xmax, if any
    // NOTE: For consistency with the linear system, the i-th derivative
    //       provided by the user must be multiplied by dx^i
    auto const& nbasis_proxy = ddc::discrete_space<bsplines_type>().nbasis();
    if constexpr (SBCUpper == SplineBuilderClosure::HERMITE) {
        assert(derivs_xmax->template extent<deriv_type>() == s_nbe_xmax);
        auto derivs_xmax_values = *derivs_xmax;
        auto const dx_proxy = m_dx;
        auto const odd_proxy = s_odd;
        ddc::parallel_for_each(
                "ddc_splines_hermite_compute_upper_coefficients",
                exec_space(),
                batch_domain(batched_interpolation_domain),
                KOKKOS_LAMBDA(
                        batch_domain_type<BatchedInterpolationDDom>::discrete_element_type j) {
                    for (int i = 0; i < s_nbe_xmax; ++i) {
                        spline(ddc::DiscreteElement<bsplines_type>(nbasis_proxy - s_nbe_xmax + i),
                               j)
                                = derivs_xmax_values(
                                          ddc::DiscreteElement<deriv_type>(i + odd_proxy),
                                          j)
                                  * ddc::detail::ipow(dx_proxy, i + odd_proxy);
                    }
                });
    } else if constexpr (SBCUpper == SplineBuilderClosure::HOMOGENEOUS_HERMITE) {
        ddc::DiscreteDomain<bsplines_type> const dx_splines(
                ddc::DiscreteElement<bsplines_type>(nbasis_proxy - s_nbe_xmax),
                ddc::DiscreteVector<bsplines_type>(s_nbe_xmax));
        batched_spline_domain_type<BatchedInterpolationDDom> const
                dx_spline_domain(dx_splines, batch_domain(batched_interpolation_domain));
        ddc::parallel_fill(exec_space(), spline[dx_spline_domain], 0.0);
    }
 
    // Allocate and fill a transposed version of spline in order to get dimension of interest as last dimension (optimal for GPU, necessary for Ginkgo). Also select only relevant rows in case of periodic boundaries
    auto const& offset_proxy = m_offset;
    ddc::Chunk spline_tr_alloc(
            m_label + " > spline_tr (ddc::SplineBuilder::operator())",
            batched_spline_tr_domain(batched_interpolation_domain),
            ddc::KokkosAllocator<double, memory_space>());
    ddc::ChunkSpan const spline_tr = spline_tr_alloc.span_view();
    ddc::parallel_for_each(
            m_label + " > ddc_splines_transpose_rhs",
            exec_space(),
            batch_domain(batched_interpolation_domain),
            KOKKOS_LAMBDA(
                    batch_domain_type<BatchedInterpolationDDom>::discrete_element_type const j) {
                for (std::size_t i = 0; i < nbasis_proxy; ++i) {
                    spline_tr(ddc::DiscreteElement<bsplines_type>(i), j)
                            = spline(ddc::DiscreteElement<bsplines_type>(i + offset_proxy), j);
                }
            });
    // Create a 2D Kokkos::View to manage spline_tr as a matrix
    Kokkos::View<double**, Kokkos::LayoutRight, exec_space> const bcoef_section(
            spline_tr.data_handle(),
            static_cast<std::size_t>(spline_tr.template extent<bsplines_type>()),
            batch_domain(batched_interpolation_domain).size());
    // Compute spline coef
    m_matrix->solve(bcoef_section, false);
    // Transpose back spline_tr into spline.
    ddc::parallel_for_each(
            m_label + " > ddc_splines_transpose_back_rhs",
            exec_space(),
            batch_domain(batched_interpolation_domain),
            KOKKOS_LAMBDA(
                    batch_domain_type<BatchedInterpolationDDom>::discrete_element_type const j) {
                for (std::size_t i = 0; i < nbasis_proxy; ++i) {
                    spline(ddc::DiscreteElement<bsplines_type>(i + offset_proxy), j)
                            = spline_tr(ddc::DiscreteElement<bsplines_type>(i), j);
                }
            });
 
    // Duplicate the lower spline coefficients to the upper side in case of periodic boundaries
    if (bsplines_type::is_periodic()) {
        ddc::parallel_for_each(
                m_label + " > ddc_splines_periodic_rows_duplicate_rhs",
                exec_space(),
                batch_domain(batched_interpolation_domain),
                KOKKOS_LAMBDA(
                        batch_domain_type<BatchedInterpolationDDom>::discrete_element_type const
                                j) {
                    if (offset_proxy != 0) {
                        for (int i = 0; i < offset_proxy; ++i) {
                            spline(ddc::DiscreteElement<bsplines_type>(i), j) = spline(
                                    ddc::DiscreteElement<bsplines_type>(nbasis_proxy + i),
                                    j);
                        }
                        for (std::size_t i = offset_proxy; i < bsplines_type::degree(); ++i) {
                            spline(ddc::DiscreteElement<bsplines_type>(nbasis_proxy + i), j)
                                    = spline(ddc::DiscreteElement<bsplines_type>(i), j);
                        }
                    }
                    for (std::size_t i(0); i < bsplines_type::degree(); ++i) {
                        ddc::DiscreteElement<bsplines_type> const i_start(i);
                        ddc::DiscreteElement<bsplines_type> const i_end(nbasis_proxy + i);
 
                        spline(i_end, j) = spline(i_start, j);
                    }
                });
    }
}
 
template <
        class ExecSpace,
        class MemorySpace,
        class BSplines,
        class InterpolationDDim,
        ddc::SplineBuilderClosure SBCLower,
        ddc::SplineBuilderClosure SBCUpper,
        SplineSolver Solver>
template <class OutMemorySpace>
std::tuple<
        ddc::Chunk<
                double,
                ddc::DiscreteDomain<
                        ddc::Deriv<typename InterpolationDDim::continuous_dimension_type>>,
                ddc::KokkosAllocator<double, OutMemorySpace>>,
        ddc::Chunk<
                double,
                ddc::DiscreteDomain<InterpolationDDim>,
                ddc::KokkosAllocator<double, OutMemorySpace>>,
        ddc::Chunk<
                double,
                ddc::DiscreteDomain<
                        ddc::Deriv<typename InterpolationDDim::continuous_dimension_type>>,
                ddc::KokkosAllocator<double, OutMemorySpace>>>
SplineBuilder<ExecSpace, MemorySpace, BSplines, InterpolationDDim, SBCLower, SBCUpper, Solver>::
        quadrature_coefficients() const
{
    // Compute integrals of bsplines
    ddc::Chunk integral_bsplines(spline_domain(), ddc::KokkosAllocator<double, MemorySpace>());
    ddc::integrals(ExecSpace(), integral_bsplines.span_view());
 
    // Remove additional B-splines in the periodic case (cf. UniformBSplines::full_domain() documentation)
    ddc::ChunkSpan const integral_bsplines_without_periodic_additional_bsplines
            = integral_bsplines[spline_domain().take_first(
                    ddc::DiscreteVector<bsplines_type>(m_matrix->size()))];
 
    // Allocate mirror with additional rows (cf. SplinesLinearProblem3x3Blocks documentation)
    Kokkos::View<double**, Kokkos::LayoutRight, MemorySpace> const
            integral_bsplines_mirror_with_additional_allocation(
                    m_label + " > integral_bsplines_mirror_with_additional_allocation",
                    m_matrix->required_number_of_rhs_rows(),
                    1);
 
    // Extract relevant subview
    Kokkos::View<double*, Kokkos::LayoutRight, MemorySpace> const integral_bsplines_mirror
            = Kokkos::
                    subview(integral_bsplines_mirror_with_additional_allocation,
                            std::
                                    pair {static_cast<std::size_t>(0),
                                          integral_bsplines_without_periodic_additional_bsplines
                                                  .size()},
                            0);
 
    // Solve matrix equation A^t*X=integral_bsplines
    Kokkos::deep_copy(
            integral_bsplines_mirror,
            integral_bsplines_without_periodic_additional_bsplines.allocation_kokkos_view());
    m_matrix->solve(integral_bsplines_mirror_with_additional_allocation, true);
    Kokkos::deep_copy(
            integral_bsplines_without_periodic_additional_bsplines.allocation_kokkos_view(),
            integral_bsplines_mirror);
 
    // Slice into three ChunkSpan corresponding to lower derivatives, function values and upper derivatives
    ddc::ChunkSpan const coefficients_derivs_xmin
            = integral_bsplines_without_periodic_additional_bsplines[spline_domain().take_first(
                    ddc::DiscreteVector<bsplines_type>(s_nbv_xmin))];
    ddc::ChunkSpan const coefficients = integral_bsplines_without_periodic_additional_bsplines
            [spline_domain()
                     .remove_first(ddc::DiscreteVector<bsplines_type>(s_nbv_xmin))
                     .take_first(
                             ddc::DiscreteVector<bsplines_type>(
                                     ddc::discrete_space<bsplines_type>().nbasis() - s_nbv_xmin
                                     - s_nbv_xmax))];
    ddc::ChunkSpan const coefficients_derivs_xmax
            = integral_bsplines_without_periodic_additional_bsplines
                    [spline_domain()
                             .remove_first(
                                     ddc::DiscreteVector<bsplines_type>(
                                             s_nbv_xmin + coefficients.size()))
                             .take_first(ddc::DiscreteVector<bsplines_type>(s_nbv_xmax))];
 
    // Multiply derivatives coefficients by dx^n
    auto const dx_proxy = m_dx;
    auto const odd_proxy = s_odd;
    ddc::parallel_for_each(
            exec_space(),
            coefficients_derivs_xmin.domain(),
            KOKKOS_LAMBDA(ddc::DiscreteElement<bsplines_type> i) {
                coefficients_derivs_xmin(i) *= ddc::detail::
                        ipow(dx_proxy,
                             static_cast<std::size_t>(get<bsplines_type>(
                                     (i - coefficients_derivs_xmin.domain().front()) + odd_proxy)));
            });
    ddc::parallel_for_each(
            exec_space(),
            coefficients_derivs_xmax.domain(),
            KOKKOS_LAMBDA(ddc::DiscreteElement<bsplines_type> i) {
                coefficients_derivs_xmax(i) *= ddc::detail::
                        ipow(dx_proxy,
                             static_cast<std::size_t>(get<bsplines_type>(
                                     (i - coefficients_derivs_xmax.domain().front()) + odd_proxy)));
            });
 
    ddc::DiscreteElement<deriv_type> const first_deriv(s_odd);
    // Allocate Chunk on deriv_type and interpolation_discrete_dimension_type and copy quadrature coefficients into it
    ddc::Chunk coefficients_derivs_xmin_out(
            ddc::DiscreteDomain<
                    deriv_type>(first_deriv, ddc::DiscreteVector<deriv_type>(s_nbv_xmin)),
            ddc::KokkosAllocator<double, OutMemorySpace>());
    ddc::Chunk coefficients_out(
            interpolation_domain().take_first(
                    ddc::DiscreteVector<interpolation_discrete_dimension_type>(
                            coefficients.size())),
            ddc::KokkosAllocator<double, OutMemorySpace>());
    ddc::Chunk coefficients_derivs_xmax_out(
            ddc::DiscreteDomain<
                    deriv_type>(first_deriv, ddc::DiscreteVector<deriv_type>(s_nbv_xmax)),
            ddc::KokkosAllocator<double, OutMemorySpace>());
    Kokkos::deep_copy(
            coefficients_derivs_xmin_out.allocation_kokkos_view(),
            coefficients_derivs_xmin.allocation_kokkos_view());
    Kokkos::deep_copy(
            coefficients_out.allocation_kokkos_view(),
            coefficients.allocation_kokkos_view());
    Kokkos::deep_copy(
            coefficients_derivs_xmax_out.allocation_kokkos_view(),
            coefficients_derivs_xmax.allocation_kokkos_view());
    return std::make_tuple(
            std::move(coefficients_derivs_xmin_out),
            std::move(coefficients_out),
            std::move(coefficients_derivs_xmax_out));
}
 
template <
        class ExecSpace,
        class MemorySpace,
        class BSplines,
        class InterpolationDDim,
        ddc::SplineBuilderClosure SBCLower,
        ddc::SplineBuilderClosure SBCUpper,
        SplineSolver Solver>
template <class KnotElement>
void SplineBuilder<
        ExecSpace,
        MemorySpace,
        BSplines,
        InterpolationDDim,
        SBCLower,
        SBCUpper,
        Solver>::
        check_n_points_in_cell(int const n_points_in_cell, KnotElement const current_cell_end_idx)
{
    if (n_points_in_cell > BSplines::degree() + 1) {
        KnotElement const rmin_idx = ddc::discrete_space<BSplines>().break_point_domain().front();
        int const failed_cell = (current_cell_end_idx - rmin_idx).value();
        throw std::runtime_error(
                "The spline problem is overconstrained. There are "
                + std::to_string(n_points_in_cell) + " points in the " + std::to_string(failed_cell)
                + "-th cell.");
    }
}
 
template <
        class ExecSpace,
        class MemorySpace,
        class BSplines,
        class InterpolationDDim,
        ddc::SplineBuilderClosure SBCLower,
        ddc::SplineBuilderClosure SBCUpper,
        SplineSolver Solver>
void SplineBuilder<
        ExecSpace,
        MemorySpace,
        BSplines,
        InterpolationDDim,
        SBCLower,
        SBCUpper,
        Solver>::check_valid_grid()
{
    std::size_t const n_interp_points = interpolation_domain().size();
    std::size_t const expected_npoints
            = ddc::discrete_space<BSplines>().nbasis() - s_nbe_xmin - s_nbe_xmax;
    if (n_interp_points != expected_npoints) {
        throw std::runtime_error(
                "Incorrect number of points supplied to NonUniformInterpolationPoints. "
                "(Received : "
                + std::to_string(n_interp_points)
                + ", expected : " + std::to_string(expected_npoints));
    }
    int n_points_in_cell = 0;
    auto current_cell_end_idx = ddc::discrete_space<BSplines>().break_point_domain().front() + 1;
    ddc::host_for_each(interpolation_domain(), [&](auto idx) {
        ddc::Coordinate<continuous_dimension_type> const point = ddc::coordinate(idx);
        if (point > ddc::coordinate(current_cell_end_idx)) {
            // Check the points found in the previous cell
            check_n_points_in_cell(n_points_in_cell, current_cell_end_idx);
            // Initialise the number of points in the subsequent cell, including the new point
            n_points_in_cell = 1;
            // Move to the next cell
            current_cell_end_idx += 1;
        } else if (point == ddc::coordinate(current_cell_end_idx)) {
            // Check the points found in the previous cell including the point on the boundary
            check_n_points_in_cell(n_points_in_cell + 1, current_cell_end_idx);
            // Initialise the number of points in the subsequent cell, including the point on the boundary
            n_points_in_cell = 1;
            // Move to the next cell
            current_cell_end_idx += 1;
        } else {
            // Indicate that the point is in the cell
            n_points_in_cell += 1;
        }
    });
    // Check the number of points in the final cell
    check_n_points_in_cell(n_points_in_cell, current_cell_end_idx);
}
 
} // namespace ddc