Implement batched serial tbsv

[yasahi-hpc]

This is a draft PR to discuss the implementation strategy of tbsv function.

If I understand correctly, a minimal set for each kernel includes

1. KokkosBatched_Tbsv_Serial_Impl.hpp: Internal interfaces
2. KokkosBatched_Tbsv_Serial_Internal.hpp: Implementation details
3. KokkosBatched_Tbsv_Decl.hpp: APIs
4. Test_Batched_SerialTbsv.hpp: Unit tests for that

Detailed description

It solves the equation Ax = b.
Here, the matrix has the following shape.

• A: (batch_count, lda, n)
n by n unit or non-unit, upper or lower triangular band matrix with (k+1) diagonals.
• x: (batch_count, n)
  Before entry, the incremented array x must contain the n element right-hand side vector b.

Example of a single batch of matrix A with n = 10 and k = 3.

upper
1 1 1 1 0 0 0 0 0 0 
0 1 1 1 1 0 0 0 0 0 
0 0 1 1 1 1 0 0 0 0 
0 0 0 1 1 1 1 0 0 0 
0 0 0 0 1 1 1 1 0 0 
0 0 0 0 0 1 1 1 1 0 
0 0 0 0 0 0 1 1 1 1 
0 0 0 0 0 0 0 1 1 1 
0 0 0 0 0 0 0 0 1 1 
0 0 0 0 0 0 0 0 0 1 

lower
1 0 0 0 0 0 0 0 0 0 
1 1 0 0 0 0 0 0 0 0 
1 1 1 0 0 0 0 0 0 0 
1 1 1 1 0 0 0 0 0 0 
0 1 1 1 1 0 0 0 0 0 
0 0 1 1 1 1 0 0 0 0 
0 0 0 1 1 1 1 0 0 0 
0 0 0 0 1 1 1 1 0 0 
0 0 0 0 0 1 1 1 1 0 
0 0 0 0 0 0 1 1 1 1 

upper_banded
0 0 0 1 1 1 1 1 1 1 
0 0 1 1 1 1 1 1 1 1 
0 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 

lower_banded
1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 0 
1 1 1 1 1 1 1 1 0 0 
1 1 1 1 1 1 1 0 0 0 

Parallelization would be made in the following manner. This is efficient only when
A is given in LayoutLeft for GPUs and LayoutRight for CPUs (parallelized over batch direction).

Kokkos::parallel_for('tbsv', 
    Kokkos::RangePolicy<execution_space, ParamTagType> policy(0, n),
    [=](const ParamTagType &, const int k) {
        auto aa = Kokkos::subview(_a, k, Kokkos::ALL(), Kokkos::ALL());
        auto bb = Kokkos::subview(_b, k, Kokkos::ALL());

        KokkosBatched::SerialTbsv<
            typename ParamTagType::uplo, typename ParamTagType::trans,
            typename ParamTagType::diag, AlgoTagType>::invoke(aa, bb, k, incx);
    });

Tests

Firstly, construct an Upper or Lower banded diagonal matrix A. Then, convert it to the banded storage Ab. Solving Ax = b with trsv and (Ab)x = b with tbsv. tbsv result is compared with trsv result.

[yasahi-hpc]

For the moment, a parameter incx is ignored for simplicity

[yasahi-hpc]

TODO tasks

1. Development

• Public interface should be KokkosBatched_Tbsv.hpp not KokkosBatched_Tbsv_Decl.hpp. This is a recent trend.
• Check the matrix format. It should be a rank-2 matrix (already confirmed) with the appropriate storage scheme as described in banded storage format.
• Enable incx parameter

2. Tests

• Add a simple and small analytical test, i.e. chose A and b and compute by hand x, then check that it matches with your routine.
• Do some corner cases, what happens when you have an empty matrix? Zero upper/lower diagonals, zeros on the diagonal, etc...
• Generate random matrices for A and x, compute b with a matvec and then check that solving Ax=b returns the original x vector. (This is partially done. Tbsv is compared with Trsv)
• Ideally if you have any assertions/runtime checks that can be done, you want to check that they are catching what they are supposed to catch
• Finally you want the test to cover all backends for which the algorithms is supposed to run. (Serial, OpenMP, Thread, CUDA, HIP can be tested)

[yasahi-hpc]

Most of the modifications are made.