Thread Indexing and Grid-Stride Loops

Inside a kernel you are one thread among many, and the only thing distinguishing you from your neighbours is a set of built-in variables describing where you sit. Turning those into a single index into your data is the first line of almost every kernel ever written.

int i = blockIdx.x * blockDim.x + threadIdx.x;
            

All four have .x, .y and .z members. If you launched a one-dimensional grid, only .x is meaningful and the others are 1.

Grids are built out of whole blocks, so unless your data size happens to be an exact multiple of the block size, the last block contains threads with no work to do. Their computed index runs off the end of your array.

int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < n) {
    c[i] = a[i] + b[i];
}
            

The block count itself uses a round-up so that no element is left uncovered:

int numBlocks = (n + threadsPerBlock - 1) / threadsPerBlock;
            

Block size is one of the few tuning knobs available at launch, and there are only a few sensible answers.

Make it a multiple of the execution group size, which is 32 on Nvidia and 64 on AMD CDNA. 128, 256 and 512 are divisible by both and are the usual defaults; 256 is a reasonable first guess for almost anything. A block of 32 is a full warp on Nvidia but only half a wavefront on CDNA, which is the classic portability trap.

Beyond that, larger blocks mean fewer blocks and more threads sharing one compute unit's registers and scratchpad, which can reduce how many blocks stay resident. There is no formula that beats measuring, and the answer changes with the kernel. Start at 256, try 128 and 512, keep whichever wins.

For naturally two-dimensional data, launching a 2D grid keeps the index arithmetic readable:

__global__ void transpose(const float *in, float *out, int width, int height)
{
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if (x < width && y < height) {
        out[x * height + y] = in[y * width + x];
    }
}

// launched with, say:
dim3 block(16, 16);
dim3 grid((width + 15) / 16, (height + 15) / 16);
            

One thread per element seems natural and has a weakness: the grid size is dictated by the data size. Launch a billion-element problem and you launch a billion threads, most of which do one addition and retire, paying scheduling cost for almost no work.

The alternative is to launch a grid sized to the machine rather than the problem, and have each thread loop:

__global__ void vecAdd(const float *a, const float *b, float *c, int n)
{
    int stride = blockDim.x * gridDim.x;
    for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += stride) {
        c[i] = a[i] + b[i];
    }
}
            

It works for any n. The same kernel handles one element or a billion, with no assumption that the grid covers the data.

You can size the grid deliberately. Launching roughly enough blocks to fill the device, rather than however many the data implies, gives more predictable occupancy and lets you reuse a tuned launch configuration across problem sizes.

It is debuggable. Launch it with a single block of one thread and the loop degenerates to a sequential pass over the whole array, which is an easy way to check correctness independently of the parallel decomposition.

The bounds check is still there, folded into the loop condition. It never goes away; it just stops needing a separate line.