Working With Multidimensional Kernels

Many important examples in parallel computing are naturally multidimensional: images, matrices, heat diffusion, cellular automata, and finite-difference stencils. alpaka supports multi-dimensional kernels and memory. That’s why it’s usually better to implement these types of problems using multiple dimensions rather than flattening everything into a single linear index.

Choose the Kernel Shape From the Data

If the data is naturally a matrix or image, it is convenient to use two-dimensional extents and frames. This avoids hand-written index decoding and makes boundary conditions easier to read.

onHost::concepts::FrameSpec auto frameSpec = onHost::getFrameSpec(device, exec::anyExecutor, problemExtents);

We recommend that the frame shape should follow the logical shape of the work. Typical use cases for different dimensions are:

  • 1D frames for flat vectors and simple reductions.

  • 2D frames for images, matrices, and most stencil codes.

  • 3D frames for volumetric problems, such as the position of particles in a space.

  • ND frames for algorithms that have more than three natural dimensions, such as position coordinates and discretized time, e.g., to build a layered detector within our particle simulation.

Keep in mind that the rightmost index, usually x, is the fastest varying dimension in alpaka buffers. This is important if you are using hand written for-loops to iterate over the data instead of makeIdxMap.

A Small 2D Stencil Example

The following kernel performs one five-point average step on a small 2D grid [1]. This introduces three important ideas at once:

  • Iterating over multidimensional buffers.

  • Handling boundaries explicitly.

  • Reading neighboring cells by using relative offsets.

    struct FivePointAverageKernel
{
    ALPAKA_FN_ACC void operator()(
        onAcc::concepts::Acc auto const& acc,
        concepts::IMdSpan auto out,
        concepts::IDataSource auto const& in) const
    {
        auto extents = out.getExtents();
        ALPAKA_ASSERT_ACC(extents == in.getExtents());
        constexpr auto xDir = CVec<uint32_t, 0u, 1u>{};
        constexpr auto yDir = CVec<uint32_t, 1u, 0u>{};

        for(auto idx : onAcc::makeIdxMap(acc, onAcc::worker::threadsInGrid, IdxRange{extents}))
        {
            /* Copy only values in the halo region
            idx.y() == 0u -> top halo, idx.x() == 0u -> left halo
            idx.y() + 1u == extents.y() -> bottom halo, idx.x() + 1u == extents.x() -> right halo */
            if(idx.y() == 0u || idx.x() == 0u || idx.y() + 1u == extents.y() || idx.x() + 1u == extents.x())
            {
                out[idx] = in[idx];
                continue;
            }

            out[idx] = (in[idx] + in[idx - yDir] + in[idx + yDir] + in[idx - xDir] + in[idx + xDir]) / 5;
        }
    }
};

The structure is still the same as in the one-dimensional tutorial:

  • Ask the output buffer for its extents.

  • Build IdxRange{extents} to describe the full valid multidimensional index domain.

  • Map worker threads via makeIdxMap to the range.

  • Guard the boundary cells, to avoid out-of-bounds accesses.

  • The own element is updated depending on the neighbor data and maybe the own element. Indices are relative calculated compare to the index of the own output coordinate.

  • Update neighbor locations by adding or subtracting direction vectors from the current Vec index.

This is the natural alpaka style for stencil code. The project examples, such as the heat-equation stencil, operate on the multidimensional index directly and move to neighbors with vector offsets instead of splitting x and y into separate scalars and rebuilding indices.

At the beginning, the range check of data access should be done within the kernel. While this approach does not offer the best performance, it is easier to implement correctly. Subsequently, the range check should be moved from the kernel to before the kernel launch, and the kernel should operate on the assumption that the range is correct in order to improve performance.

If you want to optimize the boundary handling later, you can write a separate kernel for the border and launch it with a different frame shape that only covers the halo/guard region.

Launching the 2D Kernel

The host-side launch is unchanged except that the problem extents are vectors now, and the chosen FrameSpec is 2-dimensional as well.

static_assert(frameSpec.dim() == 2);
queue.enqueue(frameSpec, KernelBundle{FivePointAverageKernel{}, outBuffer, inBuffer});

