Software Stack & Performance — From CUDA to LLM Throughput

Topics We'll Cover

The capstone of the architecture series. Up to this point we have inspected the silicon — SMs, tensor cores, memory hierarchy, NVLink, every architecture from Ampere to Blackwell. None of that runs an LLM by itself. The software stack is what turns a transistor count into a tok/s number.

• The CUDA Stack — Layered Cake
• cuBLAS / cuBLASLt — Matrix Multiplication
• cuDNN — Deep Neural Networks
• CUTLASS — Templates for Custom Matmul
• Transformer Engine — Auto FP8/FP4
• NCCL — Collective Comms
• Triton — Python-First Kernel DSL
• TensorRT-LLM + vLLM + SGLang — The Inference Layer
• Performance Model — Decode-Bound vs Prefill-Bound
• Tensor Parallel Scaling — How Architectures Differ
• Interactive: End-to-End LLM Performance Model

The CUDA Stack — Layered Cake

Every call from model.generate(...) down to a tensor-core instruction passes through six layers. Knowing where each layer sits is the difference between debugging a NCCL timeout and hunting a kernel-mode firmware bug.

CUDA versions vs minimum hardware

cuBLAS / cuBLASLt — Matrix Multiplication

The whole transformer is, fundamentally, a stack of matrix multiplications interleaved with elementwise ops. Whatever framework you use, almost every FLOP eventually lands in cuBLAS or its successor cuBLASLt.

FP8 matmul on Hopper, with FP32 accumulator

// Descriptors set: A=E4M3 row, B=E4M3 col, scaleType=FP32, accumType=FP32, output=BF16
cublasLtMatmulDesc_t op;
cublasLtMatmulDescCreate(&op, CUBLAS_COMPUTE_32F, CUDA_R_32F);
cublasLtMatmulDescSetAttribute(op, CUBLASLT_MATMUL_DESC_A_SCALE_POINTER, &a_scale, sizeof(void*));
cublasLtMatmulDescSetAttribute(op, CUBLASLT_MATMUL_DESC_B_SCALE_POINTER, &b_scale, sizeof(void*));

// Heuristic asks the library what kernel + workspace fit this problem
cublasLtMatmulHeuristicResult_t heur[8];
int n;
cublasLtMatmulAlgoGetHeuristic(handle, op, ad, bd, cd, dd, pref, 8, heur, &n);

// Three-line execution: A * B + bias, with chosen algo and workspace
cublasLtMatmul(handle, op,
    &alpha, A, ad, B, bd, &beta, C, cd, D, dd,
    &heur[0].algo, workspace, ws_bytes, stream);

Where cuBLASLt earns its keep

• Mixed precision in one call. FP8 weights, BF16 activations, FP32 accumulator, BF16 output — all expressed in descriptors, no manual upcasting.
• Fused epilogues. Bias add, GELU, ReLU, dropout mask all run in the same kernel, no extra HBM round-trip.
• Per-tensor or per-channel scales. The scaling factors that FP8 / MX-FP4 need are first-class arguments.
• Automatic algo selection. No more "I tuned this kernel for B=32 and now B=33 is 4× slower."

cuDNN — Deep Neural Networks

cuDNN is the layer above cuBLAS that knows what a neural network is rather than just a matrix. Convolutions, attention, normalization, activations, pooling, RNN cells — the abstraction is "an op fused into a graph", not "a single GEMM".

cuDNN 9's flash-attention-style attention

# Python via cudnn-frontend bindings — simplified
import cudnn

graph = cudnn.pygraph(io_data_type=cudnn.data_type.HALF)

q = graph.tensor_like(Q)   # [B, H, S, D]
k = graph.tensor_like(K)
v = graph.tensor_like(V)

o, stats = graph.sdpa(
    name="sdpa", q=q, k=k, v=v,
    is_inference=True,
    attn_scale=1.0/math.sqrt(D),
    use_causal_mask=True,
    use_alibi=False,
)
o.set_output(True).set_dim([B, H, S, D])

graph.build([cudnn.heur_mode.A])  # picks best kernel for shape
graph.execute({q: Q, k: K, v: V, o: O}, workspace)

CUTLASS — Templates for Custom Matmul

cuBLAS is a black box; you call it, it returns a result. CUTLASS is the opposite: a header-only C++ template library that exposes every level of the matmul hierarchy — tile shape, warp tile, instruction tile, pipeline depth, epilogue — for you to compose. NVIDIA writes its own kernels with it, and so does every research lab.

When canned cuBLAS isn't right

• Custom epilogue. "Multiply, add bias, then RoPE-rotate the output." cuBLASLt has many fused epilogues but not every one.
• Sparsity layouts. 2:4 structured sparsity, block-sparse, MoE-style block-diagonal — CUTLASS exposes the iterator primitives.
• Quantization shapes. INT4 with GPTQ groupwise scales, AWQ packed layouts, NVFP4 / MX-FP4 microblock metadata.
• Research kernels. Streamed attention variants, mixture-of-experts grouped GEMM, distributed weight-streaming matmul.

// Pick the architecture-specific atoms
using ArchTag    = cutlass::arch::Sm90;
using ElementA   = cutlass::float_e4m3_t;   // FP8 E4M3
using ElementB   = cutlass::float_e4m3_t;
using ElementAcc = float;                  // FP32 accumulate
using ElementC   = cutlass::bfloat16_t;     // BF16 output

using TileShape = Shape<_128, _128, _64>;  // CTA-level tile
using ClusterShape = Shape<_2, _1, _1>;    // thread-block cluster

// Compose mainloop + epilogue
using Mainloop = cutlass::gemm::collective::CollectiveBuilder<
    ArchTag, OperatorClass, ElementA, LayoutA, AlignA,
    ElementB, LayoutB, AlignB, ElementAcc,
    TileShape, ClusterShape,
    cutlass::gemm::collective::StageCountAuto,
    cutlass::gemm::KernelTmaWarpSpecializedPingpong >::CollectiveOp;

Transformer Engine — Auto FP8/FP4

Decks 03 and 07 covered the silicon: per-tensor FP8 scales on Hopper, per-microblock FP4 (NVFP4 and MX-FP4) on Blackwell. Deck 09 covered the second-gen Transformer Engine on Blackwell. Transformer Engine the library is the piece that drives the silicon.

The Python integration

import transformer_engine.pytorch as te
from transformer_engine.common import recipe

# A transformer block, with TE-aware Linear and LayerNorm
class Block(nn.Module):
    def __init__(self, d, hidden):
        super().__init__()
        self.norm = te.RMSNorm(d)
        self.fc1  = te.Linear(d, hidden)        # <-- not torch.nn.Linear
        self.fc2  = te.Linear(hidden, d)

# Per-tensor scaling recipe (Hopper); MX recipe on Blackwell
fp8_recipe = recipe.DelayedScaling(
    margin=0, fp8_format=recipe.Format.HYBRID,
    amax_history_len=16, amax_compute_algo="max")

# Forward inside the autocast block runs FP8 in tensor cores
with te.fp8_autocast(enabled=True, fp8_recipe=fp8_recipe):
    y = block(x)

NCCL — Collective Comms

Once a model crosses one GPU, every step gathers, all-reduces, or reduce-scatters tensors across ranks. NCCL (NVIDIA Collective Communications Library) is the implementation. It is in every distributed framework: PyTorch DDP, FSDP, vLLM TP, Megatron, JAX. Bad NCCL configuration is the most common cause of "8 GPUs do less than 4".

Collectives that matter for transformers

Algorithms NCCL picks between

Critical environment knobs

# Pick the algorithm explicitly — default heuristics are usually right but not always
export NCCL_ALGO=Ring                     # or Tree, or NVLS, or NVLSTree
export NCCL_PROTO=Simple                  # Simple | LL | LL128

# Enable in-network reductions on NVSwitch (Hopper/Blackwell HGX, NVL72)
export NCCL_NVLS_ENABLE=1

# Diagnostics
export NCCL_DEBUG=INFO
export NCCL_TOPO_DUMP_FILE=/tmp/topo.xml

# Turn off PXN if you are debugging weird stalls; turn it back on for perf
export NCCL_PXN_DISABLE=0

# Force IB over a particular HCA, exclude management interfaces
export NCCL_IB_HCA=mlx5_0,mlx5_1
export NCCL_SOCKET_IFNAME=^lo,docker0

Triton — Python-First Kernel DSL

Until Triton, writing a custom CUDA kernel meant C++, CUTLASS templates, and very specific knowledge of warp scheduling. Triton (originally OpenAI) provides a Python-like language where you express tile-level ops on tensors, and a compiler handles vectorisation, shared-memory tiling, and warp-level scheduling automatically.

A Triton matmul, end to end

import torch, triton
import triton.language as tl

@triton.jit
def matmul_kernel(A, B, C, M, N, K,
                  sa0, sa1, sb0, sb1, sc0, sc1,
                  BM: tl.constexpr, BN: tl.constexpr, BK: tl.constexpr):
    pid_m = tl.program_id(0)
    pid_n = tl.program_id(1)
    rm = pid_m * BM + tl.arange(0, BM)
    rn = pid_n * BN + tl.arange(0, BN)
    rk = tl.arange(0, BK)

    a_ptr = A + rm[:, None] * sa0 + rk[None, :] * sa1
    b_ptr = B + rk[:, None] * sb0 + rn[None, :] * sb1
    acc   = tl.zeros([BM, BN], dtype=tl.float32)

    for k in range(0, K, BK):
        a = tl.load(a_ptr, mask=(rm[:, None] < M) & (rk[None, :] + k < K))
        b = tl.load(b_ptr, mask=(rk[:, None] + k < K) & (rn[None, :] < N))
        acc += tl.dot(a, b)            # maps to WMMA / WGMMA / Blackwell MMA
        a_ptr += BK * sa1
        b_ptr += BK * sb0

    c_ptr = C + rm[:, None] * sc0 + rn[None, :] * sc1
    tl.store(c_ptr, acc.to(tl.float16),
             mask=(rm[:, None] < M) & (rn[None, :] < N))

def matmul(a, b, BM=128, BN=128, BK=32):
    M, K = a.shape; _, N = b.shape
    c = torch.empty((M, N), device=a.device, dtype=torch.float16)
    grid = (triton.cdiv(M, BM), triton.cdiv(N, BN))
    matmul_kernel[grid](a, b, c, M, N, K,
                        a.stride(0), a.stride(1),
                        b.stride(0), b.stride(1),
                        c.stride(0), c.stride(1),
                        BM=BM, BN=BN, BK=BK)
    return c

TensorRT-LLM + vLLM + SGLang — The Inference Layer

Above all the kernels sits the inference engine — the thing that schedules requests, manages KV cache, batches prompts, and turns a model checkpoint into an HTTP-served chat completion. Three engines dominate in 2026.

What TensorRT-LLM's engine actually contains

• Model graph — the transformer layers as a fixed DAG, partitioned for TP / PP.
• Pre-selected kernels — one CUTLASS / cuBLASLt / Flash-Attn kernel per layer, tuned for the chosen GPU and batch profile.
• KV-cache manager — paged blocks, with sizes baked in.
• Scheduler — in-flight batching of prefill + decode at token granularity.
• Plugins — custom ops the graph compiler couldn't fuse on its own (e.g. attention).

Performance Model — Decode-Bound vs Prefill-Bound

An LLM serving workload has two distinct phases, with different bottlenecks. Get this right and the rest of capacity planning falls into place.

The KV-cache cost

Every new decode token writes a fresh row into the K and V caches:

kv_per_token = 2 × num_layers × num_kv_heads × head_dim × bytes

Tensor Parallel Scaling — How Architectures Differ

Tensor parallel splits each weight matrix across N GPUs and inserts an all-reduce after every attention output and every MLP output. Whether N× throughput is realised depends entirely on how fast that all-reduce is relative to the GEMM it follows.

The compute / bandwidth ratio

The shape of TP scaling is not "NVLink = good, PCIe = bad" in isolation. It is the ratio of (GEMM time) to (all-reduce time). FP8 / FP4 halve or quarter the GEMM time without changing the all-reduce. So as you adopt lower precisions, the relative cost of comms goes up, and the gap between an NVLink box and a PCIe box widens. Hopper and Blackwell only really shine when the link beneath them is fast enough to feed the new tensor cores.

Interactive: End-to-End LLM Performance Model

Combine everything from the earlier slides into a single calculator. Pick a model, a precision, a GPU, a batch size and a sequence length, and get a quick verdict on whether you are decode-bound, prefill-bound, or out of VRAM.