GPU CUDA Programming — Kernels, Threads, Blocks, and Grids

What the model actually asks of you

The CUDA programming model is small. You write a function marked __global__, you call it with an extra pair of triple angle brackets saying how many copies to run, and you write the body as if it handles exactly one element of the problem. Which physical core runs which copy, in what order, and how many run at once are the runtime's business - and deliberately unspecified.

That last part is the contract, and it is what makes the model scale. When you launch a grid of blocks you are promising that the blocks are mutually independent: any block may run before, after, or simultaneously with any other, on any SM, and a correct kernel produces the same answer under every schedule. In exchange, one binary saturates a laptop GPU with twenty SMs and a datacenter part with a hundred and thirty-two without a recompile. Break the promise - have block 7 spin waiting on something block 9 writes - and you get a kernel that passes on one device and deadlocks on another, because nothing guarantees block 9 is ever resident at the same time.

The hardware that the contract is written against - streaming multiprocessors, warps, the memory hierarchy, occupancy - is the GPU architecture overview. This article is the other side of the boundary: what you actually type, and why the syntax has the shape it does.

Host, device, and what a launch really is

A .cu file is compiled twice. nvcc splits it, hands the host half to the system compiler (MSVC, gcc, clang) and compiles the device half itself. Four qualifiers decide which half a function lands in:

• __global__ - a kernel. Called from the host, runs on the device, must return void. The only thing launchable with <<< >>>.
• __device__ - callable only from device code, and aggressively inlined into whatever kernel calls it.
• __host__ - an ordinary CPU function. This is the default, so you rarely write it.
• __host__ __device__ - compiled into both binaries from one body, which is how utility math gets shared between a reference implementation and a kernel.

The launch configuration has four slots, not two: kernel<<<gridDim, blockDim, dynamicSharedBytes, stream>>>(args...). The first two are dim3, a three-component struct that converts implicitly from an integer, which is why <<<n, 256>>> is legal and means dim3(n,1,1) and dim3(256,1,1). The third reserves dynamic shared memory; the fourth picks a stream. Both default to zero, and zero for the stream means the default stream, which has semantics worth knowing about before you rely on them.

Kernel arguments are copied by value into a small dedicated parameter space - a few kilobytes on current hardware - at launch time. Pass device pointers, not structures containing arrays. Handing a kernel a host pointer compiles without a murmur and faults on the first dereference; it is the most common first-day CUDA bug, and the error will not be reported at the launch line.

The launch is asynchronous. kernel<<<...>>>() appends a command to a queue and returns, typically in a few microseconds, and the CPU carries on. The mandatory void return is not an oversight - there is nothing to return to, because the caller left. Everything a kernel produces it writes into memory you allocated in advance.

Indexing - mapping threads onto data

Inside a kernel, four built-ins tell a thread who it is. threadIdx is its position within its block, blockIdx the block's position within the grid, blockDim the block's shape, gridDim the grid's shape. All four are dim3, all four are read-only, and every thread executes the identical code with different values in them. That is the whole trick: one program, several million distinct index values.

__global__ void saxpy(int n, float a, const float* x, float* y) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) y[i] = a * x[i] + y[i];
}

// host side
int block = 256;
int grid  = (n + block - 1) / block;     // ceiling division
saxpy<<<grid, block>>>(n, a, d_x, d_y);

Two details carry the weight. The ceiling division means the final block is partly idle whenever n is not a multiple of the block size, and the if (i < n) guard is the only thing stopping those surplus threads from writing past the end of the array. Omitting the guard rarely crashes - it quietly corrupts whichever allocation happens to sit next in device memory, and the symptom appears in an unrelated tensor.

The choice of x is not cosmetic. Threads are bundled into warps in linear order with the x component varying fastest, so adjacency in threadIdx.x is adjacency in warp lane, which is what turns x[i] into a single wide memory transaction. For two-dimensional data, put the contiguous axis on x:

int col = blockIdx.x * blockDim.x + threadIdx.x;   // fastest-varying axis
int row = blockIdx.y * blockDim.y + threadIdx.y;
if (row < h && col < w) {
    int idx = row * pitch + col;
    out[idx] = f(in[idx]);
}

Swap x and y here and the kernel stays correct and gets several times slower, because each warp now walks a column and touches a different cache line per lane. That mechanism, and the transaction arithmetic behind it, is coalesced memory access. The shape limits are worth memorising: a block holds at most 1024 threads with the z extent capped at 64, while a grid's x extent runs to about two billion blocks and its y and z extents to 65535 each.

