Thread hierarchy: index arithmetic for grid / block / thread

CUDA's parallel model is a two-level hierarchy: a grid is made of blocks, a block is made of threads. The split is not cosmetic. It maps onto hardware: a block is assigned as a whole to one SM. Threads in a block can share shared memory and wait for each other with __syncthreads(). Threads in different blocks share nothing.

                        grid
  ┌───────────────┬───────────────┬───────────────┐
  │   block 0     │   block 1     │   block 2     │
  │ t0 t1 t2 t3   │ t0 t1 t2 t3   │ t0 t1 t2 t3   │
  └───────────────┴───────────────┴───────────────┘
    0  1  2  3      4  5  6  7      8  9 10 11     ← global index

  global = blockIdx.x * blockDim.x + threadIdx.x
             ↑             ↑              ↑
        which block    threads/block   thread in block

1D: vector add

1__global__ void vecAdd(const float* a, const float* b, float* c, int n) {2    int i = blockIdx.x * blockDim.x + threadIdx.x;3    if (i < n) {              // bounds check: do not skip this4        c[i] = a[i] + b[i];5    }6}

Why must you write if (i < n)? Block thread counts are usually a power of two such as 128 or 256, while n is an arbitrary value. To cover every element, grid size must round up, so some threads in the last block have an index past the array. Skip that check and you get an out-of-bounds device-memory write: wrong results in the mild case, a smashed neighboring data structure in the bad case.

1int n = 1'000'000;2int threads = 256;3int blocks = (n + threads - 1) / threads;   // = ceil(n / threads) = 39074 5vecAdd<<<blocks, threads>>>(d_a, d_b, d_c, n);6// 3907 * 256 = 1'000'192 threads. The last 192 are stopped by the if

grid-stride loop: a more reliable pattern

One thread per element is a poor fit when the data is very large or very small. A grid-stride loop decouples grid size from data size: no matter how many threads you launch, each thread walks the whole array with a fixed stride. The same kernel can then run any size of input with one grid size tuned to the hardware.

1__global__ void vecAddStride(const float* a, const float* b, float* c, int n) {2    int stride = gridDim.x * blockDim.x;         // total threads in the grid3    for (int i = blockIdx.x * blockDim.x + threadIdx.x;4         i < n;5         i += stride) {6        c[i] = a[i] + b[i];7    }8}

• ▸Size independent: the same grid configuration finishes correctly whether n is 1000 or a billion.
• ▸Accesses stay coalesced: neighboring threads still touch neighboring addresses in each iteration, so you do not break memory coalescing (covered in the next stage).
• ▸Easier to debug: set the grid to <<<1, 1>>> and it collapses to serial execution, so you can check results against a known reference.
• ▸Reuse registers: compute constants once outside the loop and spread that cost across many elements.

2D: images and matrices

A 2D grid is the natural fit for matrices. One convention matters: threadIdx.x should map to the contiguous direction inside a row (the column index). Neighboring threads in a warp are numbered .x first, so coalescing only happens when .x walks contiguous memory. Reverse that mapping and performance can drop by 5x or more.

1__global__ void scaleMatrix(float* m, int rows, int cols, float k) {2    // x → column (contiguous in memory), y → row3    int col = blockIdx.x * blockDim.x + threadIdx.x;4    int row = blockIdx.y * blockDim.y + threadIdx.y;5 6    if (row < rows && col < cols) {7        m[row * cols + col] *= k;     // row-major: one warp covers contiguous cols8    }9}10 11// host side12dim3 block(32, 8);                                    // 256 threads, x dimension is exactly one warp13dim3 grid((cols + block.x - 1) / block.x,14          (rows + block.y - 1) / block.y);15scaleMatrix<<<grid, block>>>(d_m, rows, cols, 2.0f);

How to pick a block size