Batches and layouts

Three systems with the same solution, laid out as each snippet says. The batch is filled above the shown code.

\[\begin{split} A^{(s)} = A + s\,I, \quad b^{(s)} = b + s\,x, \quad s = 0, 1, 2, \quad A = \left(\begin{array}{cc|cc} 4 & 1 & 0 & 2 \\ 1 & 5 & 1 & 0 \\ \hline 0 & 1 & 6 & 1 \\ 2 & 0 & 1 & 7 \end{array}\right), \quad b = \begin{pmatrix} 14 \\ 14 \\ 24 \\ 33 \end{pmatrix}, \quad x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \end{pmatrix} \end{split}\]

Residencies

The same solve on local arrays, on a slice of an SoA batch, and with the matrix in the batch while the pivot and the right-hand side stay local. Element k of system s sits at:

\[ \text{index}(s, k) = 3\,k + s \]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

int piv[4];

// internal: every operand is a caller-local array, the strides are ignored
const bool ok_local = Solver::solve_inplace<true, true, true>(A_local, 1, piv, 1, y_local, 1);

// external: every operand is system 1 of the SoA batch, walked with the batch stride 3
const bool ok_batch =
    Solver::solve_inplace<false, false, false>(A_soa + 1, 3, piv_soa + 1, 3, y_soa + 1, 3);

// mixed: the matrix of system 2 stays in the batch, its pivot and right-hand side are local
const bool ok_mixed =
    Solver::solve_inplace<true, true, false>(A_soa + 2, 3, piv, 1, y_mixed, 1);

AoS batch

Array of structures: the objects of a system are contiguous, the stride is 1. Element k of system s sits at:

\[ \text{index}(s, k) = 16\,s + k \]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverDynamic<double, config>;

// the objects of system s are contiguous, at A + s * 16 and y + s * 4
for (int s = 0; s < count; ++s) {
    int piv[4];
    if (!Solver::solve_inplace(4, A + s * 16, 1, piv, 1, y + s * 4, 1)) return 1;
}

SoA batch

Structure of arrays: consecutive systems touch consecutive addresses, the stride is the batch size. Element k of system s sits at:

\[ \text{index}(s, k) = 3\,k + s \]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverDynamic<double, config>;

// element k of system s at A + k * 3 + s: the stride is the batch size
for (int s = 0; s < count; ++s) {
    int piv[4];
    if (!Solver::solve_inplace(4, A + s, 3, piv, 1, y + s, 3)) return 1;
}

AoSoA batch

Array of structures of arrays: blocks of W = 2 interleaved systems, the stride is W. The batch is padded to 4 systems so that every block is full. Element k of system s sits at:

\[ \text{index}(s, k) = 32 \left\lfloor s / 2 \right\rfloor + 2\,k + (s \bmod 2) \]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverDynamic<double, config>;

// block s / 2, slot s % 2: the stride is the width 2
for (int s = 0; s < count; ++s) {
    int piv[4];
    double* A_s = A + (s / 2) * 16 * 2 + s % 2;
    double* y_s = y + (s / 2) * 4 * 2 + s % 2;
    if (!Solver::solve_inplace(4, A_s, 2, piv, 1, y_s, 2)) return 1;
}

An OpenMP loop over a batch

One iteration per system of an SoA batch, in a function shown whole. Element k of system s sits at:

\[ \text{index}(s, k) = \text{count} \cdot k + s \]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

// one iteration per system of the SoA batch
void solve_batch(const int count, double* A, double* y, int* ok) {
#pragma omp parallel for
    for (int s = 0; s < count; ++s) {
        // matrix and right-hand side in the batch, pivot local to the iteration
        int piv[4];
        ok[s] =
            Solver::solve_inplace<false, true, false>(A + s, count, piv, 1, y + s, count) ? 1 : 0;
    }
}

Inside a GPU kernel

One thread per system of an SoA batch, in the common CUDA/HIP dialect: the kernel is shown, the host side allocates, launches and copies back. Element k of system s sits at:

\[ \text{index}(s, k) = \text{count} \cdot k + s \]

using Solver =
    tdls::TiledLUppSolverStatic<double, 4, tdls::TiledLUppConfig<double>{.tile_size = 2}>;

// one thread per system of the SoA batch
__global__ void solve_batch(const int count, double* A, double* y, int* ok) {
    const int s = static_cast<int>(blockIdx.x) * blockDim.x + threadIdx.x;
    if (s >= count) return;

    // matrix and right-hand side in the batch, pivot in registers
    int piv[4];
    ok[s] = Solver::solve_inplace<false, true, false>(A + s, count, piv, 1, y + s, count) ? 1 : 0;
}