The grid-stride loop

The saxpy above hard-wires one thread to one element, which means the grid size is a function of the input size. The idiom that library kernels use instead decouples them:

__global__ void saxpy_gs(int n, float a, const float* x, float* y) {
    int stride = blockDim.x * gridDim.x;          // total threads in the grid
    for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += stride)
        y[i] = a * x[i] + y[i];
}

Each thread walks the array in steps equal to the grid's total thread count. The pattern buys several things at once. The launch geometry becomes a tuning parameter chosen for the device - blocks per SM times the SM count - rather than dictated by the data, so a 300-element input and a 300-million-element input use the same well-sized grid. The kernel is correct for any n, including an n smaller than one block, with no special case. Per-thread setup - address computation, loading a constant, initialising an accumulator - is amortised over several elements instead of paid once per element. And because the stride is the whole grid width, lanes within a warp still touch consecutive addresses on every iteration, so coalescing survives; a naive alternative in which each thread takes a contiguous run of elements would destroy it.

It also makes debugging tractable. Launch the same kernel as <<<1, 1>>> and it degenerates to a correct serial loop over the whole array, which is a fast way to separate an indexing bug from a race.

The costs are a loop counter, a comparison, and a register or two. For a one-shot kernel over a fixed-size buffer, the direct form is fine. For anything reusable, write the grid-stride version.

Choosing a block size

Block size is the first number people guess at and the one with the most constrained answer. Four rules cover most cases.

Make it a multiple of 32. The hardware allocates threads to warps 32 at a time. A block of 100 threads occupies four warp slots and leaves four lanes of the last warp permanently masked off - you pay for 128 threads' worth of scheduling resources and use 100.

Start at 128 or 256. These sit in the sweet spot for most kernels: large enough that per-block overhead and any shared-memory tile amortise, small enough that the block's register and shared-memory footprint does not lock out co-resident blocks.

1024 is a ceiling, not a goal. The whole block's registers must fit in the SM's register file simultaneously, so a register-hungry kernel launched with 1024 threads per block can fail outright with cudaErrorLaunchOutOfResources - a launch-time error that returns nothing unless you check for it. Bigger blocks are also not more parallelism: total threads is grid times block, and moving threads from the grid dimension into the block dimension only makes the scheduler's job harder.

Let the runtime pick when you do not care. cudaOccupancyMaxPotentialBlockSize takes a kernel and returns the block size and minimum grid size that maximise theoretical occupancy for the device it runs on, which is a better default than a hard-coded 256 in portable code. In the other direction, __launch_bounds__(256) on the kernel promises the compiler you will never exceed that block size, letting it cap register allocation accordingly.

The grid size has its own trap. If a kernel launches 140 blocks on a 132-SM device and each SM hosts one at a time, the first wave uses every SM and the second wave uses eight - the kernel takes twice as long as the work implies while 124 SMs idle. This wave quantisation is another argument for grid-stride kernels, where you can round the grid to a multiple of the resident block capacity. How registers, shared memory and slot limits combine into that capacity is the occupancy article; how register demand gets away from you is GPU register pressure.

__syncthreads() and the uniformity rule

__syncthreads() is a barrier across one thread block. It gives two guarantees: no thread proceeds past it until every thread in the block has arrived, and every shared- and global-memory write issued before it is visible to every thread in the block after it. Those two together are what make a load-into-shared-memory, then-read-shared-memory tiling pattern correct.

The rule that breaks programs is that all threads in the block must reach the same barrier. A __syncthreads() inside divergent control flow is undefined behaviour, and the practical outcome ranges from correct-by-luck to a hang.

// WRONG - threads with tid >= n never reach the barrier
if (tid < n) {
    tile[tid] = in[tid];
    __syncthreads();
    out[tid] = tile[tid] + tile[tid ^ 1];
}

// RIGHT - the barrier is unconditional, the work is guarded
if (tid < n) tile[tid] = in[tid];
__syncthreads();
if (tid < n) out[tid] = tile[tid] + tile[tid ^ 1];

The same failure hides inside loops with data-dependent trip counts: if the barrier sits in a loop body and different threads run different numbers of iterations, the barriers no longer line up. Hoist the loop bound to something uniform across the block. compute-sanitizer --tool synccheck finds these mechanically and is worth running on any kernel that uses shared memory.

