Tiled Stencils and Shared Memory Halos

A kernel that applies a five-point stencil to every element of an image reads five values per output and writes one. If each thread loads all five from global memory, the interior of the image is fetched five times across the grid — once by itself and once by each of four neighbours. The memory system will absorb much of that in the L2, but not all of it, and the traffic ratio of five loads per output is the first thing a tiling scheme is designed to fix. The second thing it fixes is the latency structure: a tile loaded once into shared memory can be re-read many times at the latency and bandwidth of shared memory rather than global. Wider stencils — a nine-point box, or a 25-point approximation of a Laplacian — make the same argument more strongly, because the ratio grows with the tap count while the halo grows only with the radius.

The tile is the work, the halo is the boundary

A block is assigned a rectangular region of the output, say $T_x \times T_y$ elements. To compute that region, the block needs the input region enlarged by the stencil radius on each side: $T_x + 2r$ by $T_y + 2r$ elements. The extra $2r$ rows and columns are the halo, and they are loaded by the same block that owns the interior, from global memory, into the shared-memory tile. No block communicates with another; the halo is fetched redundantly by each block that needs it.

For a 5-point stencil ($r = 1$) and a 32×8 tile, the shared array is 34×10 = 340 floats for 256 outputs, and the halo is 84 of those 340, a quarter of the loaded data. For a 32×32 tile the array is 34×34 = 1156 for 1024 outputs and the halo is 132, just under an eighth. That ratio is the whole argument for larger tiles, and it is why 2-D tiling has such a large effect on stencil throughput: global loads per output are $(T_x+2)(T_y+2)/(T_x T_y) = 1 + 2/T_x + 2/T_y + 4/(T_x T_y)$, so a 32×32 tile gives about 1.13 loads per output against 5 for the naive loop, which reads each of the five taps separately.

The costs that scale the other way are occupied shared memory per block — which limits how many blocks fit on a multiprocessor — and the barrier that separates the load from the compute. That barrier is the same block-wide one that cannot be reached by only part of the block, as described in block barriers and divergent synchronization. The tile cannot be larger than the shared memory available per block, and the limit is typically tens of kilobytes, so the tile size is chosen to fill that budget without stranding the block’s occupancy.

Loading the tile with a flat thread map

The tile has more elements than the block has threads, so the load is a loop, and the natural index is linear over the tile array. Flattening is what makes the load coalesced: consecutive threads take consecutive positions in the shared array, and if the shared array is laid out row-major, consecutive positions in one row are also consecutive in global memory.

#define TILE_X 32
#define TILE_Y 8
#define RADIUS 1

// Launch with blockDim == (TILE_X, TILE_Y); the compute phase assumes it.
__global__ void stencil5(const float* in, float* out, int w, int h) {
    __shared__ float tile[TILE_Y + 2][TILE_X + 2];
    const int ox = blockIdx.x * TILE_X;
    const int oy = blockIdx.y * TILE_Y;
    const int n = (TILE_X + 2) * (TILE_Y + 2);

    // One flat loop for the padded tile; each thread may store several cells.
    for (int t = threadIdx.y * blockDim.x + threadIdx.x; t < n;
         t += blockDim.x * blockDim.y) {
        const int sx = t % (TILE_X + 2);
        const int sy = t / (TILE_X + 2);
        // Clamp instead of testing: the halo edge must still hold a value.
        const int gx = min(max(ox + sx - RADIUS, 0), w - 1);
        const int gy = min(max(oy + sy - RADIUS, 0), h - 1);
        tile[sy][sx] = in[gy * w + gx];
    }
    __syncthreads();

    const int x = threadIdx.x;
    const int y = threadIdx.y;
    const int gx = ox + x;
    const int gy = oy + y;
    if (gx < w && gy < h) {
        // tile[y][x] is global (oy+y-1, ox+x-1): the halo shifts the interior by one.
        float s = tile[y + 1][x + 1] * 4.0f;
        s -= tile[y][x + 1] + tile[y + 2][x + 1]
           + tile[y + 1][x] + tile[y + 1][x + 2];
        out[gy * w + gx] = s;
    }
}

Three details in that listing are load-bearing.

The modulo and division by TILE_X + 2 compute a 2-D position in the padded tile. Because the padded width is 34 rather than a power of two, the division is an integer divide; the alternative is to iterate over rows and columns separately, which is faster to compute and slightly harder to make uniform across the block. This is the sort of overhead that disappears once the compute per element is more than a few operations.

The clamp on the global coordinate is the boundary condition. Clamping gives a zero-flux boundary; a constant-fill boundary would test the coordinate against the buffer and substitute a literal. Either works, but the load must not be skipped, because the barrier requires every thread in the block to reach it and a returned thread never arrives.

The write is guarded but the load is not. That asymmetry is deliberate: the load is what the barrier blocks on, and the store is not. The write guard also handles the last row and column of the image, where the tile extends past the output and threads compute values that nothing consumes. The fill-and-guard pattern is the only form that survives a tile extending past the buffer, for the reasons set out in tail handling and bounds checks.

