Robust Grid Stride Loops

Welcome back! So far, we’ve used a "one-thread-per-element" approach, which assumes our grid is always large enough to cover the data. But what happens if the dataset is massive or the GPU hardware limits our grid size?

To build professional, "battle-ready" CUDA code, we use the Grid Stride Loop. This pattern allows threads to process multiple elements by "stepping" through the array. By the end of this lesson, your kernels will be hardware-agnostic—capable of handling any data size correctly, regardless of how many threads you launch. Let’s upgrade our indexing for maximum robustness!

So far, we have aimed for “one thread per element.” That works well, but it assumes we launch enough threads to cover all $N$ items. However, relying on a fixed 1:1 mapping between threads and data points is not always feasible or optimal in production. What if your dataset size exceeds the hardware limits of the GPU, or you want to process millions of elements using only a modest number of threads for better cache reuse?

A grid stride loop addresses this by changing the mindset: instead of a thread handling exactly one element, each thread handles element i, then jumps ahead by a fixed stride and handles i + stride, i + 2*stride, and so on, until we run out of data.

This pattern is conceptually similar to a modulo operation, where each thread $T$ takes responsibility for all data indices $i$ such that $i \pmod{\text{grid\_stride}} = T$.

It helps in two common cases:

• When $N$ is much larger than the total threads we launch.
• When we intentionally limit grid size for better performance tuning later.

Now we add the key idea of this lesson. We compute the same starting index, but then we loop forward by a stride of gridDim.x * blockDim.x.

__global__ void vectorAdd(const float* a, const float* b, float* c, int n) {
    int i = blockDim.x * blockIdx.x + threadIdx.x;
    int stride = blockDim.x * gridDim.x;

    for (; i < n; i += stride) {
        c[i] = a[i] + b[i];
    }
}

Why this is more robust:

• Even if the grid has fewer threads than n, threads “come back around” to handle later elements.
• The loop condition i < n becomes our safety check, so every write is still in bounds.
• We can tune blocksPerGrid for performance later without breaking correctness.

Choosing the right pattern depends on your goals:

While the simple approach is great for learning, the grid stride loop is the industry standard because it ensures your kernel won't crash when ported to a different GPU with different hardware limits.

Next, the program creates a large problem size and prepares host and device memory. Notice how we compute bytes once, since CUDA allocation and copies use byte counts.

int main() {
    const int N = 50000;
    size_t bytes = N * sizeof(float);

    std::vector<float> h_a(N, 1.0f), h_b(N, 2.0f), h_c(N, 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));

Here we allocate:

• h_a, h_b as inputs, filled with constants.
• h_c as output, initialized to zero.
• d_a, d_b, d_c as device arrays with matching sizes.

We copy inputs to the GPU, then compute launch dimensions using a ceiling division. This ensures full coverage when using the baseline one-index-per-thread style.

    CUDA_CHECK(cudaMemcpy(d_a, h_a.data(), bytes, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), bytes, cudaMemcpyHostToDevice));

    int threadsPerBlock = 256;
    int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;

    std::cout << "Launching " << N << " operations using "
              << blocksPerGrid << " blocks." << std::endl;

The key formula is:

With grid stride loops, this calculation is still fine, but it is no longer the only way to stay correct, because the loop can cover the “missing” work.

Now we launch the kernel, check for launch errors, and synchronize. The same launch works whether we use the baseline kernel or the grid stride version.

    vectorAdd<<<blocksPerGrid, threadsPerBlock>>>(d_a, d_b, d_c, N);
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

Two important checks happen here:

• cudaGetLastError() catches configuration and launch issues.
• cudaDeviceSynchronize() ensures the kernel is finished before we copy results back.

Finally, we copy d_c back, verify every element is 3.0f, print success or failure, and free device memory.

    CUDA_CHECK(cudaMemcpy(h_c.data(), d_c, bytes, cudaMemcpyDeviceToHost));

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

    std::cout << "Final Result: " << (success ? "SUCCESS" : "FAILURE") << std::endl;

    CUDA_CHECK(cudaFree(d_a));
    CUDA_CHECK(cudaFree(d_b));
    CUDA_CHECK(cudaFree(d_c));
    return success ? 0 : 1;
}

This verification step performs real work for us: it confirms that our indexing logic and our launch setup actually produced correct results for all 50,000 elements.

When everything is configured correctly, the program prints the launch size and the final verification result.

Launching 50000 operations using 196 blocks.
Final Result: SUCCESS

This tells us two things: we computed a grid large enough to cover N, and the CPU-side check confirmed that every output element matched the expected sum.

We now have a stronger mental model for 1D CUDA work: global index gives each thread a starting point, and a grid stride loop lets that thread safely handle more elements by stepping forward with a fixed stride. This pattern makes kernels correct across many sizes, and it keeps working even when we change the grid size for tuning.

Next, we will move into practice problems where we apply grid stride loops and launch math ourselves, and you will see how quickly this becomes a dependable habit.