GEMM on NVIDIA GPUs: From Naive Kernels to Tiled Pipelines

Abstract

GEMM, short for General Matrix Multiplication, multiplies two matrices and optionally adds an existing matrix. Linear layers, attention projections, and the main computations in an MLP all map to GEMM. To make sense of deep-learning workload performance, it helps to know where GEMM gets fast—and where it gets stuck.

The first thing to examine is not the multiply-add itself, but data movement. A naive kernel assigns one output to each thread. In doing so, many threads fetch the same elements of A and B from global memory again and again. A tiled kernel turns those redundant transfers into reuse in shared memory and registers. As arithmetic intensity rises, the bottleneck moves away from HBM bandwidth and toward on-chip bandwidth, instruction throughput, and scheduling.

We will start with the mathematical definition of GEMM and the traffic generated by a naive implementation. From there, we will add shared-memory tiling, register blocking, and a Tensor Core pipeline, then move on to shape-aware scheduling and fusion. cuBLAS, cuBLASLt, CUTLASS, and CuTe come at the end, once there is a concrete reason for each layer to exist.

The API and architecture descriptions were checked against the official CUDA 13.3 and CUTLASS 4.6.1 documentation in August 2026. Supported dtypes, epilogues, schedulers, and instruction paths vary by toolkit and GPU generation, so confirm the documentation and API queries for the deployment environment before applying them.

1. GEMM Definition

The general GEMM computes

$$ D=\alpha\,\operatorname{op}(A)\operatorname{op}(B)+\beta C. $$

Here op means that a matrix is used as-is or transposed. Without transposition, the shapes are

$$ A\in\mathbb{R}^{M\times K},\qquad B\in\mathbb{R}^{K\times N},\qquad C,D\in\mathbb{R}^{M\times N}, $$

and one output element is

$$ D_{ij}=\alpha\sum_{k=0}^{K-1}A_{ik}B_{kj}+\beta C_{ij}. $$

There are $M\times N$ dot products, each of length $K$. Counting a multiply and an add as one FLOP each, the main matrix product requires approximately

$$ 2MNK\quad\text{FLOPs}. $$

This approximation excludes the operations that apply alpha and beta.

The property to notice here is reuse. One $A_{ik}$ contributes to $N$ outputs in the same row. One $B_{kj}$ contributes to $M$ outputs in the same column. If each input were fetched once and reused everywhere it is needed, the amount of computation per byte would grow with the matrix dimensions.

The equation may contain reuse, but a kernel does not get it for free. A fast GPU implementation must decide where those values live in the memory hierarchy and which threads share them.

2. Naive GEMM Bottlenecks

The obvious parallelization assigns one output element to one CUDA thread. The example below computes $D=AB$, with $\alpha=1$, $\beta=0$, and no transposition, and assumes row-major contiguous buffers for A, B, and D.

__global__ void naive_gemm(
    const float* A,
    const float* B,
    float* D,
    int M,
    int N,
    int K)
{
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row >= M || col >= N) {
        return;
    }

    float acc = 0.0f;
    for (int k = 0; k < K; ++k) {
        acc += A[row * K + k] * B[k * N + col];
    }
    D[row * N + col] = acc;
}

The kernel is correct and exposes plenty of output parallelism. At first glance, it looks reasonable. The problem is that every thread computes its dot product on its own.

Missing Input Reuse

Computing one output reads $K$ elements from A and $K$ elements from B. If the $MN$ outputs are all computed independently, the logical input-read count is

$$ 2MNK\quad\text{elements}. $$

With $b$ bytes per input element and $b_D$ bytes per output element, the traffic is approximately

$$ \text{bytes}_{\text{naive}} \approx 2bMNK+b_DMN, $$

and an arithmetic intensity of

$$ I_{\text{naive}} \approx \frac{2MNK}{2bMNK+b_DMN} \xrightarrow[K\to\infty]{} \frac{1}{b}. $$

This works out to only about 0.25 FLOP/byte for FP32 inputs and 0.5 FLOP/byte for FP16 or BF16 inputs.

That is not the same as literal DRAM traffic. Threads in a warp can receive an A element through a broadcast. Adjacent B elements can be coalesced, and the L1 and L2 caches absorb some duplicate requests. What the calculation shows is simpler: when every output issues its own loads, the kernel creates no explicit reuse of its own.

Cache hits and efficient transactions make actual HBM traffic lower than this estimate. Once the working set exceeds cache capacity, or many CTAs compete for that capacity, the same A and B tiles travel repeatedly from L2 and HBM. This does not mean that every naive GEMM reaches peak DRAM bandwidth. It does identify the likely bottleneck.

Naive GEMM does not own the reuse that exists in the mathematics. On large problems, repeated global-memory loads make operand delivery to the arithmetic units a likely bottleneck.

