GPU Kernels, Explained: CUDA, Triton, cuBLAS, torch.compile and KernelBench

When you run a neural network in PyTorch, Python is only the coordinator. Behind every tensor operation — from a simple activation to a multi-head attention block — sits one or more GPU kernels: specialized functions executed by tens of thousands of hardware threads in parallel. Heavy matrix multiplications call into NVIDIA’s hand-tuned cuBLAS library, convolutions route to cuDNN, and custom elementwise ops run either as kernels built into PyTorch or as JIT-compiled code generated by torch.compile. In practice, whether your training loop takes one second per step or three seconds is decided almost entirely inside these kernels.

In this deep dive, we look at how GPU kernels actually work under the hood. We will build the same ReLU kernel in raw CUDA and in Triton, use the classic roofline model to see why most deep-learning operations are choked by memory bandwidth rather than arithmetic, and walk through four core optimizations: operator fusion, streaming online softmax, tiled matrix multiplies with Tensor Cores, and GEMM epilogues. Finally, we compare the four ways kernels reach hardware in modern PyTorch pipelines and examine KernelBench — the benchmark testing whether LLMs can write production kernels that beat PyTorch.

Note on numbers: Latency figures in this post are theoretical lower bounds derived from NVIDIA H100 SXM datasheet specifications to illustrate hardware limits, rather than live empirical benchmarks.

Start here

What a GPU kernel is and how it maps work onto threads; the same kernel in CUDA and Triton; the roofline model and why fusion, tiling and Tensor Cores matter; what cuBLAS, cuDNN, eager PyTorch and torch.compile each do; how KernelBench scores a kernel and how to compare its numbers.

• Kernel: A function that runs on the GPU. The CPU launches it once; the GPU runs many copies of it in parallel, each on a different part of the data.
• Thread, block, grid: A thread runs one copy of the kernel. Threads are grouped into blocks that share fast on-chip memory; all blocks of one launch form the grid.
• SM: Streaming Multiprocessor: one of the GPU's processor cores. An H100 SXM has 132 of them; each runs several blocks at a time.
• HBM: High-bandwidth memory: the GPU's main memory (80 GB on an H100). Large, but slow compared with registers and on-chip shared memory.
• Memory-bound / compute-bound: A kernel is memory-bound when moving its data takes longer than its arithmetic, and compute-bound in the opposite case.
• Fusion: Doing several operations in one kernel so intermediate results stay on chip instead of making a round trip through HBM.
• Tensor Core: A matrix-multiply unit inside each SM, many times faster than ordinary arithmetic units, at reduced input precision (TF32, BF16, FP8).

1. What a kernel is

A common mental model is that the GPU somehow runs your Python program. In reality, Python stays entirely on the CPU. The CPU acts as an orchestrator: PyTorch figures out which math operations are needed, and the CUDA driver launches a compiled C++ or PTX function onto the GPU, specifying a grid configuration that dictates how many parallel copies to spawn.

That GPU function is the kernel. Each running instance is a thread. Threads are grouped into thread blocks (up to 1,024 threads per block), and the hardware executes them in lockstep groups of 32 known as warps. Every thread in a warp executes the exact same instruction at the same time across different data elements. (If you want a broader overview of how CPUs and accelerators divide work, our guide to CPUs, GPUs and clusters introduces the hardware layout).

To understand why kernel performance varies so wildly, you have to look at the physical silicon. An NVIDIA H100 SXM chip contains 132 Streaming Multiprocessors (SMs). Each SM has its own fast registers and configurable on-chip shared memory (SRAM), while all SMs share a 50 MB L2 cache and 80 GB of High-Bandwidth Memory (HBM3). As you move down this memory hierarchy, capacity grows, but latency and energy costs multiply:

An H100 SXM in the numbers that matter for kernels. The CPU only launches the kernel; blocks of threads are spread over the 132 SMs, and every value they use travels up a memory hierarchy that gets smaller and faster closer to the arithmetic.

The central tension in GPU programming is the imbalance between raw compute power and memory bandwidth. An H100 SXM can churn through 67 TFLOPS of standard FP32 math on its CUDA cores and nearly 1,000 TFLOPS of BF16 matrix arithmetic on its Tensor Cores. Yet its memory bus caps out at 3.35 TB/s.

Do the math: to keep standard CUDA cores 100% saturated, a kernel must perform roughly 20 arithmetic operations for every single byte loaded from HBM. For Tensor Cores, that requirement shoots up to nearly 300 operations per byte. Yet everyday deep-learning operations — activations, additions, LayerNorms — perform only 1 or 2 operations per byte. The GPU’s massive compute engines end up starved, waiting on memory traffic. Once you see this gap, modern GPU kernel design suddenly makes sense: almost every trick is about moving less data across HBM.

2. A CUDA kernel

