Kernel

A kernel is the unit of CUDA code that programmers typically write and compose, roughly the GPU counterpart of a procedure or function in CPU languages. The difference is that a kernel is called ("launched") once and returns once, but is executed many times: once by each of many threads.

How a kernel executes#

The executions of a kernel are generally concurrent (their order is non-deterministic) and parallel (they happen at the same time on different execution units).

All the threads running a kernel form a kernel grid, also called a thread block grid. This is the top level of the thread hierarchy. A grid runs across multiple Streaming Multiprocessors (SMs), so it operates at the scale of the whole GPU. The matching level of the memory hierarchy is global memory.

In CUDA C++, the host invokes a kernel with pointers to global memory on the device. The kernel returns nothing; it just mutates memory.

Example: matrix multiplication#

The "hello world" of CUDA kernels is multiplying two square matrices, A and B. Below are two implementations that differ in how they map the textbook algorithm onto the thread hierarchy and memory hierarchy.

Naive version: one thread per output element#

This kernel is inspired by the first matmul kernel in Programming Massively Parallel Processors (4th edition, Figure 3.11). Each thread computes one element of the output matrix on its own: it loads each element of one row of A and one column of B into registers, multiplies the pairs, sums the results and writes the sum back to global memory.

__global__ void mm(float* A, float* B, float* C, int N) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < N && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < N; k++) {
            sum += A[row * N + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}

Each thread does one floating point operation (FLOP) per read from global memory: a multiply and an add for a load from A and a load from B. That will never use the whole GPU, because the arithmetic bandwidth of the CUDA Cores in FLOPs/s is far higher than the memory bandwidth between GPU RAM and the SMs.

Tiled version: using shared memory#

We can raise the ratio of FLOPs to memory operations by mapping the work more carefully onto the hierarchies. This "tiled" kernel, inspired by Figure 5.9 of the same book, loads submatrices of A and B into shared memory and assigns the computation of submatrices of C to thread blocks.

#define TILE_WIDTH 16

__global__ void mm(float* A, float* B, float* C, int N) {

    // declare variables in shared memory ("smem")
    __shared__ float As[TILE_WIDTH][TILE_WIDTH];
    __shared__ float Bs[TILE_WIDTH][TILE_WIDTH];

    int row = blockIdx.y * TILE_WIDTH + threadIdx.y;
    int col = blockIdx.x * TILE_WIDTH + threadIdx.x;

    float c_output = 0;
    for (int m = 0; m < N/TILE_WIDTH; ++m) {

        // each thread loads one element of A and one of B from global memory into smem
        As[threadIdx.y][threadIdx.x] = A[row * N + (m * TILE_WIDTH + threadIdx.x)];
        Bs[threadIdx.y][threadIdx.x] = B[(m * TILE_WIDTH + threadIdx.y) * N + col];

        // we wait until all threads in the 16x16 block are done loading into smem
        // so that it contains two 16x16 tiles
        __syncthreads();

        // then we loop over the inner dimension,
        // performing 16 multiplies and 16 adds per pair of loads from global memory
        for (int k = 0; k < TILE_WIDTH; ++k) {
            c_output += As[threadIdx.y][k] * Bs[k][threadIdx.x];
        }
        // wait for all threads to finish computing
        // before any start loading the next tile into smem
        __syncthreads();
    }
    C[row * N + col] = c_output;
}

Each outer-loop iteration loads two elements, then the thread runs 16 inner-loop iterations of a multiply and an add: 16 FLOPs per global memory read.

Key idea Kernel performance depends heavily on how work maps onto the thread and memory hierarchies; reusing data from shared memory raises FLOPs per global memory read.

Going further#

This is still far from a fully optimized matmul. This worklog by Si Boehm of Anthropic walks through optimizations that push the FLOP-to-memory-read ratio higher and map the algorithm more tightly onto the hardware. The two kernels here resemble his Kernel 1 and Kernel 3; the worklog covers ten.

Both the worklog and this article only target the CUDA Cores. The fastest matrix multiplication kernels instead run on Tensor Cores, which have much higher arithmetic bandwidth.

Related#

• Thread Block Grid: the set of threads a kernel launch creates
• Shared Memory: the key to the tiled kernel
• Arithmetic Intensity: the FLOPs-per-byte ratio the tiling improves
• CUDA C++: the language kernels are usually written in