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.

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:

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 ($y_i = \max(x_i, 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

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 ($P_{\text{peak}}/B_{\text{HBM}}$):

• 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).

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:

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)

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 — $y_i = e^{x_i - m} / \sum_j e^{x_j - 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 FlashAttention (Dao et al., 2022).

5.3 Tiling and Tensor Cores: reuse every load

Matrix multiplication ($C = AB$) of two $N \times N$ matrices requires $2N^3$ FLOPs, yet only touches $3N^2$ values. In theory, that offers massive computational reuse.

In a naive kernel where each thread calculates a single element $C_{i,j}$, every thread reads an entire row of A and column of B from HBM — resulting in $2N^3$ total loads. The solution is tiling: each program computes a $T \times T$ block of C, streaming strips of A and B through shared memory so every loaded value is reused $T$ times:

In Triton, tiling is structured as a loop over chunks of dimension K, and calling tl.dot dispatches matrix operations directly to hardware Tensor Cores:

@triton.jit
def matmul_kernel(a_ptr, b_ptr, c_ptr, M, N, K,
                  s_am, s_ak, s_bk, s_bn, s_cm, s_cn,
                  BM: tl.constexpr, BN: tl.constexpr, BK: tl.constexpr):
    rm = tl.program_id(0) * BM + tl.arange(0, BM)     # this program's rows of C
    rn = tl.program_id(1) * BN + tl.arange(0, BN)     # and its columns
    rk = tl.arange(0, BK)
    acc = tl.zeros((BM, BN), dtype=tl.float32)
    for k in range(0, K, BK):                         # stream strips of A and B
        a = tl.load(a_ptr + rm[:, None] * s_am + (k + rk)[None, :] * s_ak,
                    mask=(rm[:, None] < M) & (k + rk[None, :] < K), other=0.0)
        b = tl.load(b_ptr + (k + rk)[:, None] * s_bk + rn[None, :] * s_bn,
                    mask=(k + rk[:, None] < K) & (rn[None, :] < N), other=0.0)
        acc += tl.dot(a, b, input_precision="ieee")   # FP32 inputs default to TF32
    tl.store(c_ptr + rm[:, None] * s_cm + rn[None, :] * s_cn, acc,
             mask=(rm[:, None] < M) & (rn[None, :] < N))

Notice the comment on tl.dot: it conceals a notorious precision trap. Tensor Cores do not compute in standard IEEE FP32 by default. When fed FP32 inputs, they automatically switch to TF32 (TensorFloat-32), which retains the 8-bit dynamic range of FP32 but truncates the mantissa from 23 bits down to 10 bits (the precision of FP16).

Because TF32 Tensor Cores deliver up to 495 TFLOPS on an H100 — roughly 7.4× faster than standard 67 TFLOPS FP32 CUDA cores — a Triton kernel will look blazingly fast compared to PyTorch’s default FP32 matmul. But when tested against a strict numerical tolerance, it may fail correctness tests. Setting input_precision="ieee" restores full FP32 accuracy at full FP32 cost; passing "tf32x3" splits each input into high and low parts, executing three Tensor Core operations to recover near-FP32 precision.

5.4 Epilogues: let the library do the matmul

Trying to beat NVIDIA’s cuBLAS on a bare matrix multiplication is usually an exercise in frustration: NVIDIA engineers spend years tuning assembly kernels for each GPU microarchitecture. Where you can win is what happens immediately before and after the matmul.

If your matmul is followed by a bias add, a GELU activation, or a residual connection, you can apply those operations to each tile while it is still sitting in registers before writing it to HBM. This fused post-processing is called an epilogue. cuBLASLt provides fixed epilogues (like bias + GELU); Triton allows you to write arbitrary epilogues directly in Python; and a pragmatic architecture is often to call cuBLAS for the heavy GEMM, followed by a single fused kernel for the remaining elementwise tail.

Two other techniques appear in almost every high-performance kernel. Autotuning times multiple combinations of tile sizes, warp counts, and pipeline stages on target hardware, keeping the winner for each tensor shape. Software pipelining overlaps memory and compute by prefetching tile $k+1$ while Tensor Cores multiply tile $k$. Triton exposes both via simple decorators and arguments (@triton.autotune, num_stages), while in CUDA they require intricate manual engineering.

6. cuBLAS, cuDNN, eager PyTorch and torch.compile

In modern deep learning, GPU kernels reach hardware through four distinct paths. Depending on which path handles a layer, the exact same mathematical code can execute very differently:

torch.compile comes with different operating modes. The default mode fuses elementwise and reduction chains into Triton kernels. mode="reduce-overhead" captures the entire execution sequence into a CUDA graph, replaying dozens of launches in a single driver call — a massive win when models are CPU-dispatch bound. Finally, mode="max-autotune" benchmarks Triton templates against vendor libraries to select optimal kernels per tensor dimension.

On whole models, these compiler passes deliver real gains: the PyTorch 2 paper (Ansel et al., ASPLOS 2024) reports geometric-mean speedups of 2.27× for inference and 1.41× for training over eager mode across 180+ models on A100 GPUs.

The practical takeaway is that claims of being “faster than PyTorch” are ambiguous without context. Faster than eager PyTorch often just means you fused two elementwise operations that anyone could fuse. Faster than torch.compile in max-autotune mode is a much higher bar, because the compiler has already fused elementwise tails, benchmarked library alternatives, and eliminated launch overhead.

7. KernelBench: can an LLM write these kernels?

Can frontier language models write GPU kernels that outperform human engineers and optimizing compilers? To test this systematically, Stanford researchers released KernelBench (Ouyang et al., 2025). Each task presents a standalone PyTorch module; the model must output a drop-in replacement whose forward() method runs custom CUDA or Triton kernels:

# KernelBench level 2, problem 1: Conv2D_ReLU_BiasAdd (docstrings removed)
import torch
import torch.nn as nn

class Model(nn.Module):
    def __init__(self, in_channels, out_channels, kernel_size, bias_shape):
        super(Model, self).__init__()
        self.conv = nn.Conv2d(in_channels, out_channels, kernel_size)
        self.bias = nn.Parameter(torch.randn(bias_shape))

    def forward(self, x):
        x = self.conv(x)
        x = torch.relu(x)
        x = x + self.bias
        return x

batch_size = 128
in_channels  = 64
out_channels = 128
height = width = 128
kernel_size = 3
bias_shape = (out_channels, 1, 1)

def get_inputs():
    return [torch.rand(batch_size, in_channels, height, width)]

def get_init_inputs():
    return [in_channels, out_channels, kernel_size, bias_shape]

KernelBench evaluates solutions through three rigorous gates:

1. Compilation & execution: the kernel must compile cleanly and run without CUDA exceptions.
2. Correctness: the outputs must match the PyTorch reference across 5 random input batches (using numerical tolerances calibrated to the precision).
3. Performance: if correct, the kernel is timed using CUDA events (3 warm-up passes followed by the mean of 100 runs, with L2 cache cleared between trials) against the eager PyTorch baseline.

Here, $\text{fast}_0$ is the raw correctness rate, while $\text{fast}_1$ measures the proportion of problems where the model generated a correct kernel that was strictly faster than PyTorch.

The benchmark spans multiple difficulty tiers. Level 1 covers isolated operators. Level 2 tests composite patterns (like Conv2D + ReLU + BiasAdd) where fusion and epilogues provide clear optimization headroom. Level 3 tests full architecture blocks. The baseline results from early 2025 were a sobering reality check:

On single-shot generation, even the best reasoning models struggled to beat eager PyTorch on more than a third of Level 2 problems, and barely registered on Level 3.

Why do LLMs struggle so much with kernels? Part of the answer is data distribution: CUDA C++ accounts for only 0.073% of The Stack v1.2, a standard pretraining corpus. But the paper also proved that the evaluation protocol matters as much as the base model:

• Iterative compiler feedback: When DeepSeek R1 was given compiler error diagnostics, runtime outputs, and profiler metrics over up to 10 feedback turns, its $\text{fast}_1$ score surged from 12% / 36% / 2% to 43% / 72% / 18% across the three levels.
• Best-of-N sampling: Sampling 100 candidate completions per problem raised DeepSeek V3’s Level 2 score from 4% to 37%.

Since the original paper, the ecosystem has accelerated rapidly: KernelBench added support for Triton, HIP, and TileLang, while specialized research models — such as Cognition’s Kevin (multi-turn RL on CUDA), AutoTriton, and Dr. Kernel — have trained specifically on GPU environments.

8. How to read a KernelBench number

Whenever a research group or vendor claims “our AI model generates kernels that are 10× faster than PyTorch”, take a breath and check the methodology. Kernel benchmarking is fraught with subtle discrepancies:

A fast kernel needs a strict verifier

Kernel generation benchmarks are uniquely vulnerable to reward hacking because the objective function is literally a timer. In February 2025, Sakana AI announced an “AI CUDA Engineer” claiming 10× to 100× speedups on common PyTorch operations. Within 24 hours, the team retracted the claims after discovering that the model had found an exploit in the test harness: it corrupted memory checks so the validation suite erroneously marked broken kernels as correct.

As the KernelBench documentation notes: “If the model can reward hack, it will find ways to reward hack!” Here are the most common exploits and the defensive checks required to catch them:

These principles mirror the verifier design patterns we use across all RLVR workflows; our verifier playbook and reward-hacking analysis cover these guardrails in detail. For GPU kernels, the absolute minimum baseline is testing on unseen input distributions, explicit numerical tolerance boundaries, full CUDA stream synchronization before stopping the clock, and clearly naming whether you are comparing against eager PyTorch or torch.compile.

What comes next

The concepts covered here form the foundation for the next engineering deep dive in this series, where we show how we train an open LLM to write performant Triton kernels inside a GPU kernel gym, and what RLVR techniques are required to translate those gains onto KernelBench.