3D Stencil: Stage the Plane, Pay for the Depth

A seven-point stencil over an N x N x N cube reads each point's six face neighbours: four in the same z plane, one behind, one ahead. The starter fetches all seven values from global-memory request per output; the fix stages each z plane in a shared tile as the block sweeps, which serves five of the seven reads on-chip and honestly cannot serve the other two. GPUFlight judges issued global-memory sectors.

Languages: CUDA C++, Rust

Lesson

The computation: for every interior point of the cube, output = 0.25 x centre + 0.125 x (left + right + up + down + behind + ahead), where behind and ahead are the same (x, y) position one plane away in z. Any point on a face of the cube is copied through unchanged. The profile case is a 256-cube, 16.8 million points, the same element count as the 2D labs, so the sector totals across this whole family sit on one axis.

A 32 x 16 block owns one column of the cube and sweeps it plane by plane. The harness launches a fixed 3 z-chunks, so the sweep depth depends on N and must be derived as planesPerBlock = (N + gridDim.z - 1) / gridDim.z; a hardcoded depth leaves planes unwritten and fails against the poisoned output.

In three dimensions you cannot stage the whole neighbourhood, so you stage what you can: each step, the block loads the current plane into a static 18 x 34 shared tile. The four in-plane neighbours now come from shared memory. The two depth neighbours cannot: they live in different planes, so they stay full global-memory request reads.

Where each read is served from:
Starter, per output: centre and all six neighbours from global-memory request, seven reads.
Solution, per output: centre, left, right, up, down from the shared tile; behind global-memory request; ahead global-memory request.

That is the honest accounting, and it is the lesson: staging the plane removes five of seven reads but only about 1.8x of the traffic, not the 4x the shape of the code suggests, because the two depth reads are whole sweeps of the cube and they are untouched. The disappointment is deliberate. It is the setup for the next lab, which removes exactly those two reads.

Why a 32-wide tile and not a cube-shaped one: an 8 x 8 x 8 brick stages about 1,000 cells to produce 512 outputs, and its 10-float rows waste three quarters of every 32-byte sector they touch. A 32 x 16 plane tile stages 612 cells for 512 outputs with whole-sector rows. Fewer staged cells per output AND fewer sectors per staged cell, which is why the harness pins a full-warp fast axis: 32 x 16 x 1, blockZ 1, because z is walked by the loop, not by threads.

The tile is static, so sharedMemBytes stays 0, and the harness derives the 3-D grid itself, so the launch panel's Grid fields are inert. Measured: the starter issues about 1.77 billion sectors, the plane-tiled version about 977 million, and the gate at 1.15 billion 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.

Plane staging reuses x/y neighbors but still requests depth neighbors. Check all six boundary faces; cache hits determine how many requests reach device memory.

What you will learn

Stage one z plane of a 3D stencil in a shared tile with a one cell halo on all four in-plane sides.
Derive the sweep depth from gridDim.z instead of hardcoding it, so the kernel holds at every N the harness tests.
Place two barriers per z step so a reused tile is never overwritten while it is still being read.
Account honestly for what tiling did NOT remove, and read the measured sector drop rather than the expected one.

Starter code

    // Seven-point stencil over a cube. This block owns a 32 x 16 column and one
    // chunk of z planes, and reads every neighbour from global memory, including
    // the four in-plane neighbours its own block already holds.
    const int x = blockIdx.x * blockDim.x + threadIdx.x;
    const int y = blockIdx.y * blockDim.y + threadIdx.y;

    // Derive the chunk depth from the grid the harness launched, so the block
    // covers exactly its own planes for every N.
    const int planesPerBlock = (N + gridDim.z - 1) / gridDim.z;
    const int z0 = blockIdx.z * planesPerBlock;
    if (z0 >= N || x >= N || y >= N) {
        return;
    }
    const int zEnd = (z0 + planesPerBlock) < N ? (z0 + planesPerBlock) : N;
    const int plane = N * N;

    for (int z = z0; z < zEnd; ++z) {
        const int o = (z * N + y) * N + x;
        if (x == 0 || y == 0 || z == 0 || x == N - 1 || y == N - 1 || z == N - 1) {
            output[o] = input[o];
            continue;
        }
        output[o] = 0.25f * input[o]
                + 0.125f * (input[o - 1] + input[o + 1]
                        + input[o - N] + input[o + N]
                        + input[o - plane] + input[o + plane]);
    }