AXPY, GEMV, GEMM — The Only Three Ops

BLAS Levels 1, 2, 3 — and which one earns its GPU keep

Once you have vectors and matrices, the operations BLAS organizes everything around are just three:

Op	Level	Signature	Kernel names	Typical use
AXPY	BLAS-1	`y ← α·x + y`	SAXPY, VectorAdd	Bias add, residual skip connection
GEMV	BLAS-2	`y ← A·x + y`	MatVec, DenseInference	Single-token autoregressive inference
GEMM	BLAS-3	`C ← α·A·B + β·C`	matmul, mma, wmma	Self-attention, batched inference, training

Why does BLAS Level 3 matter so much more than 1 and 2? Arithmetic intensity. AXPY does ~1 FLOP per byte loaded — every operation requires a fresh DRAM trip. GEMV is similar. GEMM, by contrast, uses each loaded element O(k) times, so it has the FLOPs/byte ratio that GPUs are built to crush.

Roofline rule of thumb:

Low FLOPs/byte (< 4) — you're bandwidth-bound. Your kernel waits on DRAM. AXPY, GEMV, element-wise ops live here.
High FLOPs/byte (> 8) — you're compute-bound. ALUs / Tensor Cores are saturated. Tiled GEMM lives here.

This is why a single-token decode (GEMV-shaped: activation × weight) on a 7B model runs at like 50 tokens/s on an RTX 4090, while a batch of 64 tokens (GEMM-shaped) runs at 800. The math is the same, the FLOPs/byte ratio is what changed.

Code

The three ops, side by side, in NumPy·python

import numpy as np

# AXPY — vector scale + add (BLAS-1)
alpha = 0.5
x = np.random.randn(1_000_000)
y = np.random.randn(1_000_000)
y = alpha * x + y          # one FLOP per element ≈ 1 FLOP/byte

# GEMV — matrix × vector (BLAS-2)
A = np.random.randn(4096, 4096)
x = np.random.randn(4096)
y = A @ x                  # ~ 2N² FLOPs over ~N² bytes ≈ 2 FLOPs/byte

# GEMM — matrix × matrix (BLAS-3)
A = np.random.randn(4096, 4096)
B = np.random.randn(4096, 4096)
C = A @ B                  # ~ 2N³ FLOPs over ~N² bytes ≈ N FLOPs/byte

Arithmetic intensity calculator (back-of-envelope)·python

def intensity(flops, bytes_):
    return flops / bytes_

# 4096 × 4096 GEMM, FP32 (4 bytes/element)
N = 4096
flops = 2 * N**3
bytes_ = 3 * N * N * 4         # A, B, C each N² elements
print(f'GEMM intensity: {intensity(flops, bytes_):.0f} FLOPs/byte')
# ~ 1365  → very compute-bound, GPUs love this

# Same shape but vector add
flops = N * N
bytes_ = 3 * N * N * 4
print(f'AXPY intensity: {intensity(flops, bytes_):.2f} FLOPs/byte')
# ~ 0.08 → severely bandwidth-bound

Exercise

Run the intensity calculator above for 1024×1024 and 8192×8192 GEMM. Notice that intensity grows linearly with N (because compute is N³ and bytes are N²). This is the algebraic reason large GEMMs are easier to make compute-bound than small ones — and why padding small matrices up to a 'good' size sometimes wins despite the extra zero work.

AXPY, GEMV, GEMM — The Only Three Ops

BLAS Levels 1, 2, 3 — and which one earns its GPU keep

Code

External links

Exercise

Progress

Comments 0