Below the block: warps are no longer implicitly in lockstep

Older CUDA code omitted barriers once the active thread count dropped to 32, on the reasoning that a warp executes in lockstep anyway. Since the Volta generation, independent thread scheduling gives every thread its own program counter, and that assumption is no longer safe. The replacements are explicit: __syncwarp(mask) for a warp-level barrier, and the _sync family of intrinsics - __shfl_down_sync, __ballot_sync, __any_sync - which take an explicit mask of participating lanes. Warp-synchronous code without those masks is not merely fragile, it is wrong on current hardware.

Above the block: there is no grid barrier

Nothing in the ordinary model synchronises blocks with each other. The end of a kernel is the grid-wide barrier, and the standard way to express a two-phase algorithm is two kernel launches. Cooperative groups add an explicit grid.sync(), but only for kernels started with a cooperative launch, which caps the grid at the number of blocks that fit resident on the device simultaneously - a real constraint, and the reason most production code still uses two launches.

Five address spaces, from the programmer side

CUDA does not hide the memory hierarchy behind a cache the way a CPU does. Where a variable lives is something you declare, and the declaration determines its scope, its lifetime, and its speed.

Local memory is the misleading one. Despite the name it is not on-chip; it is a per-thread slice of device DRAM, cached like anything else, with global-memory latency. Two things put variables there: register spills when the kernel's demand exceeds its allocation, and any local array indexed with a value the compiler cannot resolve at compile time. float acc[8] touched only by a fully-unrolled loop stays in registers; the identical array indexed by a runtime variable silently moves to local memory and the kernel gets several times slower with no warning in the source. Compile with --ptxas-options=-v and it will tell you the spill bytes.

Shared memory is the on-chip scratchpad you manage yourself. Static __shared__ arrays size at compile time; a single extern __shared__ char buf[]; declaration takes its size from the third launch parameter, which is how one kernel serves several tile sizes. Blocks get roughly 48 KB by default, and going beyond that requires opting in explicitly with cudaFuncSetAttribute(kernel, cudaFuncAttributeMaxDynamicSharedMemorySize, bytes) before the launch. Because the array is physically banked, the index expression you write determines whether a warp's access completes in one cycle or serialises - that whole subject, including padding and swizzling, is the shared memory article.

Constant memory is a 64 KB read-only region with a broadcast-optimised path: when every lane of a warp reads the same address it costs about one register read, and when lanes read different addresses it serialises. That makes it right for coefficients, small lookup tables and configuration structs, and wrong for anything indexed by thread. For large read-only arrays the equivalent hint is marking the pointer const T* __restrict__, which lets the compiler route loads through the read-only path. Relative latency and bandwidth for each of these levels is the memory hierarchy article.

Unified memory versus explicit copies

The explicit path is cudaMalloc plus cudaMemcpy: two allocations, two pointers, and you decide when bytes cross the link. The managed path is cudaMallocManaged, which returns a single pointer valid on both sides; the driver migrates pages on demand when the other side touches them.

Unified memory is the right tool in three situations. When you are bringing up a port and want the thing to run before it runs fast. When the data structure is pointer-based - a tree, a graph, anything where an explicit copy means a deep copy with pointer fixups. And when the working set genuinely exceeds device memory and you would rather the program degrade than fail, since oversubscription is handled by eviction rather than an allocation error.

It is the wrong default in a steady-state loop. The first device touch of a managed page is a fault serviced across the interconnect, and a loop that alternates host and device access to the same buffer migrates it every iteration - correct, and dramatically slower than one explicit copy plus a kernel. The middle ground is managed memory plus hints: cudaMemPrefetchAsync(ptr, bytes, device, stream) moves pages ahead of the kernel that needs them, converting a storm of faults into one bulk transfer, and cudaMemAdvise with cudaMemAdviseSetReadMostly or cudaMemAdviseSetPreferredLocation stops the driver bouncing a page that both sides read.

One detail belongs to the explicit path regardless. Pageable host memory cannot be DMA'd, so cudaMemcpyAsync from an ordinary malloc buffer stages through an internal pinned buffer and is not truly asynchronous - the call blocks for part of the transfer and no overlap happens. cudaMallocHost or cudaHostAlloc gives page-locked memory the copy engine can read directly, and that is the prerequisite for overlapping transfer with compute. Pinned memory is a scarce system-wide resource, though: pin a large fraction of host RAM and you degrade the whole machine, including processes that have nothing to do with the GPU.