The mapping above is row-major. When a warp processes adjacent columns of one output row, its B accesses are coalesced. At a fixed k, all lanes request the same A address. Hardware broadcast cuts the transaction count, but the warp may use only a small part of the cache sector. There is no guarantee that the line needed at the next k will still be resident. Other output rows and CTAs reread the same B tile as well. NVIDIA's CUDA C++ Best Practices Guide walks through a matrix-multiplication example that removes these redundant transfers with shared memory.

GEMM Arithmetic Intensity

If A and B could each be read once and D written once, the idealized traffic for the whole operation with $\beta=0$ would be roughly

$$ b(MK+KN)+b_DMN. $$

When $\beta\ne0$, reading C adds $b_CMN$ bytes. As $M$, $N$, and $K$ grow together, FLOPs grow cubically while the input and output sizes grow quadratically. GEMM is capable of high arithmetic intensity. The naive mapping simply fails to use that advantage.

From a Roofline perspective,

$$ P_{\text{attainable}} \leq \min\left(P_{\text{compute peak}},\ I\cdot B_{\text{memory}}\right). $$

When $I$ is low, as in the naive implementation, the memory roof limits performance first. The FLOP count is not what needs to shrink. The next step is to raise $I$ by moving fewer global-memory bytes for the same work.

3. Tiled GEMM and Data Reuse

With tiling, one CTA computes a group of nearby outputs. The pieces of A and B needed for that work are loaded into shared memory once.

HBM / global memory
  cooperative load of A and B tiles
          ↓
Shared memory
  reuse by multiple threads in the block
          ↓
Registers
  accumulation of per-thread partial outputs
          ↓
HBM / global memory
  one store of the completed output

A simplified square-tile kernel looks like this.

template <int TILE>
__global__ void tiled_gemm(
    const float* A,
    const float* B,
    float* D,
    int M,
    int N,
    int K)
{
    __shared__ float As[TILE][TILE];
    __shared__ float Bs[TILE][TILE];

    int tx = threadIdx.x;
    int ty = threadIdx.y;
    int row = blockIdx.y * TILE + ty;
    int col = blockIdx.x * TILE + tx;
    float acc = 0.0f;

    for (int k0 = 0; k0 < K; k0 += TILE) {
        As[ty][tx] = (row < M && k0 + tx < K)
            ? A[row * K + k0 + tx] : 0.0f;
        Bs[ty][tx] = (k0 + ty < K && col < N)
            ? B[(k0 + ty) * N + col] : 0.0f;

        __syncthreads();
        for (int k = 0; k < TILE; ++k) {
            acc += As[ty][k] * Bs[k][tx];
        }
        __syncthreads();
    }

    if (row < M && col < N) {
        D[row * N + col] = acc;
    }
}

The first barrier prevents computation from starting before the tile is ready. The second keeps the next K tile from overwriting shared memory while the current computation is still using it. In production kernels, each thread may load several elements of a larger tile before asynchronous pipelines and register tiles are added.

Global-Memory Traffic Reduction

Suppose a CTA computes a $B_M\times B_N$ output tile and traverses the K dimension in units of $B_K$. One K tile requires

$$ \begin{aligned} \text{A tile} &: B_MB_K,\\ \text{B tile} &: B_KB_N,\\ \text{compute} &: 2B_MB_NB_K\ \text{FLOPs}. \end{aligned} $$

One A element feeds $B_N$ outputs in the CTA. One B element feeds $B_M$ outputs. If each input is read once from HBM, the input-only arithmetic intensity is

$$ I_{\text{tile,input}} \approx \frac{2B_MB_NB_K} {b(B_MB_K+B_KB_N)} = \frac{2B_MB_N}{b(B_M+B_N)}. $$

For a square tile with $B_M=B_N=T$,

$$ I_{\text{tile,input}}\approx\frac{T}{b}. $$

Doubling the tile edge doubles input reuse in this simple model. The CTA now manages the reuse that the naive mapping left to cache behavior. If cooperative loads follow contiguous addresses, the global accesses are coalesced as well.

On-Chip Bottlenecks

A larger tile is not automatically faster. It increases several costs at once.

• shared-memory capacity and bandwidth
• per-thread accumulator registers and register pressure
• block barriers and pipeline state
• inactive lanes and tail waste at problem boundaries
• lower residency caused by greater per-CTA resource use

The equation above counts input traffic only. It leaves out the D store, the read of $C$ when beta is nonzero, alignment and padding, and redundant loads between CTAs. Output and epilogue traffic take a larger share when K is short. Small M or N can also leave too few tiles to fill the GPU.

Tiling is not a complete cure for a memory-bound kernel. It removes one source of repeated traffic between HBM and shared memory. Once that traffic falls, the next limit becomes visible: shared-memory access, register reuse, compute instructions, or synchronization.

4. Hierarchical Tiling and the Tensor Core Pipeline

One CTA tile is not enough for a high-performance GEMM. The output and K reduction are divided again to match the CUDA execution hierarchy and the memory hierarchy.

Device GEMM
  ↓
CTA tile: global memory → shared memory, grid scheduling
  ↓
Warp / warp-group tile: partitioning of shared-memory tiles
  ↓
MMA instruction tile: Tensor Core or CUDA Core operation
  ↓
Register / architecture-specific state: operand and accumulator reuse

Register Blocking

After shared-memory tiling reduces HBM traffic, shared-memory bandwidth may show up as the next limit. Suppose each thread or warp owns a small output tile instead of one scalar. An A fragment fetched from shared memory can feed several column accumulators. Likewise, one B fragment can feed several row accumulators.

Register blocking increases the number of FMAs per shared-memory load and keeps partial sums close to the arithmetic units throughout the mainloop. More accumulators also mean more register pressure. The number of resident warps may fall, and forcing the register count too low can spill values into local memory, bringing back traffic that tiling had removed.

Tensor Core Pipeline

A Tensor Core instruction performs

$$ D_{\text{frag}}\leftarrow A_{\text{frag}}B_{\text{frag}}+D_{\text{frag}} $$

on small matrix tiles. That instruction alone is not a GEMM kernel. Global-memory loads, shared-memory layouts, synchronization, tail handling, the epilogue, and output stores all remain. The higher the peak Tensor Core throughput, the more important the tiling and pipeline that deliver operands on time.

The operand path varies by architecture.

• In the warp-level MMA path widely used from Volta through Ampere, shared-memory operands move into per-thread register fragments and accumulators remain in registers.
• Hopper WGMMA references B through a shared-memory descriptor. Depending on the configuration, A comes from shared memory or registers, while accumulators remain in registers. PTX ISA — WGMMA
• Blackwell SM100 tcgen05.mma stores accumulators in Tensor Memory (TMEM). A may come from shared memory or TMEM, and B from shared memory. NVIDIA CUTLASS — tcgen05 MMA Programming Guide

For that reason, a single “per-thread register fragment” model does not cover every generation.

Load–Compute Overlap

Tiling reduces the number of transferred bytes. Software pipelining hides the latency of the transfers that remain.

load tile 0
wait tile 0

compute tile 0  ||  load tile 1
compute tile 1  ||  load tile 2
compute tile 2  ||  load tile 3

Double buffering and multi-stage pipelines prefetch the next tiles. They also consume more shared memory and pipeline state. A pipeline that is too deep lowers occupancy; on short-K problems, it may add setup cost without enough steady-state work to recover it. Ampere asynchronous copies, TMA from Hopper onward, and generation-specific MMA and barriers take different forms, but all are used to overlap loads with computation.

Epilogue

When the mainloop finishes, the accumulators are written to global memory in the output layout. If the combination is supported, the epilogue handles alpha, beta, bias, activation, clamping, and dtype conversion in the same stage.

A supported epilogue removes the intermediate round trip in which GEMM stores a result and another kernel reads it back. The tradeoff changes if irregular indexing or a heavy transform slows the mainloop. In that case, a tuned GEMM followed by a small, separate kernel may be faster.

5. Bottlenecks After Tiling

No tile and pipeline works best for every GEMM shape. A large square GEMM has plenty of output tiles and a long steady-state mainloop. Small-M, short-K, and tail-heavy problems run into different limits.

Split-K divides the K range of one output tile among several workers.

$$ P^{(s)}_{ij}=\sum_{k\in K_s}A_{ik}B_{kj},\qquad P_{ij}=\sum_sP^{(s)}_{ij}. $$

K-direction parallelism increases, but the partial sums must be stored and combined. If the existing M/N tiles already fill the GPU, that reduction overhead makes the kernel slower. The Stream-K family spreads work more evenly and pays for scheduler and fix-up work in return.

In a small-M GEMM, the number of tiles along M falls. So does the opportunity to reuse the same weights B across output rows. Nominal TFLOPS from a large-GEMM tile says little here. Absolute latency, weight bytes, and active CTAs are more useful. Batching requests improves B reuse and parallelism, but adds queueing latency.

For tile tails, look at the remainder. If $M=130$ and the CTA tile size along M is 64, the third tile uses only two rows. In that case, boundary waste matters more than the peak throughput of the larger tile.

Record the input, compute or accumulator, and output dtypes separately. BF16 input, FP32 accumulation, and BF16 output make up one combination. FP8 and block-scaled formats also need a scale dtype and granularity. If Split-K, a scheduler, or an epilogue changes the reduction order or a rounding boundary, the same GEMM equation can produce different bits.

