2D Stencil: Roll the Row Window Through Registers

Lesson

The computation is unchanged from the previous lab: output = 0.5 x centre + 0.125 x (left + right + up + down), borders copied through. What changes is who owns what. A 256-thread block owns one column-strip of the grid and walks it top to bottom, one output row per step.

The harness launches a fixed 16 row chunks, so the strip height depends on N: about 256 rows per block at the profile size, and odd heights at the tail sizes. Your kernel must derive it as rowsPerBlock = (N + gridDim.y - 1) / gridDim.y. A hardcoded height under-covers the grid, and the harness pre-poisons the output buffer, so the missing rows fail verification loudly.

Now the waste. At each step the starter reads row y-1, row y and row y+1 from global-memory request. One step later, two of those three are needed again, and it reads them again. Across the strip, every row crosses global-memory request three times: once as y+1, once as y, once as y-1.

The fix is to notice the window is the same data sliding by one: what was row y+1 this step IS row y next step. Values a thread will need again belong in its registers, not in repeated global loads.

Where each read is served from:
Starter, per step: row y-1 global-memory request, row y global-memory request, row y+1 global-memory request.
Solution, per step: row y-1 register, row y register, row y+1 global-memory request, the one new read.
Then the window rolls: y+1 becomes y, y becomes y-1.

Registers are private to a thread, so they can only carry the vertical neighbours, which live in the same column. The left and right neighbours live in other threads, and that is all shared memory still does here: one 258-float row, published each step, with a barrier on each side of the reuse. Shared memory shrank from a 2-D tile to a single row precisely because the y direction moved into registers.

Two honest notes. A rolling three-row shared buffer matches this traffic and also passes the gate; the register form is taught because its shared footprint does not grow with the strip height. And this lab's starter costs more than the previous lab's solution, about 853 million sectors against 578 million, because the strip contract re-reads rows before the fix lands; the path is not going backwards, it is setting up the fix. Measured, the solution issues about 435 million sectors, roughly 2x under the starter, and the gate at 530 million sits between them. The block is pinned at 256 x 1 x 1 and the shared row is static, so sharedMemBytes stays 0.

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.

Rolling rows trade fewer shared/global requests for register pressure. Compare registers, occupancy and elapsed time together.

What you will learn

• Recognize a rolling window: a stencil neighborhood that shifts by one row per step already holds two of its three values.
• Carry the y neighbors in registers across loop iterations and roll them forward at the end of each step.
• Keep shared memory to the neighbors that genuinely live in other lanes, which for this stencil is the x direction only.
• Derive the row chunk height from gridDim.y so the block covers exactly its own rows for every N the harness tests.

Starter code

    // Strip contract: this block owns blockDim.x columns and one chunk of rows,
    // and walks the chunk one row at a time. The x neighbours come from a shared
    // row; the y neighbours are re-read from global on every single step.
    __shared__ float row[258];

    const int x = blockIdx.x * blockDim.x + threadIdx.x;
    const int cx = x < N ? x : N - 1;
    const int sx = threadIdx.x + 1;
    const int xl = x > 0 ? x - 1 : 0;
    const int xr = (x + 1) < N ? x + 1 : N - 1;

    // Derive the chunk height from the grid the harness launched, so the block
    // covers exactly its own rows for every N instead of a hardcoded count.
    const int rowsPerBlock = (N + gridDim.y - 1) / gridDim.y;
    const int y0 = blockIdx.y * rowsPerBlock;
    if (y0 >= N) {
        return;
    }
    const int yEnd = (y0 + rowsPerBlock) < N ? (y0 + rowsPerBlock) : N;

    for (int y = y0; y < yEnd; ++y) {
        const float above = (y > 0) ? input[(y - 1) * N + cx] : 0.0f;
        const float centre = input[y * N + cx];
        const float below = ((y + 1) < N) ? input[(y + 1) * N + cx] : 0.0f;

        row[sx] = centre;
        if (threadIdx.x == 0) {
            row[0] = input[y * N + xl];
        }
        if (threadIdx.x == blockDim.x - 1) {
            row[257] = input[y * N + xr];
        }
        __syncthreads();

        if (x < N) {
            const int o = y * N + x;
            if (x == 0 || y == 0 || x == N - 1 || y == N - 1) {
                output[o] = centre;
            } else {
                output[o] = 0.5f * centre
                        + 0.125f * (row[sx - 1] + row[sx + 1] + above + below);
            }
        }
        __syncthreads();
    }