Streams, events, and honest timing

Every operation you enqueue - kernel, copy, memset, event - goes into a stream, an ordered queue. Operations in the same stream execute in issue order. Operations in different streams have no ordering relative to each other and may overlap, which is the entire basis of copy/compute overlap.

Omit the stream argument and the work goes to the default stream. Under legacy semantics the default stream implicitly synchronises with every other blocking stream in the process, which is why a program that has been carefully converted to use streams but still leaves half its calls on the default stream shows no overlap at all in a profile. Compiling with --default-stream per-thread changes this, giving each host thread its own non-synchronising default stream, and it is usually what you want in multi-threaded code.

Timing is the other place the asynchrony bites. A CPU-side timer wrapped around a launch measures the enqueue - a few microseconds, unrelated to how long the kernel takes - and adding a synchronise to fix that changes what you are measuring by draining the pipeline. The correct instrument is a pair of events:

cudaEvent_t beg, end;
cudaEventCreate(&beg); cudaEventCreate(&end);

cudaEventRecord(beg, stream);
my_kernel<<<grid, block, 0, stream>>>(args);
cudaEventRecord(end, stream);

cudaEventSynchronize(end);                 // wait for the marker, not the CPU
float ms; cudaEventElapsedTime(&ms, beg, end);

Events are markers placed in the stream, so the interval they measure is device time between two points in the same queue - which is what you wanted. Record both events in the same stream, and remember that the first launch of any kernel includes module load and possibly JIT compilation, so discard a warm-up iteration before believing a number.

The patterns built on top of this - multi-stream double buffering, expressing cross-stream dependencies with events, and capturing a repeated launch sequence into a graph so the CPU stops paying per-kernel launch cost - are CUDA streams and graphs. Stream priorities, MPS and how the hardware arbitrates between concurrent work are GPU scheduling.

Error handling: why cudaGetLastError after a launch is not optional

A kernel launch has no return value to check. Errors arrive by two different routes, and conflating them is why CUDA bugs get reported at the wrong line.

Synchronous launch errors are detected while the command is being enqueued: an invalid block dimension, a request for more dynamic shared memory than the device allows, a kernel with no compiled image for this device, insufficient registers for the requested block size. These are recorded immediately and returned by the next CUDA call you make - or by cudaGetLastError() if you ask right away.

Asynchronous execution errors are raised while the kernel is running: an illegal address, a misaligned access, a failed assert. The CPU is long past the launch line by then, so these surface at the next synchronising call, which may be several kernels and a hundred lines later.

#define CUDA_CHECK(x) do {                                              \
    cudaError_t err__ = (x);                                            \
    if (err__ != cudaSuccess) {                                         \
        fprintf(stderr, "%s:%d %s\n", __FILE__, __LINE__,               \
                cudaGetErrorString(err__));                             \
        abort();                                                        \
    } } while (0)

my_kernel<<<grid, block>>>(args);
CUDA_CHECK(cudaGetLastError());        // launch configuration errors, here
CUDA_CHECK(cudaDeviceSynchronize());   // execution errors (debug builds)

cudaGetLastError() returns the stored error and clears it; cudaPeekAtLastError() returns it without clearing. That clearing behaviour is why the check has to go after every launch rather than once at the end of a phase: an uncleared error from an earlier call is what the next check will report, and you will spend an afternoon on an innocent kernel.

Some errors are sticky. After cudaErrorIllegalAddress the CUDA context is destroyed; every subsequent API call in the process returns the same error and no recovery is possible short of exiting. So a stack trace pointing at some harmless cudaMemcpy means the real fault happened in an earlier kernel. The tool for pinning it down is the environment variable CUDA_LAUNCH_BLOCKING=1, which makes every launch synchronous so errors are reported at the launch that caused them. It is a debugging setting only - it serialises everything and destroys performance.

From source to SASS - nvcc, PTX, and fat binaries

The device half of your source goes through two compilers, and knowing which is which explains most deployment failures.

The front end emits PTX, a virtual instruction set targeted at a virtual architecture such as compute_80. PTX is stable and forward compatible. Then ptxas compiles PTX into SASS, the real machine encoding for a real architecture such as sm_80. SASS is where register allocation and instruction scheduling actually happen, and it is not portable across GPU generations.

-gencode arch=compute_80,code=sm_80 embeds SASS for one architecture. Repeat the flag for several architectures and you get a fat binary carrying multiple SASS images; add code=compute_80 and it carries the PTX as well. At load time the driver looks for SASS matching the device. If it finds none but there is embedded PTX, it JIT-compiles that PTX to SASS on the spot - the forward-compatibility path that lets a binary built before a GPU shipped still run on it, at the cost of a compile on first load. The result is cached on disk, controlled by CUDA_CACHE_PATH and disabled with CUDA_CACHE_DISABLE=1. If the binary has neither matching SASS nor PTX you get no kernel image is available for execution on the device, which is nearly always a missing -gencode for the deployment GPU rather than anything wrong with the code.

Compute capability is not just a version stamp. Language features are gated on it: __syncwarp and the _sync intrinsics need 7.0 or later, and the asynchronous-copy and newer tensor-core paths need later generations still. The -arch you pass decides which intrinsics compile at all, not merely how fast the result is.

Two flags earn their place in any build you intend to tune. --ptxas-options=-v prints registers, shared memory and spill bytes per kernel - the numbers every occupancy discussion starts from. And -rdc=true enables relocatable device code so __device__ functions can be called across translation units; it costs inlining opportunities and usually a few registers, so leave it off unless you need device-side linking.

Debugging entry points before you reach a profiler

When a kernel produces wrong numbers, correctness tools come before performance tools. compute-sanitizer runs the real binary under instrumentation and has four modes worth knowing:

• the default memcheck catches out-of-bounds and misaligned device accesses, with a source line if you compiled with -lineinfo;
• racecheck catches shared-memory races - almost always a missing or misplaced __syncthreads();
• initcheck catches reads of device global memory that was never written;
• synccheck catches illegal barrier and warp-intrinsic usage, including the non-uniform __syncthreads() above.

It is slow - an order of magnitude is normal - but it converts "wrong sometimes" into a file and a line, which no amount of staring at index arithmetic will.

printf works inside kernels, with caveats. Output goes to a fixed-size circular buffer flushed at the next synchronisation, so lines can be dropped silently, and interleaving across thousands of threads makes ordering meaningless. Guard it - if (blockIdx.x == 0 && threadIdx.x == 0) - or you will produce a hundred thousand useless lines. Device-side assert also works and raises a sticky error, so a failed assertion takes the context down with it.

cuda-gdb gives real breakpoints inside kernels with the ability to focus on a specific thread or warp. It needs -G, which turns off device optimisation and changes register allocation entirely - so a -G build tells you nothing about performance. For anything you intend to measure, build with -lineinfo instead: optimisations stay on and Nsight Compute can still attribute stalls back to source lines. Where to go from there - timeline capture with Nsight Systems, per-kernel counters with Nsight Compute, and how to read them without jumping to conclusions - is GPU profiling.

When to write a kernel at all

Most people should not. A matrix multiply from cuBLAS, a convolution from cuDNN, or a CUTLASS template instantiation carries years of tuning and sits close to hardware peak; a hand-written version will lose, and lose again on the next architecture. Compilers cover much of the rest: torch.compile, XLA and TensorRT fuse elementwise chains and matmul epilogues automatically, which removes the most common reason to write a kernel by hand - see kernel fusion.

The cases that still justify hand-written kernels are specific. The operation is memory-bound and the compiler will not fuse it the way the algorithm needs. The data layout is unusual - a custom sparsity pattern, a packed quantisation format - so no library kernel matches. You need an epilogue nobody exposes. Or the operations are tiny and numerous, and launch overhead dominates everything else.

Triton is the useful middle rung. You write a Python function over blocks rather than threads, and the compiler takes care of index arithmetic, coalescing and pipeline scheduling. For elementwise, reduction and attention-shaped kernels it reaches a large fraction of hand-tuned performance at a fraction of the surface area, and it is where most new custom kernels in the ML stack now get written. Dropping to CUDA C++ or CUTLASS is what you do when you need explicit control over tensor-core fragments, asynchronous copy pipelines, or warp specialisation.

The order that wastes the least time: profile first and confirm the kernel is actually the bottleneck; try a library; try a compiler; write it in Triton; and only then reach for CUDA C++.