6. Abstraction Layers

So far, we have followed kernel optimization from the bottom up. Product implementation usually goes the other way. Start at the highest layer that can express the requirement.

Mathematical operation
  GEMM semantics and shapes
        ↓
Library / API
  cuBLAS → cuBLASLt
        ↓
Kernel construction
  CUTLASS device operators → CuTe components
        ↓
Custom implementation
  CUDA C++ / architecture-specific instructions
        ↓
Hardware
  scheduler, SM, Tensor Core, register, shared memory, HBM

cuBLASLt is not a higher version of cuBLAS. It is a separate API for describing layouts, algorithms, heuristics, and epilogues more flexibly. A kernel written with CUTLASS or CuTe is not automatically faster than either library.

To understand the mechanism, start with a naive kernel and work upward. To build the product, start with cuBLAS and move downward only as far as the requirement demands.

Shape Manifest

Collect the shapes from the production path before choosing an abstraction layer.

name, M, N, K, batch/group, transA, transB,
A/B/C/D layout and dtype, compute type,
alpha, beta, epilogue, alignment,
workspace limit, frequency or probability

An average shape can hide tails, small-M and short-K regimes, and heterogeneous groups. For a standard operation, use cuBLAS as the baseline and compare cuBLASLt heuristic candidates under the same conditions. Add CUTLASS or a custom kernel when the API cannot express the requirement, or when the same bottleneck keeps appearing on important shapes.

Keep the following information next to the timing results.

• For a small GEMM, prioritize median and tail absolute latency over TFLOPS.
• For a large GEMM, record latency and achieved FLOP/s, then inspect SM and Tensor Core utilization and DRAM and L2 traffic in Nsight Compute.
• For a custom tile, record register count, spills, shared-memory use, occupancy, and tail utilization.
• For a fused epilogue or quantized decode, measure the end-to-end interval around the call rather than pure GEMM alone.
• Move allocation and copies outside the measured interval. If production is not cache-hot, rotate through multiple buffers. NVIDIA CUTLASS — GEMM Performance Measurement Methodology

Define a reference and tolerance for the input, compute, and output dtypes as well. Byte-exact reproducibility may require fixing the toolkit, GPU architecture, algorithm, and workspace configuration together. If that level of reproducibility is unnecessary, use an error bound that is meaningful to the application. NVIDIA cuBLAS — Results Reproducibility

GEMM Optimization Criteria

Thread count is not the reason naive GEMM is slow. The problem is that each output calculation keeps moving the same operands. Shared-memory tiling reuses A and B within a CTA and cuts global-memory traffic. Register blocking applies the same idea at the warp and thread levels. Once that dataflow is in place, Tensor Cores and software pipelines can raise compute throughput and hide latency.

At that point, shape changes the answer. Small-M, short-K, tails, batching, quantization, and epilogues expose different bottlenecks and call for different schedulers. cuBLAS, cuBLASLt, CUTLASS, CuTe, and custom CUDA let us make those choices at different levels. There is no reason to drop to the lowest layer by default. Use the layer that exposes enough control to remove the bottleneck in front of you.

References

1. NVIDIA, cuBLAS 13.3 Documentation: https://docs.nvidia.com/cuda/cublas/
2. NVIDIA, cuBLASLt API: https://docs.nvidia.com/cuda/cublas/#using-the-cublaslt-api
3. NVIDIA, CUDA C++ Best Practices Guide — Shared Memory in Matrix Multiplication: https://docs.nvidia.com/cuda/cuda-c-best-practices-guide/#shared-memory-in-matrix-multiplication-c-ab
4. NVIDIA, CUTLASS 4.6.1 Documentation: https://docs.nvidia.com/cutlass/latest/overview.html
5. NVIDIA, Efficient GEMM in CUDA: https://docs.nvidia.com/cutlass/latest/media/docs/cpp/efficient_gemm.html
6. NVIDIA, CUTLASS GEMM API: https://docs.nvidia.com/cutlass/latest/media/docs/cpp/gemm_api.html
7. NVIDIA, GEMM Performance Measurement Methodology Guidelines: https://docs.nvidia.com/cutlass/latest/media/docs/cpp/gemm_performance_measurement_methodology_guidelines.html
8. NVIDIA, PTX ISA — WGMMA: https://docs.nvidia.com/cuda/parallel-thread-execution/#asynchronous-warpgroup-level-matrix-multiply-accumulate-instructions-wgmma
9. NVIDIA, tcgen05 MMA Programming Guide: https://docs.nvidia.com/cutlass/latest/media/docs/pythonDSL/mma_docs/tcgen05_programming.html
10. NVIDIA, Grouped Kernel Schedulers: https://docs.nvidia.com/cutlass/latest/media/docs/cpp/grouped_scheduler.html