Why the interior offsets need no checks

Once the tile is in shared memory, the compute phase indexes tile[y][x+1] and tile[y+2][x+1] with no bounds test. The reason is that the padding makes those indices legal for every thread that reaches them: a thread at y = 0 reads row 0, which is the halo row holding the image row above the tile, and a thread at the block’s last y reads row TILE_Y + 1, the halo row below. The array is declared with TILE_Y + 2 rows, so both exist. The +1 shift is the only place the halo layout is visible, and getting it wrong — indexing the interior as if the tile began at tile[0][0] — produces an output displaced by one element in both axes, which looks like a blur rather than a bug.

The corner case is what makes the shift worth stating. A thread at y = 0, x = 0 reads tile[0][0], tile[0][1], tile[1][0], tile[2][0] and tile[0][2] — all inside the array, all holding the correct clamped values because the halo was filled by the same clamp. Nothing underflows. An implementation that stored the interior without a halo and then subtracted one from the row index would read outside the array at the first row of the block, and whether that is a fault or silent garbage depends on what the shared-memory allocator placed there.

The same reasoning explains why the halo must be padded on both axes rather than just one, and why a 1-D shared array with a hand-computed pitch works only if the pitch includes the pad. A pitch of TILE_X with a halo in the row direction and none in the column direction leaves tile[y][-1] reading the previous row’s last element, which is a valid address holding the wrong value — the linear-index hazard in its tile-local form.

Shared-memory bank conflicts in the tile

A shared-memory array is banked by address. On the common architectures there are 32 banks of 4 bytes each, addressed modulo the bank count. A 2-D array with a row length that is a multiple of the bank count has a specific pathology: the column-major accesses tile[y+k][x] for fixed x and varying y all land in the same bank, giving a 32-way conflict and serialising the read. With a row of 32 floats, column access collides; with a row of 33 floats, it does not.

The tile above is 34 floats wide, so stepping one row advances the bank index by 2. Thirty-two consecutive rows therefore land on 16 distinct banks, each hit by two threads: a 2-way conflict, tolerable. A tile of width 32 would advance by 0 and put all 32 rows on one bank, a 32-way conflict. The universal fix is to add a pad of one element to the row stride, which is why shared arrays are so often declared one wider than their logical width — a power-of-two width is the worst case, and a width one greater than a power of two is the best.

Register blocking on top of the tile

The shared-memory tile fixes global traffic. The next layer fixes shared-memory traffic, which becomes the limit once the tile is large: each output reads four neighbours and itself, and adjacent outputs read overlapping sets. A thread that computes a $2 \times 2$ or $4 \times 4$ patch of outputs can load the surrounding $(p+2) \times (p+2)$ values into registers once and compute $p^2$ outputs from them, cutting shared-memory reads per output from five to roughly $(p+2)^2 / p^2$.

For $p = 4$ that is 36 reads for 16 outputs, or 2.25 per output, against 5. The cost is registers: the patch plus the accumulator is 36 values held live, which is a fifth of a 192-register budget and pushes the block’s occupancy down. The occupancy loss must be weighed against the reduced shared-memory traffic, and the trade favours register blocking when the block is already limited by something other than registers.

Register blocking also fixes the number of threads relative to the tile. A thread mapping of 32×8 with $p = 4$ means the block’s threads correspond to a $32 \times 4$ grid of patches covering a $128 \times 32$ output region, which requires a shared tile of $130 \times 34$ floats, about 17 KB. That fits, but the tile size now exceeds the block size by a factor of 16, so the load loop runs 16 iterations per thread. The load is more work per thread, and the halo is a smaller fraction of the tile — both good — but the barrier now separates a long load from a short compute, so the block’s warps are more likely to be unevenly paced at the barrier.

Double buffering and the serialisation of load and compute

A single-tile kernel alternates: all warps load, all warps wait, all warps compute. Nothing overlaps. On a machine that hides latency by having several blocks resident per multiprocessor, the idle time during a block’s barrier is partly filled by another block’s compute, which is the main reason small blocks with modest tiles still achieve good throughput.

The explicit alternative is double buffering: two tiles in shared memory, with the block loading tile $k+1$ while computing on tile $k$. This requires two barriers per iteration (one after the load is visible, one after the compute has finished reading the buffer about to be overwritten) and doubles the shared-memory footprint. It is worth it when the problem is one long strip being processed in stages, where there is no other block to fill the gap; it is usually not worth it for a single-pass image transform where the grid already provides independent work.

The rule of thumb is that any mechanism whose purpose is to overlap a block’s own load and compute is competing with the hardware’s inter-block overlap, and the hardware’s version is free. Adding the software version costs shared memory and barriers, so it should be adopted only after measuring that the blocks are not filling the machine.