3D Stencil: Rotate the Depth Window Through Registers

Lesson

Start from where you left off. This lab's starter IS the previous lab's solution, on the same harness, the same fixture and the same profile size, so its measured cost, about 977 million sectors, is the number you already earned. The question is where the remaining traffic hides.

Follow one plane of the cube through the starter's sweep. The block fetches it from global-memory request to stage the tile. One step later the block fetches it again, as the behind neighbour. One step earlier it had already fetched it, as the ahead neighbour. Every plane of the cube crosses global-memory request three times, and two of those three crossings are pure re-reads of data the block held moments ago.

The fix is the row-window idea from the 2D register lab, applied to the axis the block already walks in order. Keep three per-thread values: behind, current, ahead. Each step reads only the plane entering the window, uses all three, then rotates them: current becomes behind, ahead becomes current. The shared tile does not go away, it is simply stored from the current register instead of re-read from global-memory request, and the in-plane neighbours still cross it.

Where each value lives, per z step:
Starter: the tile staged from global-memory request; behind read from global-memory request; ahead read from global-memory request.
Solution: behind held in a register; the tile stored from the current register; ahead read from global-memory request, the only read.
Then the window rotates: ahead becomes current, current becomes behind.

Every input element is now fetched from global-memory request once, which is the floor for any kernel that reads its input at all. This is 2.5D blocking, the same structure production 3D stencil kernels use, and note what it did not cost: shared memory did not grow by a byte, and the barrier count per step did not change.

Two honest notes. A rolling three-plane shared buffer matches this traffic and also passes; the register form is what scales, because its footprint is independent of how deep the block sweeps. And a correct submission may trip the register-pressure insight, which suggests capping registers with launch bounds. Here that advice points the wrong way: the extra registers are the optimization. Measured, the rotation drops the total from about 977 million sectors to about 567 million, roughly 1.7x, and the gate at 660 million sits between them.

Metric contract: algorithmic reads/writes count source-level sweeps. GLOBAL_MEMORY_SECTORS counts SASS-derived issued/requested sectors in the selected profiling scope, including repeated launches. It is cache-blind: it is neither L2 traffic nor device-memory bytes. Trace elapsed time is a separate measurement; a global load may hit cache.

Register rotation reuses interior depth values. Halo overlap still duplicates requests between tiles; it does not guarantee exactly one device-memory fetch per plane.

What you will learn

• Recognize that a depth neighbor in one iteration is the working plane of another, and count global traffic per element rather than per thread.
• Carry a rolling three-plane window in registers and rotate it once per z step so each plane is fetched exactly once.
• Fill a shared tile interior from a register while the halo ring keeps coming from global memory.
• Judge a register-blocking change by issued global-memory sectors, and recognize when generic register-pressure advice does not apply.

Starter code

    // The four in-plane neighbours overlap between threads of this block, so
    // stage one 34 x 18 plane per z step and read them from shared memory. The
    // two depth neighbours belong to other planes and are still read globally,
    // which is why this only removes two of the six neighbour reads.
    __shared__ float tile[18][34];

    const int x = blockIdx.x * blockDim.x + threadIdx.x;
    const int y = blockIdx.y * blockDim.y + threadIdx.y;
    const int tx = threadIdx.x + 1;
    const int ty = threadIdx.y + 1;

    // Clamp rather than return early: every launched lane has to reach both
    // barriers below, so an out-of-range lane stages a duplicated edge value.
    const int cx = x < N ? x : N - 1;
    const int cy = y < N ? y : N - 1;
    const int xl = x > 0 ? x - 1 : 0;
    const int xr = (x + 1) < N ? x + 1 : N - 1;
    const int yu = y > 0 ? y - 1 : 0;
    const int yd = (y + 1) < N ? y + 1 : N - 1;

    const int planesPerBlock = (N + gridDim.z - 1) / gridDim.z;
    const int z0 = blockIdx.z * planesPerBlock;
    if (z0 >= N) {
        return;
    }
    const int zEnd = (z0 + planesPerBlock) < N ? (z0 + planesPerBlock) : N;
    const int plane = N * N;

    if (threadIdx.x == 0 && threadIdx.y == 0) {
        // Defined but never read: a seven-point rule has no use for the four
        // in-plane corners, so they cost four shared stores and no traffic.
        tile[0][0] = 0.0f;
        tile[0][33] = 0.0f;
        tile[17][0] = 0.0f;
        tile[17][33] = 0.0f;
    }

    for (int z = z0; z < zEnd; ++z) {
        const int base = z * plane;

        // Barrier first: the previous step's readers must finish before this
        // step overwrites the tile they are still reading.
        __syncthreads();
        tile[ty][tx] = input[base + cy * N + cx];
        if (threadIdx.x == 0) {
            tile[ty][0] = input[base + cy * N + xl];
        }
        if (threadIdx.x == blockDim.x - 1) {
            tile[ty][33] = input[base + cy * N + xr];
        }
        if (threadIdx.y == 0) {
            tile[0][tx] = input[base + yu * N + cx];
        }
        if (threadIdx.y == blockDim.y - 1) {
            tile[17][tx] = input[base + yd * N + cx];
        }
        __syncthreads();

        if (x >= N || y >= N) {
            continue;
        }
        const int o = base + y * N + x;
        const float c = tile[ty][tx];
        if (x == 0 || y == 0 || z == 0 || x == N - 1 || y == N - 1 || z == N - 1) {
            output[o] = c;
            continue;
        }
        output[o] = 0.25f * c
                + 0.125f * (tile[ty][tx - 1] + tile[ty][tx + 1]
                        + tile[ty - 1][tx] + tile[ty + 1][tx]
                        + input[o - plane] + input[o + plane]);
    }