onHost::memcpy(queue, hostOutput, outBuffer);
onHost::wait(queue);

What Users Usually Need To Know Early

The following habits are worth learning from the start:

  • Keep boundary handling explicit. It is fine to have branching for the border in the first implementation.

  • Use multidimensional buffers when the algorithm is multidimensional.

  • Use temporary variables to avoid reading from and writing to memory multiple times. The variable name provides a meaningful name for read and write operations, making it easier to track when data is being read and written.

  • Start with a clear kernel and a small test case before trying to optimize shared memory use or tuning the FrameSpec

Complete Source File

080_multidim.cpp
 1/* Copyright 2026 René Widera
 2 * SPDX-License-Identifier: ISC
 3 */
 4
 5#include "docsTest.hpp"
 6
 7#include <alpaka/alpaka.hpp>
 8
 9#include <catch2/catch_template_test_macros.hpp>
10#include <catch2/catch_test_macros.hpp>
11
12using namespace alpaka;
13
14struct FivePointAverageKernel
15{
16    ALPAKA_FN_ACC void operator()(
17        onAcc::concepts::Acc auto const& acc,
18        concepts::IMdSpan auto out,
19        concepts::IDataSource auto const& in) const
20    {
21        auto extents = out.getExtents();
22        ALPAKA_ASSERT_ACC(extents == in.getExtents());
23        constexpr auto xDir = CVec<uint32_t, 0u, 1u>{};
24        constexpr auto yDir = CVec<uint32_t, 1u, 0u>{};
25
26        for(auto idx : onAcc::makeIdxMap(acc, onAcc::worker::threadsInGrid, IdxRange{extents}))
27        {
28            /* Copy only values in the halo region
29            idx.y() == 0u -> top halo, idx.x() == 0u -> left halo
30            idx.y() + 1u == extents.y() -> bottom halo, idx.x() + 1u == extents.x() -> right halo */
31            if(idx.y() == 0u || idx.x() == 0u || idx.y() + 1u == extents.y() || idx.x() + 1u == extents.x())
32            {
33                out[idx] = in[idx];
34                continue;
35            }
36
37            out[idx] = (in[idx] + in[idx - yDir] + in[idx + yDir] + in[idx - xDir] + in[idx + xDir]) / 5;
38        }
39    }
40};
41
42
43TEMPLATE_LIST_TEST_CASE("tutorial multidimensional stencil kernel", "[docs]", docs::test::TestBackends)
44{
45    auto selector = onHost::makeDeviceSelector(TestType::makeDict());
46    if(!selector.isAvailable())
47        return;
48    onHost::concepts::Device auto device = selector.makeDevice(0);
49    onHost::Queue queue = device.makeQueue();
50
51    auto const problemExtents = Vec{5u, 5u};
52    auto hostInput = onHost::allocHost<int>(problemExtents);
53    auto hostOutput = onHost::allocHostLike(hostInput);
54
55    for(auto& value : hostInput)
56        value = 0;
57    hostInput[Vec{2u, 2u}] = 100;
58
59    auto inBuffer = onHost::allocLike(device, hostInput);
60    auto outBuffer = onHost::allocLike(device, hostInput);
61
62    onHost::memcpy(queue, inBuffer, hostInput);
63    onHost::memset(queue, outBuffer, 0x00);
64
65    onHost::concepts::FrameSpec auto frameSpec = onHost::getFrameSpec(device, exec::anyExecutor, problemExtents);
66
67    static_assert(frameSpec.dim() == 2);
68    queue.enqueue(frameSpec, KernelBundle{FivePointAverageKernel{}, outBuffer, inBuffer});
69
70    onHost::memcpy(queue, hostOutput, outBuffer);
71    onHost::wait(queue);
72
73    for(auto const idx : IdxRange{problemExtents})
74    {
75        auto const isCross = idx == Vec{2u, 2u} || idx == Vec{1u, 2u} || idx == Vec{2u, 1u} || idx == Vec{2u, 3u}
76                             || idx == Vec{3u, 2u};
77        CHECK(hostOutput[idx] == (isCross ? 20 : 0));
78    }
79}