2D Grid Stride Loops

Welcome back to 2D Grids and Matrix Math. We are now at Lesson 3, and that means we already know two key ideas from this course: how to map 2D thread coordinates to matrix positions, and how to use that mapping inside a real matrix kernel. In this lesson, we take the next step and make that pattern scalable.

As you may recall from previous lessons, we often launch enough threads so that each matrix element gets one thread. That is simple and clear, but it is not always enough for very large data. Here, we will build a 2D grid-stride loop so a smaller grid can still cover a much larger matrix correctly.

A standard 1-to-1 mapping fails if the grid is smaller than the matrix. While we usually aim for one thread per element, launching a smaller grid is often necessary for:

• Hardware Limits: Matrices can exceed the maximum grid dimensions allowed by the GPU.
• Portability: It ensures the kernel works regardless of the specific GPU hardware or the input size.
• Efficiency: Reusing threads can improve cache hits and reduce the overhead of managing massive numbers of blocks.

A 2D grid-stride loop fixes this by letting each thread "step" through the matrix at regular gaps:

Instead of stopping after one element, each thread moves by strideX across columns and strideY across rows until the work is done. This makes the kernel robust across different matrix sizes and GPU architectures.

Let's start by defining the kernel itself. The first part should look familiar: each thread computes its starting column and row from block and thread indices.

__global__ void matrixAdd2DStride(const float* A, const float* B, float* C,
                                  int width, int height) {
    int startCol = blockIdx.x * blockDim.x + threadIdx.x;
    int startRow = blockIdx.y * blockDim.y + threadIdx.y;

There are two important details here. First, A and B are marked const because the kernel only reads from them. Second, width and height describe the matrix shape, while startCol and startRow tell us where this thread begins. So far, this is the same coordinate mapping pattern we already know. The new idea will come from what the thread does after computing that starting point.

Here is the heart of the lesson. Instead of touching one element and stopping, each thread moves through the matrix in two dimensions. This is what allows a small launch to process a much larger matrix.

    int strideX = blockDim.x * gridDim.x;
    int strideY = blockDim.y * gridDim.y;

    for (int r = startRow; r < height; r += strideY) {
        for (int c = startCol; c < width; c += strideX) {
            int index = r * width + c;
            C[index] = A[index] + B[index];
        }
    }
}

A few ideas are packed into this small block:

• strideX is the total number of threads covering the x direction, and strideY is the same for y.
• The loops stop naturally at the matrix boundaries, so no separate bounds check is needed here.
• index = r * width + c uses the same row-major mapping from earlier lessons.
• Each visited element is just matrix addition: C[index] = A[index] + B[index].

With the kernel ready, we can set up a large input in main(). The program uses a 2048 x 2048 matrix, which is big enough to show why this technique matters.

int main() {
    const int W = 2048;
    const int H = 2048;
    const size_t bytes = W * H * sizeof(float);

    std::vector<float> h_A(W * H, 1.0f);
    std::vector<float> h_B(W * H, 2.0f);
    std::vector<float> h_C(W * H, 0.0f);

    float *d_A, *d_B, *d_C;
    CUDA_CHECK(cudaMalloc(&d_A, bytes));
    CUDA_CHECK(cudaMalloc(&d_B, bytes));
    CUDA_CHECK(cudaMalloc(&d_C, bytes));

    CUDA_CHECK(cudaMemcpy(d_A, h_A.data(), bytes, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_B, h_B.data(), bytes, cudaMemcpyHostToDevice));

This test is designed to make verification easy. Every value in h_A is 1.0f and every value in h_B is 2.0f, so every output value should become 3.0f. We also allocate three device buffers and copy only the inputs to the GPU. The output buffer d_C does not need an initial host value because the kernel will overwrite every valid position it reaches.

Now comes the interesting part: we intentionally launch a grid that is much smaller than the matrix. This would fail in a one-pass kernel, but it works here because of the 2D grid-stride loops.

    dim3 threadsPerBlock(16, 16);
    dim3 numBlocks(16, 16);

    std::cout << "Matrix Size: " << W << "x" << H
              << " (" << W * H << " elements)" << std::endl;
    std::cout << "Grid Size: " << numBlocks.x * threadsPerBlock.x << "x"
              << numBlocks.y * threadsPerBlock.y << " threads" << std::endl;
    std::cout << "Launching 2D grid-stride kernel..." << std::endl;

    matrixAdd2DStride<<<numBlocks, threadsPerBlock>>>(d_A, d_B, d_C, W, H);

    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

The printed lines help us see the mismatch clearly: the matrix has over four million elements, but the launch covers only 256 x 256 thread positions at a time. That is exactly the scenario where the stride logic earns its place. The error checks after launch are also important; they catch launch failures and then wait until the GPU has finished before we continue.

Matrix Size: 2048x2048 (4194304 elements)
Grid Size: 256x256 threads
Launching 2D grid-stride kernel...

After the kernel finishes, we copy the result back, check every element, print the status, and free device memory. This completes the full workflow from host setup to GPU execution to verification.

    CUDA_CHECK(cudaMemcpy(h_C.data(), d_C, bytes, cudaMemcpyDeviceToHost));

    bool success = true;
    for (int i = 0; i < W * H; ++i) {
        if (std::fabs(h_C[i] - 3.0f) > 1e-5f) {
            success = false;
            break;
        }
    }

    std::cout << "2D Grid-Stride verification: "
              << (success ? "SUCCESS" : "FAILURE") << std::endl;

    CUDA_CHECK(cudaFree(d_A));
    CUDA_CHECK(cudaFree(d_B));
    CUDA_CHECK(cudaFree(d_C));

    return success ? 0 : 1;
}

The verification loop checks all W * H values, not just a sample. Using std::fabs with 1e-5f is a safe floating-point habit, even though this example is very simple. If the kernel handled the whole matrix correctly, the program prints the success message below, which means every output element matched the expected value 3.0f.

2D Grid-Stride verification: SUCCESS

In this lesson, we extended the 2D indexing pattern from earlier lessons into something much more practical: a 2D grid-stride loop. We saw how a thread computes a starting row and column, how strideX and strideY let it revisit more work, and how this makes a matrix kernel scale beyond the size of the launched grid. We also walked through the host code that builds a large test, launches a deliberately small grid, verifies every result, and frees GPU memory safely.

This pattern is simple, but it is also powerful. It lets us write kernels that are less tied to one exact launch shape and more ready for real workloads. In the practice section ahead, we will put this idea to work so we can write scalable matrix kernels with real confidence.