Order 1 splines (Id matrix) is not handled like it is in SLL (avoid solver call)

[blegouix]

    if constexpr (bsplines_type::degree() == 1)
        return compute_interpolant_degree1(spline, vals);

Is not in the main

[blegouix]

With (copy-paste from SLL, it was non-batched splines):

template <class BSplines, class interpolation_mesh_type, BoundCond BcXmin, BoundCond BcXmax>
void SplineBuilder<BSplines, interpolation_mesh_type, BcXmin, BcXmax>::compute_interpolant_degree1(
        ddc::ChunkSpan<double, ddc::DiscreteDomain<bsplines_type>> const spline,
        ddc::ChunkSpan<double const, interpolation_domain_type> const vals) const
{
    for (std::size_t i = 0; i < ddc::discrete_space<BSplines>().nbasis(); ++i) {
        spline(ddc::DiscreteElement<bsplines_type>(i))
                = vals(ddc::DiscreteElement<interpolation_mesh_type>(i));
    }
    if constexpr (bsplines_type::is_periodic()) {
        spline(ddc::DiscreteElement<bsplines_type>(ddc::discrete_space<BSplines>().nbasis()))
                = spline(ddc::DiscreteElement<bsplines_type>(0));
    }
}

[blegouix]

Atm I think we build the Identity matrix and solve IX=B, which is pretty stupid. A natural way to handle this would be to build a SplinesLinearProblem (the parent class) or a SplinesLinearProblemIdentity, whose factorize and solve methods do nothing, and call a make_new_identity for degree==1 case.

The downside of it is it goes through the GPU kernel https://github.com/CExA-project/ddc/blob/ce2108ef653420f235eef5f574545a231dd2ec1d/include/ddc/kernels/splines/spline_builder.hpp#L776 which do nothing because nbc_xmin=0, and more important it transposes and transposes back at https://github.com/CExA-project/ddc/blob/ce2108ef653420f235eef5f574545a231dd2ec1d/include/ddc/kernels/splines/spline_builder.hpp#L828

Alternatively we can write the code to treat the special case and bypass everything, like it was in SLL. But it would require testing the degree=1 case so I don't think this is the simplest atm.

[blegouix]

@tpadioleau should I address this, and if yes do you know in which way ?

[tpadioleau]

Quite low priority for now, unless @EmilyBourne thinks it could be useful ?

@blegouix Can you test whether it actually works if we use this class with degree 1 splines ?

[blegouix]

It's low priority but as I am leaving in three weeks, either I address it either I unassign myself ^^

[tpadioleau]

I would say priority to check that it works correctly.

[blegouix]

I works with GINKGO, not with KK because it calls getrs with a matrix of size 0x0

Using BiCGStab to solve IX=Y is 🤘

[tpadioleau]

I works with GINKGO, not with KK because it calls getrs with a matrix of size 0x0

Can you make it work ?

[blegouix]

This commit make it work: https://github.com/CExA-project/ddc/pull/480/commits/beb1ee90e0ada1171cd299fd1705441936caf9ff

In the main it is already working (I guess the LAPACK call is fine with matrix of size 0)

[tpadioleau]

Can you identify the KK functions that do not work with size 0 ? With what kind of error ?

[blegouix]

KokkosBatched::getrs: the leading dimension of the array A must satisfy lda >= max(1, n): A: 0 x 0

Actually this is just a warning which floods the cout, the simulation runs correctly: