The "hello world" of GPU programming is a kernel that adds two vectors. It is three lines long, and nearly every line hides an idea that decides whether real kernels are fast or slow. Which element does each thread work on? What happens when the vector's length is not a multiple of the number of threads? Why does it matter that neighboring threads read neighboring addresses? Answer those three questions properly and you understand most of what CUDA asks of you.

This section teaches CUDA's execution model, the thread, the block, and the grid, by writing kernels in CUDA C++ and, because CUDA cannot run in this page, by simulating each one in NumPy with the same indices. Then it covers the two memory rules that matter most, coalescing for global memory and bank conflicts for shared memory, counted exactly with small simulators. It ends with synchronization and with calling your own kernel from PyTorch.

Threads, blocks, and grids

A CUDA kernel is a function that runs once per thread. You launch it on a grid of blocks, each made of the same number of threads. Inside the kernel, every thread can read four built-in variables that tell it who it is.

• threadIdx, its position within its block, from 0 to blockDim - 1.
• blockIdx, its block's position within the grid, from 0 to gridDim - 1.
• blockDim and gridDim, the sizes of the block and the grid.

Each is a three-component vector (.x, .y, .z), so grids and blocks can be one-, two-, or three-dimensional. That is purely for convenience, a way of naming threads that matches the shape of the data. The hardware sees a block as a flat list of threads, numbered with x fastest, then y, then z, and cut into warps of 32 (§11.1).

The hierarchy is not arbitrary. A block is the unit that is placed on one SM, so its threads can cooperate, sharing the SM's fast shared memory and waiting for each other at barriers. Different blocks may run on different SMs, in any order, possibly not at the same time, so blocks cannot wait for each other within a kernel. (Two later additions relax this in controlled ways. Hopper's thread block clusters guarantee that up to 8 blocks, or 16 with an opt-in, run at the same time on neighboring SMs, where they can synchronize and read each other's shared memory. A cooperative launch of a grid small enough to be resident all at once allows a grid-wide barrier. Ordinary kernels use neither.) A block holds at most 1,024 threads. A grid can hold billions of blocks, and the hardware streams them onto the SMs as earlier blocks finish.

In plain words, a grid is a big bag of independent blocks, and each block is a small team of threads that can talk to each other.

A vector add

Here is the whole program, kernel and host code. The kernel computes each thread's global index from its block and thread indices, checks that the index is inside the vector, and does one addition.

#include <cstdio>
#include <cstdlib>
#include <cuda_runtime.h>

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

__global__ void vector_add(const float* a, const float* b, float* out, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;   // this thread's global index
    if (i < n) {                                     // the last block can run past the end
        out[i] = a[i] + b[i];
    }
}

int main() {
    const int n = 1000;
    const size_t bytes = n * sizeof(float);
    float* h_a = (float*)malloc(bytes);
    float* h_b = (float*)malloc(bytes);
    float* h_out = (float*)malloc(bytes);
    for (int i = 0; i < n; ++i) { h_a[i] = i; h_b[i] = 2.0f * i; }

    float *d_a, *d_b, *d_out;                        // device (GPU) pointers
    CUDA_CHECK(cudaMalloc(&d_a, bytes));
    CUDA_CHECK(cudaMalloc(&d_b, bytes));
    CUDA_CHECK(cudaMalloc(&d_out, bytes));
    CUDA_CHECK(cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_b, h_b, bytes, cudaMemcpyHostToDevice));

    const int threads = 256;
    const int blocks = (n + threads - 1) / threads;  // ceil(n / threads) = 4
    vector_add<<<blocks, threads>>>(d_a, d_b, d_out, n);
    CUDA_CHECK(cudaGetLastError());                  // launch errors
    CUDA_CHECK(cudaMemcpy(h_out, d_out, bytes, cudaMemcpyDeviceToHost));  // waits for the kernel

    printf("out[999] = %f\n", h_out[999]);           // 2997.0
    cudaFree(d_a); cudaFree(d_b); cudaFree(d_out);
    free(h_a); free(h_b); free(h_out);
    return 0;
}

Read it line by line. __global__ marks a function that runs on the device and is launched from the host, while __device__ marks a helper that runs on the device and is called from kernels (the warp reduction of §11.4 is one), and __host__ __device__ compiles a function for both. The triple angle brackets <<<blocks, threads>>> set the grid and block sizes. The launch returns immediately, and the cudaMemcpy back to the host waits for the kernel to finish, because both are queued on the same stream. Errors in the launch itself (an invalid configuration, say) appear only if you ask with cudaGetLastError, and errors inside the kernel appear at the next synchronizing call. Production code checks every call.

The index formula blockIdx.x * blockDim.x + threadIdx.x numbers all threads in the grid consecutively. With 1,000 elements and blocks of 256 we need $\lceil 1000/256\rceil = 4$ blocks, 1,024 threads, and the last 24 have nothing to do. Without the if (i < n), those 24 threads would read and write past the end of the arrays. In CUDA that does not raise an error. It silently corrupts whatever lives there, or crashes some later, unrelated kernel.

The same kernel, simulated

We can run exactly the same logic in NumPy. The simulator below loops over blocks, as the hardware's block scheduler does, and handles all threads of a block at once with a vector of threadIdx values, as the SIMT hardware does.

Grid-stride loops

The vector add launches one thread per element. An alternative launches a fixed grid, sized to fill the GPU (a few blocks per SM), and lets each thread loop over the data with a stride equal to the total number of threads.

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

This grid-stride loop works for any $n$ with any launch configuration, reuses each thread's setup (handy when a thread first loads constants), and keeps every iteration coalesced, because in each pass of the loop, consecutive threads still touch consecutive elements. Notice the int64_t. A tensor with more than $2^{31}$ elements, eight GB of float32, overflows a 32-bit index, a bug that appears only on the largest inputs.

import numpy as np

def grid_stride_indices(n, grid_dim, block_dim):
    S = grid_dim * block_dim
    return [np.array([g]) if g < n else np.array([], int) for g in range(S)]   # placeholder: one pass only

Two-dimensional grids

For a matrix, a 2D grid of 2D blocks is the natural naming scheme. The universal convention is that x indexes columns and y indexes rows.

__global__ void scale_matrix(float* A, float alpha, int M, int N) {   // A is M x N, row-major
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    if (row < M && col < N) {
        A[row * N + col] *= alpha;
    }
}
// launch: dim3 block(32, 8); dim3 grid((N + 31) / 32, (M + 7) / 8);
//         scale_matrix<<<grid, block>>>(A, 2.0f, M, N);

Why columns on x? Because x is the fastest-varying index within a block, so the 32 threads of a warp have consecutive threadIdx.x and the same threadIdx.y. With x mapped to columns of a row-major matrix, a warp reads 32 consecutive elements of one row, which is the access pattern memory likes best, as we are about to count. Swap the roles and each warp reads down a column, touching 32 different rows.

The bounds check now has two parts, and both are needed whenever $M$ or $N$ is not a multiple of the block size. Getting this mapping right, with the right ceilings, is a classic interview question.

import numpy as np

def matrix_launch(M, N, bx, by):
    gx, gy = -(-M // bx), -(-N // by)             # placeholder: x and y swapped
    BY, BX, TY, TX = np.meshgrid(np.arange(gy), np.arange(gx), np.arange(by), np.arange(bx), indexing="ij")
    rows = BX * bx + TX
    cols = BY * by + TY
    return (gx, gy), rows, cols, (rows < M) & (cols < N)

Memory coalescing

When a warp executes a load, the hardware collects the 32 addresses and serves them with as few memory transactions as possible. Global memory is accessed in sectors of 32 bytes (grouped into 128-byte cache lines). What a warp pays for is the number of distinct sectors its 32 addresses fall into, not the number of bytes it asked for.

• If 32 threads read 32 consecutive float32 values starting at an aligned address, the 128 bytes fall in exactly 4 sectors. Every byte moved is used. This is a coalesced access.
• If the same 32 values start one element later, they straddle 5 sectors, and 20 percent of the traffic is wasted.
• If each thread reads every second float, the warp touches 256 bytes, 8 sectors, and half the traffic is wasted.
• If the stride is 8 floats (32 bytes) or more, every thread hits its own sector, 32 sectors for 128 useful bytes, an efficiency of 12.5 percent.

In plain words, memory is delivered in 32-byte pieces, and a warp should ask for whole pieces. Neighboring threads should read neighboring addresses.

The cell counts sectors for each pattern. Reading a column of a row-major matrix is the stride-$N$ case, and it is exactly what a naive transpose or a badly mapped 2D kernel does.

Every kernel's effective bandwidth is the peak times this efficiency. A kernel with the stride-8 pattern can never exceed one eighth of the HBM bandwidth, however well it is otherwise written. Move the sliders to see the lanes and the sectors they touch.

Two remarks complete the picture. The L1 and L2 caches can rescue a strided pattern if other warps soon read the rest of each sector, which is why measured efficiency is sometimes better than this count. And vectorized loads, where each thread reads 8 or 16 consecutive bytes with one instruction (float2, float4), keep the access coalesced while issuing a quarter as many instructions, a standard trick in memory-bound kernels.

import numpy as np

def coalescing(byte_addrs, elem_bytes, sector=32):
    n = len(byte_addrs) * elem_bytes // sector        # placeholder: assumes perfect coalescing
    return n, 1.0

Shared memory and bank conflicts

Shared memory is the SM's programmer-managed scratchpad, declared with __shared__ and visible to every thread of a block. It is how threads of a block exchange data and how kernels reuse data without going back to HBM. Physically it is split into 32 banks, each 4 bytes wide, with consecutive 4-byte words in consecutive banks. Word $w$ lives in bank $w \bmod 32$. Each bank can deliver one word per cycle, so a warp's request is served in one cycle when its 32 lanes hit 32 different banks.

If several lanes ask for different words in the same bank, those requests are serialized, and a request in which some bank receives $k$ distinct words takes $k$ cycles, a $k$-way bank conflict. Lanes that ask for the same word do not conflict, because the word is broadcast to all of them.

The pattern that bites in practice is the column of a tile. Suppose a block stores a $32\times32$ tile of floats in shared memory, float tile[32][32], and a warp reads one column, tile[lane][c]. Lane $j$ reads word $32j + c$, which is in bank $(32j + c) \bmod 32 = c$ for every lane, a 32-way conflict. The classic fix is to pad each row by one word, float tile[32][33]. Now lane $j$ reads word $33j + c$, in bank $(33j + c) \bmod 32 = (j + c) \bmod 32$, a different bank for every lane.

In plain words, a stride that is a multiple of 32 words sends every lane to the same bank, and padding the row length to 33 spreads them back out. More generally, any odd stride is conflict-free.

The rules for 8- and 16-byte accesses (a float2 or float4 per lane) are more involved, since the hardware then splits a warp's request into phases of 16 or 8 lanes, but the principle is the same. Swizzled layouts, which permute columns with an XOR of the row index instead of padding, avoid conflicts without wasting space, and are what libraries like CUTLASS and Triton use for tensor-core tiles.

import numpy as np

def bank_conflicts(word_addrs, banks=32):
    w = np.asarray(word_addrs) % banks
    return int(np.bincount(w, minlength=banks).max())    # placeholder: treats broadcasts as conflicts

Synchronization

Threads of a block cooperate through shared memory, and cooperation needs ordering. __syncthreads() is a barrier. No thread of the block passes it until every thread has reached it, and every shared-memory write before it is visible to every thread after it. Consider a kernel that reverses each block's slice of an array through shared memory.

__global__ void reverse_blocks(const float* in, float* out, int n) {   // assumes n % blockDim.x == 0
    extern __shared__ float s[];                     // size set at launch (third <<< >>> argument)
    int t = threadIdx.x;
    int i = blockIdx.x * blockDim.x + t;
    s[t] = in[i];                                    // phase 1: every thread writes one element
    __syncthreads();                                 // wait until the whole block has written
    out[i] = s[blockDim.x - 1 - t];                  // phase 2: read an element another thread wrote
}
// reverse_blocks<<<n / 256, 256, 256 * sizeof(float)>>>(in, out, n);

Each thread reads an element written by a different thread, usually in a different warp. Without the barrier, a warp that reaches phase 2 early reads slots that other warps have not filled yet. The result depends on scheduling, may be right on one run and wrong on the next, and is called a race condition. The simulation runs warps one after another, which is one legal schedule, with and without the barrier.

Without the barrier, the first half of the block reads slots no one has written yet. Three rules keep barriers safe. Every thread of the block must reach the same __syncthreads(), so never put one inside a branch that only some threads take, which can deadlock. A barrier synchronizes one block only, and in an ordinary launch there is no barrier across blocks inside a kernel (the kernel boundary is the global barrier, apart from the clusters and cooperative launches mentioned earlier). And when threads of different blocks must combine results, they use atomic operations such as atomicAdd, which make a read-modify-write indivisible. Atomics on floats make the order of additions, and so the rounding, vary from run to run, which is one source of nondeterminism in PyTorch (§0.9, and the reductions of §11.4).

Putting it together, a transpose

A matrix transpose shows all three ideas at once. Reading rows and writing columns directly makes one side uncoalesced. The fix is to stage a tile through shared memory. Each block reads a $32\times32$ tile with coalesced row reads, waits at a barrier, and writes the transposed tile with coalesced row writes, reading the tile column-wise from shared memory, which is where the padding matters.

#define TILE 32
// out = in^T, where in is M x N (row-major) and out is N x M. Launch with block(32, 8),
// grid((N + 31) / 32, (M + 31) / 32); each thread copies TILE / 8 = 4 elements.
__global__ void transpose(const float* in, float* out, int M, int N) {
    __shared__ float tile[TILE][TILE + 1];           // +1 column: no bank conflicts on column reads
    int x = blockIdx.x * TILE + threadIdx.x;         // column in `in`
    int y = blockIdx.y * TILE + threadIdx.y;         // row in `in`
    for (int j = 0; j < TILE; j += blockDim.y)
        if (x < N && y + j < M) tile[threadIdx.y + j][threadIdx.x] = in[(y + j) * N + x];
    __syncthreads();
    x = blockIdx.y * TILE + threadIdx.x;             // column in `out`
    y = blockIdx.x * TILE + threadIdx.y;             // row in `out`
    for (int j = 0; j < TILE; j += blockDim.y)
        if (x < M && y + j < N) out[(y + j) * M + x] = tile[threadIdx.x][threadIdx.y + j];
}

Both global accesses have consecutive threadIdx.x touching consecutive addresses, and the column read from tile is conflict-free thanks to the extra column. A well-written transpose runs close to the copy bandwidth of the GPU, while the naive version can be several times slower.

Calling your kernel from PyTorch

You rarely write a whole CUDA program. You write a kernel and call it on PyTorch tensors. The quickest route is torch.utils.cpp_extension.load_inline, which compiles a C++ and CUDA source string at runtime, binds the functions you name, and caches the build. This needs a CUDA GPU and the CUDA toolkit, so it is shown, not run.

A few details in that listing are easy to get wrong. The kernel is launched on PyTorch's current stream, so it is ordered correctly with the surrounding PyTorch operations. The wrapper checks the device, the dtype, and the shapes, and calls .contiguous(), because the kernel indexes memory as a dense array and would silently read the wrong elements of a transposed view (strides, §0.1). And the result is checked against PyTorch's own a + b. For anything beyond experiments, build the same code ahead of time as a CUDAExtension in setup.py, register it as a custom operator so torch.compile can see it, and give it a backward pass with torch.autograd.Function (§0.9). To time it, synchronize first, as §0.9 warned, or use CUDA events.

In PyTorch

The CPU has its own version of coalescing, the cache line. The block reads the same number of float32 values from a large matrix twice, once as a contiguous block of columns and once with a stride of 16 elements (64 bytes, one cache line per element), and times both.

Read more

• NVIDIA. CUDA C++ Programming Guide. The reference for the programming model, built-in variables, launch limits, shared memory, and synchronization.
• NVIDIA. CUDA C++ Best Practices Guide. Coalescing, shared memory, and bank conflicts, with the official rules and diagrams.
• Kirk, Hwu & El Hajj (2022). Programming Massively Parallel Processors, 4th edition. Chapters 2, 3, and 6 develop vector add, multidimensional grids, and memory coalescing carefully.
• Harris (2013). An Efficient Matrix Transpose in CUDA C/C++. The tiled, padded transpose and its measured bandwidth, step by step.
• PyTorch. Custom C++ and CUDA Extensions. Building and loading extensions, including the ahead-of-time setup.py route.
• NVIDIA. Compute Sanitizer. The tool that finds out-of-bounds accesses, races, and missing synchronization in your kernels.

Watch next

▶

▶