GPU Computing · Kernel Programming

Thread Hierarchy and Global Index Arithmetic in a Kernel Launch

How a grid of blocks and threads is flattened into buffer indices, where the arithmetic overflows, and what the launch dimensions deliberately do not guarantee.

A kernel launch is described by three numbers: how many blocks form the grid, how many threads form each block, and how many bytes of dynamic shared memory each block receives. Every index a kernel computes derives from those numbers plus a handful of built-in variables. Almost all of the practical difficulty of kernel programming starts here, because the mapping from a two- or three-dimensional launch to a one-dimensional buffer is done by hand, and hand-written index arithmetic fails in two characteristic ways: it addresses outside the buffer, or it wraps around silently and addresses the wrong place inside it.

Three coordinates and the units they name

A launch is a grid of blocks; a block is a group of threads. The built-ins gridDim, blockDim, blockIdx and threadIdx are each a three-component vector of unsigned integers, and a launch that leaves a component unset gets 1 for it. The components are independent, so a grid of 64×64 blocks of 16×16 threads is legal and is not the same shape as a grid of 4096 blocks of 256 threads, even though both cover 65536 threads in total.

Three limits shape the arithmetic. A block may contain at most 1024 threads. The x component of the grid may reach $2^{31}-1$; the y and z components are limited to 65535 each. And threads are executed in groups of 32 — a warp — so a block of 100 threads is issued as four warps, the last of which has 4 active lanes and 28 that are masked off. Nothing in the language requires a block size that is a multiple of 32, but a block size that is not wastes issue slots on every instruction that the partial warp executes.

The block, not the thread, is the unit of scheduling. A block is placed on one streaming multiprocessor and runs there until it finishes; it is never migrated or preempted mid-execution in the ordinary model. That is what makes __syncthreads() meaningful: every thread that can reach the barrier is on the same processor, so the arrival count is a hardware counter in that processor’s scheduler rather than a protocol over memory, and the memory ordering the barrier implies is scoped to that processor’s own caches. It also means threads in different blocks have no shared storage and no ordering relationship, which is the constraint that all cross-block algorithms have to work around.

Flattening a multi-dimensional problem

Buffers are linear. A two-dimensional image is a contiguous run of rows, and the element at column $x$ of row $y$ sits at offset $y \cdot W + x$ where $W$ is the row stride. The launch coordinates therefore have to be folded into that offset. The convention that costs nothing is to let threadIdx.x be the fastest-varying coordinate, because consecutive lanes then touch consecutive addresses, and the memory system services a warp’s accesses as a small number of wide transactions instead of many scattered ones.

__global__ void scale2d(float* out, const float* in,
                        int width, int height, float k) {
    // Cast before multiplying: blockIdx.x * blockDim.x overflows 32 bits
    // long before a large grid is reached.
    const size_t x = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
    const size_t y = (size_t)blockIdx.y * blockDim.y + threadIdx.y;
    if (x >= (size_t)width || y >= (size_t)height) return;

    const size_t i = y * (size_t)width + x;
    out[i] = in[i] * k;
}

The multiplication order in that first line is not stylistic. blockIdx.x is a 32-bit unsigned value, and blockDim.x is another; their product is computed in 32-bit arithmetic before the assignment to size_t would have widened anything. Casting the first operand forces the whole expression into 64-bit arithmetic. A grid of 4 million blocks of 256 threads covers about $2^{30}$ threads, and the last block’s base index is just under $2^{30}$, which fits comfortably; a grid of 17 million blocks of the same size covers more than $2^{32}$ threads, and the indices of the final blocks wrap. The failure is not a crash. The index wraps to a small value and the kernel writes to the beginning of the buffer, which is why out-of-range indexing is usually found by a memory-checking tool rather than by observing incorrect output.

The row stride is not always the width

When rows are padded — for alignment, or because a sub-rectangle of a larger image is being processed — the distance between the start of one row and the next is the pitch, and the pitch is greater than or equal to the width. The offset formula becomes $y \cdot \text{pitch} + x$, and any code that uses the logical width where the pitch belongs will shear the image by one element per row. CUDA exposes allocation and copy routines that carry a pitch alongside the pointer for exactly this case, and a kernel that receives such a pointer must receive the pitch as a separate parameter; there is no way to recover it from the pointer alone.