CUDA C++ is NVIDIA’s native dialect for programming accelerators. In CUDA, you write code from the viewpoint of a single thread. Built-in runtime variables — blockIdx, blockDim, and threadIdx — allow each thread to compute its unique position in the global array.

Here is how you write an elementwise ReLU ( yi​=max(xi​,0) ) for an array of n floating-point numbers:

__global__ void relu_kernel(const float* x, float* y, int n) { int i = blockIdx.x * blockDim.x + threadIdx.x; // this thread's element if (i < n) y[i] = fmaxf(x[i], 0.0f); // the last block runs past n } // Host side: 256 threads per block, enough blocks to cover n elements. relu_kernel<<<(n + 255) / 256, 256>>>(x, y, n); 

On the host side, the launch syntax <<<...>>> configures how many 256-thread blocks to schedule so that all n elements are covered. Because array sizes rarely align perfectly to multiples of 256, the defensive guard if (i < n) stops out-of-bounds memory writes on the final block.

Notice how thread i reads element x[i]. Because consecutive threads in a warp access adjacent memory addresses in parallel, the GPU memory controller can merge their 32 individual 4-byte reads into a single 128-byte transaction. This is called memory coalescing. If your access pattern is strided or scattered, memory throughput can drop by 5× to 10× without changing the mathematical output.

ReLU is straightforward because each data point is completely independent. Real-world machine learning kernels, however, require threads to collaborate. In a matrix multiplication, threads must stage input tiles into on-chip shared memory (a shared scratchpad inside the SM), coordinate through __syncthreads() barriers to prevent data races, issue specialized Tensor Core instructions, and carefully budget registers per thread to maintain high SM occupancy. CUDA gives you bare-metal control over all of these knobs — but it also means every misplaced barrier or register spill can silently tank your performance.

3. The same kernel in Triton

Writing raw CUDA can feel like writing C with inline assembly: you are constantly juggling thread IDs, shared memory pointers, and race conditions. This complexity is why Triton (Tillet, Kung and Cox, MAPL 2019; developed by OpenAI and the open-source community) has become so popular in the AI systems world.

In Triton, you write Python. More importantly, you stop thinking about individual scalar threads. Instead, you write from the perspective of an entire program block processing small array tiles (typically 64 to 1,024 elements) using vectorized NumPy-like operations. The Triton compiler takes care of mapping blocks to warps, vectorizing loads, ensuring coalesced access, and orchestrating shared memory:

import torch import triton import triton.language as tl @triton.jit def relu_kernel(x_ptr, y_ptr, n, BLOCK: tl.constexpr): pid = tl.program_id(0) # this program's block  offs = pid * BLOCK + tl.arange(0, BLOCK) # BLOCK positions at once  mask = offs < n # the last block runs past n  x = tl.load(x_ptr + offs, mask=mask) tl.store(y_ptr + offs, tl.maximum(x, 0.0), mask=mask) def relu(x: torch.Tensor) -> torch.Tensor: y = torch.empty_like(x) n = x.numel() relu_kernel[(triton.cdiv(n, 1024),)](x, y, n, BLOCK=1024) return y 

The same 20-element ReLU. CUDA assigns one value per thread and leaves cooperation to the programmer; Triton assigns one block per program and leaves the mapping onto threads to the compiler.

What you decide and what the toolchain decides.

What do you give up with Triton? You give up fine-grained control. Individual program instances cannot communicate with each other except through slow global HBM memory, and low-level hardware tricks (like hand-crafted warp shuffles or irregular memory layouts) remain out of reach. But for the fused, memory-bound operations that dominate modern Transformers, that trade-off is an easy win: you get 80–95% of hand-tuned CUDA speed in a few dozen lines of readable Python. That is why PyTorch chose Triton as the primary code generation target for torch.compile.

4. Memory-bound or compute-bound: the roofline

Before spending days attempting to optimize a kernel, you need to answer a fundamental question: what is actually holding it back? Is the GPU waiting on arithmetic logic units to finish crunching numbers, or are the execution units idling while waiting for tensors to travel across the memory bus?

The classic roofline model (Williams, Waterman and Patterson, 2009) formalizes this intuition into a simple formula:

Here, $I$ represents the arithmetic intensity — the ratio of floating-point math operations performed per byte of data moved across HBM. Every GPU architecture has a distinct transition point called the ridge point ( Ppeak​/BHBM​ ):

• If your kernel sits to the left of the ridge (low intensity), it is memory-bound. Faster arithmetic units will not speed it up by a single nanosecond. The only way to win is to move fewer bytes across the bus or reuse data in cache.
• If your kernel sits to the right of the ridge (high intensity), it is compute-bound. Memory bandwidth is plenty fast; your runtime is determined by how well you saturate ALUs or Tensor Cores.

On an NVIDIA H100 SXM, this ridge line sits at 20 FLOPs/byte for FP32 on CUDA cores, roughly 148 FLOPs/byte for TF32 on Tensor Cores, and 295 FLOPs/byte for BF16 Tensor Cores (using dense datasheet peaks; NVIDIA’s marketing often cites 2:4 structured sparsity numbers, which are double these).

Arithmetic intensity of common operations in FP32, counting the minimum traffic (each input read once, each output written once).

Notice where standard deep-learning building blocks fall: activations, normalizations, and softmax sit two full orders of magnitude to the left of the ridge. For these layers, optimizing math formulas is irrelevant; all that matters is memory traffic. Conversely, large matrix multiplications sit comfortably to the right: they offer abundant theoretical data reuse, so the challenge becomes organizing memory access so Tensor Cores never stall. Use the interactive explorer below to inspect these bounds across different tensor shapes:

The site has an interactive version with three tabs (fusion, roofline, matmul tiling) and sliders for element count and dtype — try it on the original post.

5. Four optimizations, from simple to hard

5.1 Fusion: stop writing intermediates to memory

Consider a familiar one-liner: y = relu(x * s + b). In standard eager PyTorch, this statement executes as three distinct kernel launches:

1. Kernel 1 (scale): reads x from HBM, computes x * s, and writes a temporary tensor tmp1 back to HBM.
2. Kernel 2 (bias add): reads tmp1 and b, computes tmp1 + b, and writes temporary tensor tmp2 back to HBM.
3. Kernel 3 (activation): reads tmp2, applies relu, and writes final tensor y back to HBM.

Notice the waste: tmp1 and tmp2 exist solely to pass numbers to the next line of code, yet they force millions of bytes across the memory bus. A fused kernel collapses all three steps into a single pass: it reads x once, keeps intermediate values in fast on-chip registers, and writes y once:

@triton.jit def scale_bias_relu(x_ptr, b_ptr, y_ptr, n, cols, s, BLOCK: tl.constexpr): offs = tl.program_id(0) * BLOCK + tl.arange(0, BLOCK) mask = offs < n x = tl.load(x_ptr + offs, mask=mask) b = tl.load(b_ptr + offs % cols, mask=mask) # one bias per column, stays in cache  tl.store(y_ptr + offs, tl.maximum(x * s + b, 0.0), mask=mask) 

Three eager kernels move six tensors through HBM; one fused kernel moves two. For 16.8M FP32 values the bandwidth bound falls from 120 µs to 40 µs, and two kernel launches disappear.

Fusion does not eliminate a single FLOP of arithmetic. What it eliminates is memory traffic and driver launch overhead. That is why operator fusion is always the first optimization compilers apply, and why cleaning up chains of small elementwise layers yields such dramatic speedups.

5.2 Softmax: fuse a reduction, then make it streaming

Softmax is a notorious bottleneck in attention layers. The standard formula — yi​=exi​−m/∑j​exj​−m with the row maximum $m$ subtracted for numerical stability — is built from two reductions (max and sum) and several elementwise operations. Implemented naively from PyTorch primitives, it requires five kernel launches and eight passes over the data:

def naive_softmax(x): # x: [rows, cols]  m = x.max(dim=1)[0] # kernel 1: read x  z = x - m[:, None] # kernel 2: read x, write z  e = torch.exp(z) # kernel 3: read z, write e  d = e.sum(dim=1) # kernel 4: read e  return e / d[:, None] # kernel 5: read e, write y 

If an entire row fits inside on-chip SRAM, one program block per row can load it once into registers, compute the maximum and the sum on the fly, and write the normalized probabilities back to HBM in a single pass. That cuts memory traffic from eight passes down to two — delivering an immediate ~4× speedup. Below is a simplified implementation adapted from Triton’s fused-softmax tutorial:

@triton.jit def softmax_kernel(x_ptr, y_ptr, cols, stride, BLOCK: tl.constexpr): row = tl.program_id(0) # one program per row  offs = tl.arange(0, BLOCK) # BLOCK = next power of 2 >= cols  mask = offs < cols x = tl.load(x_ptr + row * stride + offs, mask=mask, other=-float("inf")) x = x - tl.max(x, axis=0) # subtract the row max for stability  e = tl.exp(x) tl.store(y_ptr + row * stride + offs, e / tl.sum(e, axis=0), mask=mask) 

What happens when the sequence length is too large to fit in on-chip SRAM (like in long-context models)? You cannot compute the true row maximum before seeing the entire row — unless you use the online softmax trick (Milakov and Gimelshein, 2018). Online softmax maintains a running maximum $m$ and a running sum $d$ , dynamically rescaling previous partial sums whenever a larger element appears in the next chunk:

This streaming reduction allows the kernel to compute exact softmax statistics in a single streaming pass. Applied to attention, this exact rescaling lets a kernel process keys and values block-by-block without ever materializing the massive $N\times N$ attention matrix in HBM — the core algorithmic breakthrough behind