The same distinction appears in the other direction when a kernel processes a sub-tile: the tile’s own dimensions are small, but the stride used to step from one tile row to the next is the full image width. Confusing the two produces the classic half-image output, where the upper-left corner is correct and everything after the first row is displaced.

Where the arithmetic breaks

Three failure points account for most index defects.

The first is width. A 32-bit signed index covers about 2.1 billion elements; a buffer of 3 billion single-precision floats is 12 GB, which is an ordinary size on a device with 40 or 80 GB of memory. Any index declared int in such a kernel is wrong regardless of how carefully the launch is configured. Declaring indices size_t or unsigned long long fixes the storage, but not the intermediate products, which must also be widened.

The second is rounding up. Computing the grid size for $n$ elements with $t$ threads per block requires $\lceil n/t \rceil$, usually written (n + t - 1) / t. When $n$ is within $t-1$ of the maximum value of its type, the addition overflows and the result is a grid far too small, so the tail of the problem is never processed. The overflow-free form adds the deficit instead: n / t + (n % t != 0). This is the same rounding that governs the bounds checks discussed in tail handling.

The third is the product of two dimensions. A 4096×4096 image has 16.7 million elements, comfortably inside 32 bits; a 65536×65536 image has 4.3 billion, which is not. The per-axis coordinates fit in 32 bits while the linear offset does not, so the widening has to happen at the point of combination, not at the point of declaration.

What the launch numbers do not promise

The dimensions specify how many threads exist, not when or where they run. There is no ordering guarantee between blocks: block 7 may begin before block 3, may run concurrently with it, and may be placed on any multiprocessor. A kernel that has block $k$ write a value and block $k+1$ read it has no correctness argument at all in the ordinary model, because nothing forces the write to happen first. Algorithms that need cross-block ordering either launch separate kernels, which imposes ordering at the launch boundary, or use the explicit grid-wide mechanisms covered in grid-wide barriers and persistent work queues.

The block size is also not a purely semantic choice. Changing it changes how many warps each block occupies, how much shared memory is claimed per block, and therefore how many blocks fit on a multiprocessor. A kernel that is correct at 256 threads and correct at 1024 threads can differ substantially in achieved throughput, because occupancy and the amount of work available to hide latency both change. The index arithmetic, however, must be invariant to the block size — if a kernel only produces correct results for one particular launch configuration, the arithmetic is relying on an accident.

The relationship between the two index forms

Kernels are usually written against a global linear index rather than against the launch coordinates directly. The transformation is one line:

const size_t tid = (size_t)blockIdx.x * blockDim.x + threadIdx.x;

This is the 1-D case of the flattening above, and it is the form that makes a kernel independent of the launch shape: any grid whose total thread count covers the problem will do, provided every thread computes the same linear index from its own coordinates. That property is what allows the launch size to be tuned for the machine rather than for the data, and it is the basis of the loop structure described in grid-stride loops. The cost of the abstraction is that the kernel no longer knows its position within a row, which matters when the algorithm needs neighbours in two dimensions; those kernels keep the coordinate form and pay for the extra arithmetic.

Other programming models express the same hierarchy with different names. OpenCL calls the grid the NDRange, a block a work-group, and a thread a work-item, with get_global_id, get_local_id and get_group_id returning the corresponding coordinates; the work-group is likewise the scope in which local memory is shared and barriers are legal. WebGPU compute shaders declare the workgroup size in the shader with an attribute and read global_invocation_id and local_invocation_id built-ins. The numbers differ — WebGPU caps a workgroup at 256 invocations in its baseline limits, well under the 1024 threads a CUDA block may hold — but the arithmetic and its hazards are identical, because the underlying machine is the same shape.

References

  1. CUDA C++ Programming Guide — NVIDIA
  2. OpenCL C Specification, Version 1.2 — Work-item Functions — Khronos Group
  3. PTX ISA: Parallel Thread Execution ISA — NVIDIA