Lab · runnable experiments
How GPUs Actually Run Machine Learning
Created Sep 10, 2026 Updated Sep 11, 2026
Read the parent noteThe note asks one question: why is the same GPU sometimes a hundred times faster than a CPU and sometimes barely faster at all? This lab answers it with measurements on one real card, and every term the note introduces — warps, SIMT, latency hiding, the memory hierarchy, FLOPs, arithmetic intensity, the roofline — is pinned to a number measured here.
Part 1 reproduces the paradox: matrix multiplication and element-wise addition swept from tiny to huge, on the CPU and on the GPU, with and without the copy over PCIe. Part 2 opens the programming model with hand-written CUDA kernels: what a warp is, and what it costs when the 32 threads of a warp disagree. Part 3 is memory: access patterns, the latency of each level of the hierarchy, how the GPU hides that latency, and the bandwidth of each level. Part 4 is compute: the FLOP/s the card delivers against what the datasheet promises, and what the tensor cores add. Part 5 puts both roofs on one chart and walks a single kernel from memory-bound to compute-bound, then places real PyTorch operations on it. Part 6 closes the loop: a three-number model of each machine, asked to predict the speed-ups from Part 1.
Everything runs on one Kaggle Tesla T4. The code is the code that ran; the outputs are the outputs it printed.
Setup
Run on Kaggle with Settings → Accelerator → GPU T4 x2 (the lab uses one of the two cards) and Internet on — CuPy is installed only if the image does not already have it.
import os
os.environ.setdefault("CUDA_DEVICE_ORDER", "PCI_BUS_ID") # torch, CuPy and nvidia-smi agree on "GPU 0"
import sys, json, time, datetime, platform, subprocess, threading, traceback, contextlib
import numpy as np, pandas as pd, torch
import matplotlib.pyplot as plt
assert torch.cuda.is_available(), "No GPU: Settings → Accelerator → GPU T4 x2"
DEV = torch.device("cuda:0")
torch.cuda.set_device(DEV)
torch.backends.cuda.matmul.allow_tf32 = False # keep fp32 honest (Turing has no TF32 anyway)
torch.backends.cudnn.allow_tf32 = False
try:
import cupy as cp
except ImportError:
subprocess.run([sys.executable, "-m", "pip", "install", "-q", "cupy-cuda12x"], check=True)
import cupy as cp
cp.cuda.Device(0).use()
# CuPy kernels and PyTorch ops are timed with the same CUDA events, which is only valid if
# both libraries issue work to the same (legacy default) stream.
assert cp.cuda.get_current_stream().ptr == 0 and torch.cuda.current_stream().cuda_stream == 0
INK, LINE, CREAM = "#2a2a2a", "#d9d3c7", "#f6f2e9"
BLUE, EMBER, FOREST, GRAY, PLUM = "#3b6ea5", "#c4521e", "#2e7d5b", "#9a9384", "#7a4f8c"
plt.rcParams.update({"font.size": 9, "axes.edgecolor": LINE, "axes.spines.top": False,
"axes.spines.right": False, "figure.dpi": 120})
RESULTS = {"meta": {}, "errors": {}}
@contextlib.contextmanager
def experiment(name):
"""One experiment per cell. A failure is printed and recorded, and the run moves on."""
t0 = time.time()
try:
yield
except Exception:
RESULTS["errors"][name] = traceback.format_exc()
print(f"!! {name} failed:\n{traceback.format_exc()}")
finally:
print(f"[{name}: {time.time() - t0:.1f} s]")
print("executed:", datetime.date.today().isoformat())
print("torch", torch.__version__, "| CUDA", torch.version.cuda, "| CuPy", cp.__version__,
"| Python", platform.python_version())executed: 2026-09-11
torch 2.10.0+cu128 | CUDA 12.8 | CuPy 14.0.1 | Python 3.12.13Three stopwatches. gpu_time brackets the work with CUDA events, so it measures what the GPU did, not how long Python took to ask. wall_time is an ordinary clock, for the CPU and for anything that includes copies. sustained keeps the GPU busy for a couple of seconds while a background thread polls nvidia-smi, because a card’s clock under load is not the clock on its datasheet. All three report the median of repeated runs.
def gpu_time(fn, iters=15, warmup=3, budget=1.0):
"""Median seconds of fn(), measured on the GPU with CUDA events."""
for _ in range(warmup):
fn()
torch.cuda.synchronize()
s, e = torch.cuda.Event(enable_timing=True), torch.cuda.Event(enable_timing=True)
s.record(); fn(); e.record(); e.synchronize()
n = int(min(iters, max(3, budget / max(s.elapsed_time(e) / 1e3, 1e-6))))
ts = []
for _ in range(n):
s.record(); fn(); e.record(); e.synchronize()
ts.append(s.elapsed_time(e) / 1e3)
return float(np.median(ts))
def wall_time(fn, iters=15, warmup=1, budget=1.5, sync=False):
"""Median wall-clock seconds of fn(); sync=True waits for the GPU to finish."""
def run():
if sync:
torch.cuda.synchronize()
t0 = time.perf_counter()
fn()
if sync:
torch.cuda.synchronize()
return time.perf_counter() - t0
for _ in range(warmup):
run()
n = int(min(iters, max(3, budget / max(run(), 1e-9))))
return float(np.median([run() for _ in range(n)]))
def smi(fields="clocks.sm,clocks.mem,power.draw,temperature.gpu"):
out = subprocess.run(["nvidia-smi", "-i", "0", f"--query-gpu={fields}",
"--format=csv,noheader,nounits"], capture_output=True, text=True).stdout
res = {}
for key, val in zip(fields.split(","), (v.strip() for v in out.strip().split(","))):
try:
res[key] = float(val)
except ValueError:
res[key] = val
return res
class Sampler:
"""Polls nvidia-smi in a background thread while a workload runs."""
def __init__(self, period=0.2):
self.period, self.samples, self._stop = period, [], threading.Event()
def _loop(self):
while not self._stop.is_set():
self.samples.append(smi())
self._stop.wait(self.period)
def __enter__(self):
self._thread = threading.Thread(target=self._loop, daemon=True)
self._thread.start()
return self
def __exit__(self, *exc):
self._stop.set()
self._thread.join()
def summary(self):
df = pd.DataFrame(self.samples[1:] or self.samples) # first sample can be pre-load
return {k: float(df[k].median()) for k in df.columns if pd.api.types.is_numeric_dtype(df[k])}
def sustained(fn, seconds=2.0):
"""Seconds per call of fn() under continuous load, plus the clocks and power it ran at."""
fn(); torch.cuda.synchronize()
s, e = torch.cuda.Event(enable_timing=True), torch.cuda.Event(enable_timing=True)
with Sampler() as smp:
s.record()
n, t0 = 0, time.perf_counter()
while time.perf_counter() - t0 < seconds:
fn(); n += 1
if n % 4 == 0:
torch.cuda.synchronize() # keep the queue from running far ahead of the clock
e.record(); e.synchronize()
return s.elapsed_time(e) / 1e3 / n, smp.summary()
def jsonable(o):
if isinstance(o, np.bool_):
return bool(o)
if isinstance(o, (np.integer,)):
return int(o)
if isinstance(o, (np.floating,)):
return float(o)
if isinstance(o, np.ndarray):
return o.tolist()
return str(o)The machine, as it describes itself
Before measuring anything, the card’s own description: how many streaming multiprocessors (SMs), how many FP32 lanes each has, how big the register file, shared memory and L2 are, and the memory bus. From those, two theoretical ceilings: FP32 FLOP/s = SMs × lanes × 2 × clock (a fused multiply-add counts as two FLOPs), and DRAM bandwidth = memory clock × 2 × bus width. Everything later is measured against these.
props = cp.cuda.runtime.getDeviceProperties(0)
tprops = torch.cuda.get_device_properties(0)
def prop(key, torch_attr=None):
v = props.get(key)
if v is None and torch_attr:
v = getattr(tprops, torch_attr, None)
return v.decode() if isinstance(v, (bytes, bytearray)) else v
cc = (prop("major", "major"), prop("minor", "minor"))
FP32_LANES_PER_SM = {(7, 0): 64, (7, 5): 64, (8, 0): 64, (8, 6): 128, (8, 7): 128, (8, 9): 128,
(9, 0): 128, (10, 0): 128, (12, 0): 128}
SPEC_FP16_TENSOR = {"Tesla T4": 65e12} # vendor datasheet, dense FP16 tensor-core peak
stat = smi("name,clocks.max.sm,clocks.max.mem,power.limit,pcie.link.gen.max,"
"pcie.link.width.max,driver_version")
sms = prop("multiProcessorCount", "multi_processor_count")
mem_mhz = (prop("memoryClockRate") or 0) / 1e3 or stat["clocks.max.mem"]
bus_bits = prop("memoryBusWidth")
spec = {
"gpu": prop("name", "name"),
"compute_capability": f"{cc[0]}.{cc[1]}",
"sms": sms,
"fp32_lanes_per_sm": FP32_LANES_PER_SM.get(cc, 64),
"warp_size": prop("warpSize", "warp_size"),
"max_threads_per_sm": prop("maxThreadsPerMultiProcessor", "max_threads_per_multi_processor"),
"regs_per_sm": prop("regsPerMultiprocessor", "regs_per_multiprocessor"),
"shared_mem_per_sm_bytes": prop("sharedMemPerMultiprocessor", "shared_memory_per_multiprocessor"),
"l2_bytes": prop("l2CacheSize", "L2_cache_size"),
"dram_bytes": prop("totalGlobalMem", "total_memory"),
"mem_bus_bits": bus_bits,
"mem_clock_mhz": mem_mhz,
"max_sm_clock_mhz": stat["clocks.max.sm"],
"power_limit_w": stat["power.limit"],
"pcie": f"gen {stat['pcie.link.gen.max']:.0f} x{stat['pcie.link.width.max']:.0f}",
"driver": stat["driver_version"],
}
spec["fp32_lanes_total"] = spec["sms"] * spec["fp32_lanes_per_sm"]
spec["max_warps_per_sm"] = spec["max_threads_per_sm"] // spec["warp_size"]
spec["register_file_total_bytes"] = spec["regs_per_sm"] * 4 * spec["sms"]
spec["peak_fp32_flops_theory"] = spec["fp32_lanes_total"] * 2 * spec["max_sm_clock_mhz"] * 1e6
spec["peak_dram_bw_theory"] = mem_mhz * 1e6 * 2 * bus_bits / 8
spec["peak_fp16_tensor_flops_spec"] = SPEC_FP16_TENSOR.get(spec["gpu"])
lscpu = {}
try:
out = subprocess.run(["lscpu"], capture_output=True, text=True).stdout
lscpu = {k.strip(): v.strip() for k, v in (l.split(":", 1) for l in out.splitlines() if ":" in l)}
except Exception:
pass
cpu = {
"model": lscpu.get("Model name", platform.processor()),
"logical_cpus": os.cpu_count(),
"torch_threads": torch.get_num_threads(),
"flags_avx512": "avx512f" in lscpu.get("Flags", ""),
"flags_avx2": "avx2" in lscpu.get("Flags", ""),
}
RESULTS["spec"], RESULTS["cpu"] = spec, cpu
rows = [
("GPU", spec["gpu"]), ("compute capability", spec["compute_capability"]),
("streaming multiprocessors (SMs)", spec["sms"]),
("FP32 lanes per SM / total", f"{spec['fp32_lanes_per_sm']} / {spec['fp32_lanes_total']}"),
("warp size", spec["warp_size"]),
("max resident threads per SM (warps)", f"{spec['max_threads_per_sm']} ({spec['max_warps_per_sm']})"),
("register file per SM / whole GPU", f"{spec['regs_per_sm'] * 4 / 1024:.0f} KB / {spec['register_file_total_bytes'] / 2**20:.1f} MB"),
("shared memory per SM", f"{spec['shared_mem_per_sm_bytes'] / 1024:.0f} KB"),
("L2 cache", f"{spec['l2_bytes'] / 2**20:.1f} MB"),
("DRAM", f"{spec['dram_bytes'] / 2**30:.1f} GB, {bus_bits}-bit bus @ {mem_mhz:.0f} MHz"),
("max SM clock / power limit", f"{spec['max_sm_clock_mhz']:.0f} MHz / {spec['power_limit_w']:.0f} W"),
("theoretical FP32 peak", f"{spec['peak_fp32_flops_theory'] / 1e12:.2f} TFLOP/s"),
("datasheet FP16 tensor-core peak", f"{(spec['peak_fp16_tensor_flops_spec'] or float('nan')) / 1e12:.0f} TFLOP/s"),
("theoretical DRAM bandwidth", f"{spec['peak_dram_bw_theory'] / 1e9:.0f} GB/s"),
("host link", spec["pcie"]),
("CPU", f"{cpu['model']} · {cpu['logical_cpus']} logical CPUs · torch uses {cpu['torch_threads']} threads"),
]
print(pd.DataFrame(rows, columns=["", "value"]).to_string(index=False)) value
GPU Tesla T4
compute capability 7.5
streaming multiprocessors (SMs) 40
FP32 lanes per SM / total 64 / 2560
warp size 32
max resident threads per SM (warps) 1024 (32)
register file per SM / whole GPU 256 KB / 10.0 MB
shared memory per SM 64 KB
L2 cache 4.0 MB
DRAM 14.6 GB, 256-bit bus @ 5001 MHz
max SM clock / power limit 1590 MHz / 70 W
theoretical FP32 peak 8.14 TFLOP/s
datasheet FP16 tensor-core peak 65 TFLOP/s
theoretical DRAM bandwidth 320 GB/s
host link gen 3 x16
CPU Intel(R) Xeon(R) CPU @ 2.00GHz · 4 logical CPUs · torch uses 2 threadsPart 1 — the paradox
Two operations, every size, two machines
Matrix multiplication (n × n by n × n, FP32) and element-wise addition (x + y, FP32), each swept from tiny to as large as fits. For each size three timings: the CPU; the GPU with data already in its memory; and the GPU as a naive script would use it — copy the inputs over, compute, copy the answer back.
with experiment("paradox"):
torch.manual_seed(0)
rows = []
for n in [16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536, 2048, 3072, 4096]:
a_c, b_c = torch.randn(n, n), torch.randn(n, n)
a_g, b_g = a_c.to(DEV), b_c.to(DEV)
rows.append(dict(
op="matmul", size=n, flops=2 * n**3, bytes=3 * n * n * 4,
cpu_s=wall_time(lambda: a_c @ b_c),
gpu_s=gpu_time(lambda: a_g @ b_g),
gpu_e2e_s=wall_time(lambda: (a_c.to(DEV) @ b_c.to(DEV)).cpu(), sync=True),
))
for k in range(10, 29):
N = 2**k
x_c, y_c = torch.randn(N), torch.randn(N)
x_g, y_g = x_c.to(DEV), y_c.to(DEV)
o_c, o_g = torch.empty(N), torch.empty(N, device=DEV)
rows.append(dict(
op="add", size=N, flops=N, bytes=3 * N * 4,
cpu_s=wall_time(lambda: torch.add(x_c, y_c, out=o_c)),
gpu_s=gpu_time(lambda: torch.add(x_g, y_g, out=o_g)),
gpu_e2e_s=wall_time(lambda: torch.add(x_c.to(DEV), y_c.to(DEV)).cpu(), sync=True),
))
del x_c, y_c, x_g, y_g, o_c, o_g
torch.cuda.empty_cache()
paradox = pd.DataFrame(rows)
paradox["speedup"] = paradox.cpu_s / paradox.gpu_s
paradox["speedup_e2e"] = paradox.cpu_s / paradox.gpu_e2e_s
RESULTS["paradox"] = paradox.to_dict("records")
fmt = lambda v: f"{v:.3g}"
cols = ["size", "cpu_s", "gpu_s", "gpu_e2e_s", "speedup", "speedup_e2e"]
print("matmul (n × n):")
print(paradox[(paradox.op == "matmul") & paradox["size"].isin([32, 256, 1024, 4096])][cols]
.to_string(index=False, float_format=fmt))
print("\nadd (elements):")
print(paradox[(paradox.op == "add") & paradox["size"].isin([2**10, 2**16, 2**22, 2**28])][cols]
.to_string(index=False, float_format=fmt))
fig, axes = plt.subplots(1, 2, figsize=(9, 3.4), sharey=True)
for ax, op, xlabel in zip(axes, ["matmul", "add"],
["matrix side n (n×n @ n×n, fp32)", "elements (x + y, fp32)"]):
d = paradox[paradox.op == op]
ax.plot(d["size"], d.speedup, "-o", color=BLUE, ms=3, label="GPU, data already on the GPU")
ax.plot(d["size"], d.speedup_e2e, "-o", color=EMBER, ms=3, label="GPU, including the copies")
ax.axhline(1, color=GRAY, lw=0.8, ls="--")
ax.set_xscale("log", base=2); ax.set_yscale("log")
ax.set_xlabel(xlabel); ax.set_title(op)
axes[0].set_ylabel("speed-up over the CPU (×)")
axes[0].legend(frameon=False, fontsize=8)
plt.tight_layout(); plt.show()matmul (n × n):
size cpu_s gpu_s gpu_e2e_s speedup speedup_e2e
32 6.47e-06 3.43e-05 0.000119 0.189 0.0545
256 0.000178 7e-05 0.000402 2.54 0.442
1024 0.0106 0.000852 0.00399 12.4 2.66
4096 0.67 0.0354 0.118 18.9 5.7
add (elements):
size cpu_s gpu_s gpu_e2e_s speedup speedup_e2e
1024 2.7e-06 2.05e-05 8.44e-05 0.132 0.032
65536 1.15e-05 2.33e-05 0.000343 0.494 0.0335
4194304 0.00197 0.000219 0.0115 8.98 0.172
268435456 0.143 0.0131 1.5 11 0.0954
[paradox: 48.3 s]The floor under every tiny operation
At the small end of both sweeps the GPU is not computing — it is waiting to be told what to compute. Each PyTorch call on a CUDA tensor goes through the Python interpreter, the dispatcher and the CUDA driver before the GPU sees a kernel. Here that cost is isolated: a thousand tiny additions issued one by one from Python; the same thousand captured once as a CUDA graph and replayed as a single submission; and the same thousand on the CPU.
with experiment("launch_overhead"):
K = 1000
x_g, x_c = torch.zeros(1024, device=DEV), torch.zeros(1024)
def eager_gpu():
for _ in range(K):
x_g.add_(1.0)
def eager_cpu():
for _ in range(K):
x_c.add_(1.0)
side = torch.cuda.Stream()
side.wait_stream(torch.cuda.current_stream())
with torch.cuda.stream(side):
eager_gpu() # warm-up on a side stream, as capture requires
torch.cuda.current_stream().wait_stream(side)
graph = torch.cuda.CUDAGraph()
with torch.cuda.graph(graph):
eager_gpu()
overhead = {
"gpu_eager_per_op_s": wall_time(eager_gpu, sync=True) / K,
"gpu_graph_per_op_s": wall_time(graph.replay, sync=True) / K,
"cpu_per_op_s": wall_time(eager_cpu) / K,
"gpu_single_roundtrip_s": wall_time(lambda: x_g.add_(1.0), sync=True),
}
RESULTS["launch_overhead"] = overhead
for k, v in overhead.items():
print(f"{k:24s} {v * 1e6:8.2f} µs")gpu_eager_per_op_s 9.30 µs
gpu_graph_per_op_s 1.35 µs
cpu_per_op_s 7.02 µs
gpu_single_roundtrip_s 25.74 µs
[launch_overhead: 0.4 s]Part 2 — threads, warps, and SIMT
The kernels
From here on the GPU is programmed directly, in CUDA C, compiled at run time by NVRTC through CuPy. Every kernel used in the lab is in this one cell. Each is a function run by every thread of a grid; the grid is cut into blocks; each block is scheduled onto one SM; and the SM executes its threads in groups of 32 — warps — that issue one instruction at a time together.
SRC = r'''
// ── SIMT: what happens when the 32 threads of a warp disagree ─────────────────
extern "C" __global__ void diverge(float* out, int paths, int by_lane, int work) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
int group = by_lane ? (threadIdx.x & 31) : (i >> 5); // lane inside a warp, or the warp itself
int p = group % paths;
float x = 1.0f + 1e-7f * (float)i;
if (p == 0) { for (int k = 0; k < work; ++k) x = fmaf(x, 0.9999f, 1e-4f); }
else if (p == 1) { for (int k = 0; k < work; ++k) x = fmaf(x, 1.0001f, -1e-4f); }
else if (p == 2) { for (int k = 0; k < work; ++k) x = fmaf(x, 0.9998f, 2e-4f); }
else { for (int k = 0; k < work; ++k) x = fmaf(x, 1.0002f, -2e-4f); }
out[i] = x;
}
extern "C" __global__ void straggler(float* out, int every, int work) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
float x = 1.0f + 1e-7f * (float)i;
if (i % every == 0) { for (int k = 0; k < work; ++k) x = fmaf(x, 0.9999f, 1e-4f); }
out[i] = x;
}
// ── access pattern: a warp's 32 loads are served together, in 32-byte sectors ─
extern "C" __global__ void strided_read(const float* __restrict__ x, float* __restrict__ y,
long long n, long long mask, int stride) {
long long i = blockIdx.x * (long long)blockDim.x + threadIdx.x;
if (i < n) y[i] = x[(i * stride) & mask];
}
extern "C" __global__ void gather_read(const float* __restrict__ x, const int* __restrict__ idx,
float* __restrict__ y, long long n) {
long long i = blockIdx.x * (long long)blockDim.x + threadIdx.x;
if (i < n) y[i] = x[idx[i]];
}
// ── latency: one thread, and every load needs the result of the previous one ──
extern "C" __global__ void chase(const unsigned int* __restrict__ next, unsigned int start,
int steps, unsigned int* sink, long long* cycles) {
unsigned int p = start;
long long t0 = clock64();
for (int s = 0; s < steps; ++s) p = next[p];
long long t1 = clock64();
sink[0] = p;
cycles[0] = t1 - t0;
}
// ── bandwidth vs working-set size: the whole GPU re-reads one buffer ──────────
// Every block walks the *whole* buffer (from its own starting point), so the footprint of
// each SM is the full working set — the thing whose size the sweep is about.
extern "C" __global__ void reread(const float4* __restrict__ x, long long mask4, int reps,
float* sink) {
long long start = (long long)blockIdx.x * 7919 * blockDim.x; // spread the blocks out
float acc = 0.f;
#pragma unroll 8
for (int r = 0; r < reps; ++r) {
float4 v = x[(start + threadIdx.x + (long long)r * blockDim.x) & mask4];
acc += v.x + v.y + v.z + v.w;
}
if (acc == 1234.5f) sink[0] = acc; // never true; keeps the loads alive
}
// ── latency hiding: a streaming read with ILP independent loads per thread ────
template <int ILP>
__device__ __forceinline__ void stream_read_impl(const float4* __restrict__ x, long long n4,
float* sink) {
long long tid = blockIdx.x * (long long)blockDim.x + threadIdx.x;
long long nth = (long long)gridDim.x * blockDim.x;
float acc = 0.f;
#pragma unroll 1
for (long long base = tid; base < n4; base += nth * ILP) {
float4 v[ILP];
#pragma unroll
for (int j = 0; j < ILP; ++j) {
long long i = base + j * nth;
float4 z = {0.f, 0.f, 0.f, 0.f};
v[j] = (i < n4) ? x[i] : z;
}
#pragma unroll
for (int j = 0; j < ILP; ++j) acc += v[j].x + v[j].y + v[j].z + v[j].w;
}
if (acc == 1234.5f) sink[0] = acc;
}
extern "C" __global__ void stream_read_1(const float4* __restrict__ x, long long n4, float* s) { stream_read_impl<1>(x, n4, s); }
extern "C" __global__ void stream_read_2(const float4* __restrict__ x, long long n4, float* s) { stream_read_impl<2>(x, n4, s); }
extern "C" __global__ void stream_read_4(const float4* __restrict__ x, long long n4, float* s) { stream_read_impl<4>(x, n4, s); }
extern "C" __global__ void stream_read_8(const float4* __restrict__ x, long long n4, float* s) { stream_read_impl<8>(x, n4, s); }
// ── the compute roof: eight independent FMA chains per thread, no memory traffic ─
extern "C" __global__ void fma_peak(float* out, int iters) {
float a0 = 1e-6f * threadIdx.x, a1 = a0 + 0.1f, a2 = a0 + 0.2f, a3 = a0 + 0.3f,
a4 = a0 + 0.4f, a5 = a0 + 0.5f, a6 = a0 + 0.6f, a7 = a0 + 0.7f;
const float b = 0.99999f, c = 1e-5f;
#pragma unroll 4
for (int k = 0; k < iters; ++k) {
a0 = fmaf(a0, b, c); a1 = fmaf(a1, b, c); a2 = fmaf(a2, b, c); a3 = fmaf(a3, b, c);
a4 = fmaf(a4, b, c); a5 = fmaf(a5, b, c); a6 = fmaf(a6, b, c); a7 = fmaf(a7, b, c);
}
out[blockIdx.x * blockDim.x + threadIdx.x] = a0 + a1 + a2 + a3 + a4 + a5 + a6 + a7;
}
// ── the one-knob roofline: read an element, do k FMAs on it, write it back ────
extern "C" __global__ void intensity(const float4* __restrict__ x, float4* __restrict__ y,
long long n4, int k) {
long long tid = blockIdx.x * (long long)blockDim.x + threadIdx.x;
long long nth = (long long)gridDim.x * blockDim.x;
for (long long i = tid; i < n4; i += nth) {
float4 v = x[i];
#pragma unroll 4
for (int j = 0; j < k; ++j) {
v.x = fmaf(v.x, 0.9999f, 1e-4f); v.y = fmaf(v.y, 0.9999f, 1e-4f);
v.z = fmaf(v.z, 0.9999f, 1e-4f); v.w = fmaf(v.w, 0.9999f, 1e-4f);
}
y[i] = v;
}
}
// ── the GPU's own nanosecond clock, independent of the SM clock ───────────────
__device__ __forceinline__ unsigned long long gtimer() {
unsigned long long t;
asm volatile("mov.u64 %0, %%globaltimer;" : "=l"(t));
return t;
}
// ── a clock probe: spin for spin_ns of wall time, count the SM cycles that passed ─
extern "C" __global__ void clock_probe(long long* out, int idx, int spin_ns) {
unsigned long long t0 = gtimer(), t1 = t0;
long long c0 = clock64();
while (t1 - t0 < (unsigned long long)spin_ns) t1 = gtimer();
long long c1 = clock64();
out[3 * idx] = (long long)t0;
out[3 * idx + 1] = (long long)(t1 - t0);
out[3 * idx + 2] = c1 - c0;
}
// ── strided reads with no writes, served from L2 only (ld.global.cg skips L1) ─
extern "C" __global__ void strided_sum_cg(const float* __restrict__ x, long long mask, int stride,
int reps, float* out) {
long long tid = blockIdx.x * (long long)blockDim.x + threadIdx.x;
long long nth = (long long)gridDim.x * blockDim.x;
float acc = 0.f;
#pragma unroll 4
for (int r = 0; r < reps; ++r) acc += __ldcg(&x[((tid + (long long)r * nth) * stride) & mask]);
out[tid] = acc;
}
// ── latency under load: block 0 chases pointers while every other block streams ─
extern "C" __global__ void loaded_chase(const unsigned int* __restrict__ next, unsigned int start,
int steps, const float4* __restrict__ bg, long long bg_n4,
volatile int* done, unsigned long long* bg_loads,
long long* out) {
if (blockIdx.x == 0) {
if (threadIdx.x == 0) {
unsigned int p = start;
unsigned long long t0 = gtimer();
for (int s = 0; s < steps; ++s) p = next[p];
unsigned long long t1 = gtimer();
out[0] = (long long)(t1 - t0);
out[1] = p;
__threadfence();
*done = 1;
}
return;
}
long long nth = (long long)(gridDim.x - 1) * blockDim.x;
long long i = (long long)(blockIdx.x - 1) * blockDim.x + threadIdx.x;
float acc = 0.f;
unsigned long long n = 0;
while (*done == 0) {
#pragma unroll 8
for (int j = 0; j < 32; ++j) {
float4 v = bg[i];
acc += v.x;
i += nth;
if (i >= bg_n4) i -= bg_n4;
}
n += 32;
}
atomicAdd(bg_loads, n);
if (acc == 1234.5f) out[2] = 1;
}
'''
KERNEL_NAMES = ["diverge", "straggler", "strided_read", "gather_read", "chase", "reread",
"stream_read_1", "stream_read_2", "stream_read_4", "stream_read_8",
"fma_peak", "intensity", "clock_probe", "strided_sum_cg", "loaded_chase"]
MOD = cp.RawModule(code=SRC)
KERN = {name: MOD.get_function(name) for name in KERNEL_NAMES}
def launch(name, grid, block, *args):
KERN[name]((int(grid),), (int(block),), args)
def free_gpu():
cp.get_default_memory_pool().free_all_blocks()
torch.cuda.empty_cache()
FULL_BLOCKS = spec["sms"] * (spec["max_threads_per_sm"] // 256) # 256-thread blocks that fill every SM once
def build_chain(nbytes, line=128, order="random", seed=0):
"""A pointer chain built on the GPU: one node per `line` bytes, visited in random or
sequential order. Building it on the GPU keeps the card busy (and at working clocks)
instead of idling for seconds while the CPU shuffles a large permutation."""
words, stride = nbytes // 4, line // 4
if order == "random":
cp.random.seed(seed)
idx = cp.random.permutation(words // stride)
else:
idx = cp.arange(words // stride)
idx = (idx * stride).astype(cp.uint32)
nxt = cp.zeros(words, dtype=cp.uint32)
nxt[idx] = cp.roll(idx, -1)
return nxt, int(idx[0].get())
_WARM = cp.empty(64 * 2**20 // 4, dtype=cp.float32)
_WARM_SINK = cp.zeros(1, dtype=cp.float32)
def warm_up(seconds=0.05):
"""Keep the GPU busy for a moment so a light measurement starts at working clocks."""
t0 = time.perf_counter()
while time.perf_counter() - t0 < seconds:
launch("stream_read_4", FULL_BLOCKS, 256, _WARM, np.int64(_WARM.size // 4), _WARM_SINK)
cp.cuda.Device().synchronize()
print(f"compiled {len(KERN)} kernels; one full wave = {FULL_BLOCKS} blocks × 256 threads")compiled 15 kernels; one full wave = 160 blocks × 256 threadsDivergence: one instruction stream per warp
The diverge kernel gives each thread one of paths different loops to run. The only thing that changes between rows is who gets which loop: either whole warps agree (warp 0 takes path 0, warp 1 takes path 1, …) or the split runs through every warp (lane 0 takes path 0, lane 1 path 1, …). The total arithmetic is identical in every row.
Each configuration is measured seven times, in a shuffled order, so that a drifting clock can’t masquerade as an effect: the spread between repeats is the noise floor, measured rather than guessed.
with experiment("simt_divergence"):
n = FULL_BLOCKS * 256 * 4
out = cp.empty(n, dtype=cp.float32)
WORK = 20000
configs = [(by_lane, paths) for by_lane in (0, 1) for paths in (1, 2, 4)]
trials = {c: [] for c in configs}
rng = np.random.default_rng(0)
for rep in range(7):
for j in rng.permutation(len(configs)):
by_lane, paths = configs[j]
trials[(by_lane, paths)].append(gpu_time(
lambda: launch("diverge", n // 256, 256, out, np.int32(paths), np.int32(by_lane),
np.int32(WORK)), iters=5, warmup=1))
rows = []
for (by_lane, paths), ts in trials.items():
ts = np.array(ts) * 1e3
rows.append(dict(split="lanes inside every warp" if by_lane else "whole warps", paths=paths,
ms=float(np.median(ts)), ms_min=float(ts.min()), ms_max=float(ts.max())))
div = pd.DataFrame(rows)
base = div.ms.iloc[0]
for c in ("ms", "ms_min", "ms_max"):
div["vs_uniform" + c[2:]] = div[c] / base
RESULTS["divergence"] = div.to_dict("records")
print(div.to_string(index=False, float_format=lambda v: f"{v:.3f}"))
del out; free_gpu() split paths ms ms_min ms_max vs_uniform vs_uniform_min vs_uniform_max
whole warps 1 0.895 0.877 0.921 1.000 0.979 1.029
whole warps 2 0.885 0.872 0.903 0.989 0.974 1.009
whole warps 4 0.885 0.836 0.901 0.989 0.934 1.006
lanes inside every warp 1 0.893 0.848 0.931 0.998 0.948 1.040
lanes inside every warp 2 1.770 1.688 1.790 1.978 1.886 2.000
lanes inside every warp 4 3.511 3.301 3.594 3.923 3.688 4.016
[simt_divergence: 0.5 s]The same fact from the other side: one thread in every does a long loop and the rest do nothing. If the warp is the unit of execution, the time depends not on how many threads work but on how many warps contain at least one that does.
with experiment("simt_straggler"):
n, WORK = FULL_BLOCKS * 256 * 4, 20000
out = cp.empty(n, dtype=cp.float32)
rows = []
for every in [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048]:
t = gpu_time(lambda: launch("straggler", n // 256, 256, out, np.int32(every), np.int32(WORK)))
rows.append(dict(every=every, busy_threads=1 / every,
busy_warps=min(1.0, 32 / every), ms=t * 1e3))
strag = pd.DataFrame(rows)
strag["vs_all_busy"] = strag.ms / strag.ms.iloc[0]
# the floor: one warp alone, nothing to interleave with — the latency of its dependent chain
one_warp_ms = gpu_time(lambda: launch("straggler", 1, 32, out, np.int32(1), np.int32(WORK))) * 1e3
RESULTS["straggler"] = strag.to_dict("records")
RESULTS["straggler_one_warp_ms"] = one_warp_ms
print(strag.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
print(f"\none warp on its own: {one_warp_ms:.4f} ms "
f"(= {one_warp_ms / 1e3 / WORK * 1e9:.2f} ns per dependent FMA)")
fig, ax = plt.subplots(figsize=(6.4, 3.2))
ax.plot(strag.every, strag.vs_all_busy, "-o", color=BLUE, ms=4, label="measured time")
ax.plot(strag.every, strag.busy_threads, ":", color=GRAY, label="fraction of threads that work")
ax.plot(strag.every, strag.busy_warps, "--", color=EMBER, label="fraction of warps with ≥1 working thread")
ax.set_xscale("log", base=2); ax.set_yscale("log")
ax.set_xlabel("one working thread in every …"); ax.set_ylabel("relative to all threads working")
ax.legend(frameon=False, fontsize=8); plt.tight_layout(); plt.show()
del out; free_gpu() every busy_threads busy_warps ms vs_all_busy
1 1 1 0.8466 1
2 0.5 1 0.8451 0.9982
4 0.25 1 0.8519 1.006
8 0.125 1 0.8517 1.006
16 0.0625 1 0.8504 1.004
32 0.03125 1 0.8438 0.9967
64 0.01562 0.5 0.4825 0.57
128 0.007812 0.25 0.2956 0.3491
256 0.003906 0.125 0.241 0.2847
512 0.001953 0.0625 0.1311 0.1548
1024 0.0009766 0.03125 0.1167 0.1379
2048 0.0004883 0.01562 0.07357 0.0869
one warp on its own: 0.0696 ms (= 3.48 ns per dependent FMA)
[simt_straggler: 0.5 s]Part 3 — memory
Coalescing: a warp reads together
Each thread copies one float, but reads it from x[i * stride]. With stride 1 the 32 threads of a warp touch 128 consecutive bytes; with stride 8 or more every thread lands in its own 32-byte sector, and the memory system moves 32 bytes to deliver 4. The last row is a random gather. The first “predicted” column assumes nothing but that sector size. The second was added after the first run of this lab, when the measurements kept falling past stride 8: it assumes that DRAM is read 64 bytes at a time (two sectors), whatever the caches do.
with experiment("coalescing"):
n = 2**24
x = cp.random.random(2**28, dtype=cp.float32) # 1 GB source, far beyond L2
y = cp.empty(n, dtype=cp.float32)
rows = []
for stride in [1, 2, 4, 8, 16, 32, 64]:
t = gpu_time(lambda: launch("strided_read", n // 256, 256, x, y, np.int64(n),
np.int64(2**28 - 1), np.int32(stride)))
# bytes moved per element: the 4-byte write + whatever the read drags in
rows.append(dict(pattern=f"stride {stride}", stride=stride, useful_gb_s=8 * n / t / 1e9,
predicted_32b=8 / (4 + 4 * min(stride, 8)),
predicted_64b=8 / (4 + 4 * min(stride, 16))))
idx = cp.random.randint(0, 2**28, n, dtype=cp.int32)
t = gpu_time(lambda: launch("gather_read", n // 256, 256, x, idx, y, np.int64(n)))
rows.append(dict(pattern="random gather", stride=None, useful_gb_s=8 * n / t / 1e9,
predicted_32b=8 / (4 + 4 + 32), predicted_64b=8 / (4 + 4 + 64)))
coal = pd.DataFrame(rows)
coal["measured_fraction"] = coal.useful_gb_s / coal.useful_gb_s.iloc[0]
RESULTS["coalescing"] = coal.to_dict("records")
print(coal.to_string(index=False, float_format=lambda v: f"{v:.3f}"))
del x, y, idx; free_gpu() pattern stride useful_gb_s predicted_32b predicted_64b measured_fraction
stride 1 1.000 238.625 1.000 1.000 1.000
stride 2 2.000 167.611 0.667 0.667 0.702
stride 4 4.000 106.059 0.400 0.400 0.444
stride 8 8.000 59.086 0.222 0.222 0.248
stride 16 16.000 29.455 0.222 0.118 0.123
stride 32 32.000 27.211 0.222 0.118 0.114
stride 64 64.000 25.218 0.222 0.118 0.106
random gather NaN 14.181 0.200 0.111 0.059
[coalescing: 1.2 s]Where does a 64-byte granularity come from — the caches, or DRAM? The same strided read again, this time with no writes at all and with loads that skip L1 (__ldcg), from two sources: a 2 MB buffer that lives in L2, and a 1 GB buffer that has to come from DRAM. If the extra halving between stride 8 and stride 16 belongs to DRAM, it appears only in the second.
with experiment("coalescing_l2_vs_dram"):
nth = FULL_BLOCKS * 256
reps = 3200
out = cp.empty(nth, dtype=cp.float32)
rows = []
for label, nbytes in [("L2 (2 MB source)", 2 * 2**20), ("DRAM (1 GB source)", 2**30)]:
src = cp.random.random(nbytes // 4, dtype=cp.float32)
mask = np.int64(nbytes // 4 - 1)
for stride in [1, 2, 4, 8, 16, 32, 64]:
t = gpu_time(lambda: launch("strided_sum_cg", FULL_BLOCKS, 256, src, mask, np.int32(stride),
np.int32(reps), out), iters=7)
rows.append(dict(source=label, stride=stride, useful_gb_s=nth * reps * 4 / t / 1e9))
del src; free_gpu()
cl2 = pd.DataFrame(rows)
cl2["fraction_of_stride_1"] = cl2.useful_gb_s / cl2.groupby("source").useful_gb_s.transform("first")
cl2["predicted_32b"] = 1 / np.minimum(cl2.stride, 8)
cl2["predicted_64b"] = 1 / np.minimum(cl2.stride, 16)
RESULTS["coalescing_l2_vs_dram"] = cl2.to_dict("records")
print(cl2.to_string(index=False, float_format=lambda v: f"{v:.3f}"))
del out; free_gpu() source stride useful_gb_s fraction_of_stride_1 predicted_32b predicted_64b
L2 (2 MB source) 1 1325.137 1.000 1.000 1.000
L2 (2 MB source) 2 677.585 0.511 0.500 0.500
L2 (2 MB source) 4 330.323 0.249 0.250 0.250
L2 (2 MB source) 8 146.650 0.111 0.125 0.125
L2 (2 MB source) 16 159.970 0.121 0.125 0.062
L2 (2 MB source) 32 109.809 0.083 0.125 0.062
L2 (2 MB source) 64 110.979 0.084 0.125 0.062
DRAM (1 GB source) 1 239.403 1.000 1.000 1.000
DRAM (1 GB source) 2 125.857 0.526 0.500 0.500
DRAM (1 GB source) 4 67.160 0.281 0.250 0.250
DRAM (1 GB source) 8 34.206 0.143 0.125 0.125
DRAM (1 GB source) 16 17.088 0.071 0.125 0.062
DRAM (1 GB source) 32 15.475 0.065 0.125 0.062
DRAM (1 GB source) 64 13.389 0.056 0.125 0.062
[coalescing_l2_vs_dram: 1.7 s]From L2 the fall stops at stride 8: every thread has its own 32-byte sector and nothing more is wasted, which is the sector model exactly. From DRAM it halves once more at stride 16. The caches work in 32-byte sectors; DRAM is read 64 bytes at a time.
Latency: how long one load takes, level by level
A single thread follows a chain of pointers laid out in random order, one per 128-byte line, so every load depends on the one before it and nothing can be overlapped or prefetched. Grow the chain and it spills out of L1, then out of L2, into DRAM. The time per step is the latency of whichever level the chain fits in. clock64() counts SM cycles inside the kernel, so the ratio of cycles to nanoseconds is the clock the SM actually ran at.
The chains are built on the GPU itself, and the card is kept busy for a moment right before each measurement. (A first version built them on the CPU; for the largest chain that left the GPU idle for seconds, it dropped to idle clocks, and the measurement came out 70% slow for a reason that had nothing to do with memory.) The sweep runs to 8 GB — half the card — to see whether anything lies beyond DRAM’s own latency, such as the cost of address translation.
with experiment("latency"):
sink, cyc = cp.zeros(1, dtype=cp.uint32), cp.zeros(1, dtype=cp.int64)
STEPS = 50_000
rows = []
for k in range(12, 34):
try:
chain, start = build_chain(2**k)
except cp.cuda.memory.OutOfMemoryError:
print(f"{2**k / 2**30:.0f} GB chain does not fit — stopping the sweep here")
break
warm_up()
t = gpu_time(lambda: launch("chase", 1, 1, chain, np.uint32(start), np.int32(STEPS), sink, cyc),
iters=3, warmup=1)
cycles = int(cyc.get()[0]) / STEPS
rows.append(dict(bytes=2**k, ns=t / STEPS * 1e9, cycles=cycles,
sm_ghz=cycles / (t / STEPS * 1e9)))
del chain; free_gpu()
lat = pd.DataFrame(rows)
RESULTS["latency"] = lat.to_dict("records")
print(lat.to_string(index=False, float_format=lambda v: f"{v:.3g}"))
fig, ax = plt.subplots(figsize=(6.4, 3.2))
ax.plot(lat["bytes"], lat.ns, "-o", color=BLUE, ms=4)
ax.axvline(spec["l2_bytes"], color=EMBER, ls="--", lw=0.8)
ax.text(spec["l2_bytes"], ax.get_ylim()[1] * 0.95, " L2 size", color=EMBER, fontsize=8, va="top")
ax.set_xscale("log", base=2)
ax.set_xlabel("working set of the pointer chain (bytes)"); ax.set_ylabel("ns per dependent load")
plt.tight_layout(); plt.show() bytes ns cycles sm_ghz
4096 26.7 40.1 1.5
8192 27 40.2 1.49
16384 27.5 40.5 1.47
32768 27.8 41 1.48
65536 96.4 147 1.52
131072 150 232 1.55
262144 152 232 1.53
524288 150 232 1.55
1048576 151 232 1.54
2097152 150 232 1.55
4194304 162 249 1.53
8388608 314 491 1.56
16777216 314 489 1.56
33554432 311 486 1.56
67108864 316 492 1.56
134217728 316 493 1.56
268435456 316 493 1.56
536870912 317 494 1.56
1073741824 317 493 1.56
2147483648 317 494 1.56
4294967296 317 495 1.56
8589934592 326 491 1.51
[latency: 4.8 s]Latency hiding: Little’s law on a GPU
The chain above is the worst case: one thread, nothing else to do while it waits. A real kernel has thousands of loads in flight at once, and while one warp waits the scheduler issues another. How many must be in flight? Little’s law: bytes in flight = latency × bandwidth. Here a streaming read is launched with a controlled number of warps, each thread keeping 1, 2, 4 or 8 independent 16-byte loads outstanding (ILP). The warps are spread one per block, so they land on as many different SMs as possible. If Little’s law is right, the four curves are four different functions of warps but one function of bytes in flight.
A second series packs the same warps into a single block, so they all run on one SM. That separates the GPU-wide question (how much must be in flight?) from a per-SM one (how much can one SM keep in flight?).
with experiment("latency_hiding"):
n4 = 256 * 2**20 // 16 # 256 MB, far beyond L2
xs = cp.random.random(n4 * 4, dtype=cp.float32)
sink = cp.zeros(1, dtype=cp.float32)
resident_warps = spec["sms"] * spec["max_warps_per_sm"]
max_blocks_per_sm = prop("maxBlocksPerMultiProcessor") or (16 if cc < (8, 0) else 32)
lat_df = pd.DataFrame(RESULTS.get("latency", []))
past_l2 = lat_df[(lat_df["bytes"] >= 4 * spec["l2_bytes"]) & (lat_df["bytes"] <= 16 * spec["l2_bytes"])] if len(lat_df) else lat_df
dram_latency_s = float(past_l2.ns.median()) * 1e-9 if len(past_l2) else 500e-9
rows = []
for ilp in (1, 2, 4, 8):
for warps in [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]:
# one warp per block while blocks can still all be resident; fatter blocks beyond that
per_block = max(1, -(-warps // (spec["sms"] * max_blocks_per_sm)))
per_block = 1 << (per_block - 1).bit_length() # round up to a power of two
block, grid = 32 * per_block, warps // per_block
t = gpu_time(lambda: launch(f"stream_read_{ilp}", grid, block, xs, np.int64(n4), sink),
iters=5, warmup=1)
live = min(warps, resident_warps)
rows.append(dict(layout="spread", ilp=ilp, warps=warps, sms_used=min(grid, spec["sms"]),
resident_warps=live, bytes_in_flight=live * 32 * 16 * ilp,
gb_s=n4 * 16 / t / 1e9))
for warps in [1, 2, 4, 8, 16, 32]: # the same warps, all on one SM
n4_small = 32 * 2**20 // 16 # 32 MB: one SM is slow, keep it short
t = gpu_time(lambda: launch(f"stream_read_{ilp}", 1, 32 * warps, xs, np.int64(n4_small), sink),
iters=5, warmup=1)
rows.append(dict(layout="one SM", ilp=ilp, warps=warps, sms_used=1, resident_warps=warps,
bytes_in_flight=warps * 32 * 16 * ilp, gb_s=n4_small * 16 / t / 1e9))
hide = pd.DataFrame(rows)
spread = hide[hide.layout == "spread"]
one_sm = hide[hide.layout == "one SM"]
bw_max = spread.gb_s.max() * 1e9
hide["littles_law_gb_s"] = np.minimum(bw_max, hide.bytes_in_flight / dram_latency_s) / 1e9
RESULTS["latency_hiding"] = {"rows": hide.to_dict("records"), "dram_latency_s": dram_latency_s,
"bw_max": bw_max, "one_sm_max": one_sm.gb_s.max() * 1e9}
print(f"DRAM latency used for Little's law: {dram_latency_s * 1e9:.0f} ns; "
f"best streaming read {bw_max / 1e9:.0f} GB/s")
print(f"=> bytes that must be in flight to saturate: {bw_max * dram_latency_s / 1024:.0f} KB")
print(f"one SM on its own tops out at {one_sm.gb_s.max():.1f} GB/s "
f"({spec['sms']} of them: {spec['sms'] * one_sm.gb_s.max():.0f} GB/s)")
print("\nspread over the SMs (GB/s):")
print(spread.pivot(index="warps", columns="ilp", values="gb_s").round(1).to_string())
print("\nall on one SM (GB/s):")
print(one_sm.pivot(index="warps", columns="ilp", values="gb_s").round(1).to_string())
fig, axes = plt.subplots(1, 2, figsize=(9, 3.4), sharey=True)
colors = {1: GRAY, 2: BLUE, 4: FOREST, 8: EMBER}
for ilp, d in spread.groupby("ilp"):
axes[0].plot(d.warps, d.gb_s, "-o", ms=3, color=colors[ilp], label=f"ILP {ilp}")
axes[1].plot(d.bytes_in_flight, d.gb_s, "o", ms=3, color=colors[ilp])
ll = hide[hide.layout == "spread"].sort_values("bytes_in_flight")
axes[1].plot(ll.bytes_in_flight, ll.littles_law_gb_s, "--", color=INK, lw=0.9,
label="latency × bandwidth")
axes[0].set_xscale("log", base=2); axes[1].set_xscale("log", base=2); axes[0].set_yscale("log")
axes[0].set_xlabel("warps launched"); axes[1].set_xlabel("bytes in flight")
axes[0].set_ylabel("read bandwidth (GB/s)")
axes[0].legend(frameon=False, fontsize=8); axes[1].legend(frameon=False, fontsize=8)
plt.tight_layout(); plt.show()
del xs; free_gpu()DRAM latency used for Little's law: 314 ns; best streaming read 275 GB/s
=> bytes that must be in flight to saturate: 84 KB
one SM on its own tops out at 48.9 GB/s (40 of them: 1957 GB/s)
spread over the SMs (GB/s):
ilp 1 2 4 8
warps
1 1.5 2.7 5.1 8.6
2 3.0 5.4 10.2 17.0
4 6.0 10.9 20.2 34.0
8 11.9 21.6 39.9 66.2
16 23.8 43.0 78.6 130.2
32 47.3 84.3 149.3 236.7
64 92.8 160.5 251.1 264.1
128 174.3 253.5 255.0 266.1
256 257.0 264.3 258.0 265.5
512 264.0 265.9 253.6 263.7
1024 271.9 265.6 256.0 264.0
2048 275.2 265.5 256.2 264.3
4096 274.7 266.2 259.0 263.9
all on one SM (GB/s):
ilp 1 2 4 8
warps
1 1.5 2.6 4.9 7.2
2 3.0 5.3 9.6 15.5
4 5.9 10.4 18.6 31.3
8 11.7 19.7 35.3 46.4
16 22.7 35.3 48.9 47.8
32 39.6 48.5 48.9 48.0
[latency_hiding: 6.5 s]Latency under load
Little’s law above used the latency of an idle memory system — one pointer chase with nothing else running. Under load, requests queue, and each one should take longer. That is measured directly here: block 0 of the grid chases pointers through a 256 MB chain while every other block streams through a separate 1 GB buffer, and the number of streaming blocks is swept from none to the whole GPU. For each setting: the latency the chase saw, and the bandwidth the background got.
with experiment("loaded_latency"):
chain, start = build_chain(256 * 2**20)
bg = cp.random.random(2**30 // 4, dtype=cp.float32)
bg_n4 = np.int64(bg.size // 4)
STEPS = 20_000
rows = []
for blocks in [0, 1, 2, 4, 8, 16, 24, 32, 48, 64, 96, 128, FULL_BLOCKS - 1]:
samples = []
for _ in range(3):
done = cp.zeros(1, dtype=cp.int32)
loads = cp.zeros(1, dtype=cp.uint64)
out = cp.zeros(3, dtype=cp.int64)
warm_up()
launch("loaded_chase", 1 + blocks, 256, chain, np.uint32(start), np.int32(STEPS), bg, bg_n4,
done, loads, out)
cp.cuda.Device().synchronize()
chase_ns = int(out.get()[0])
samples.append((chase_ns / STEPS, int(loads.get()[0]) * 16 / chase_ns)) # ns per load, GB/s
lat_ns, bw = np.median([s[0] for s in samples]), np.median([s[1] for s in samples])
rows.append(dict(background_blocks=blocks, latency_ns=float(lat_ns), background_gb_s=float(bw)))
loaded = pd.DataFrame(rows)
RESULTS["loaded_latency"] = loaded.to_dict("records")
del chain, bg; free_gpu()
print(loaded.to_string(index=False, float_format=lambda v: f"{v:.1f}"))
fig, ax = plt.subplots(figsize=(6.4, 3.2))
ax.plot(loaded.background_gb_s, loaded.latency_ns, "-o", color=EMBER, ms=4)
ax.set_xlabel("bandwidth used by everything else (GB/s)"); ax.set_ylabel("latency of one load (ns)")
plt.tight_layout(); plt.show() background_blocks latency_ns background_gb_s
0 328.6 0.0
1 331.6 38.9
2 338.6 76.7
4 357.3 147.4
8 455.2 260.3
16 1165.0 270.3
24 1693.4 270.1
32 2148.1 268.4
48 3631.9 268.1
64 4582.0 267.7
96 7029.4 267.4
128 9541.0 265.4
159 9595.9 264.1
[loaded_latency: 4.6 s]Up to about 150 GB/s of background traffic the latency barely moves; near saturation it climbs to ~450 ns — the latency the Little’s-law loads actually saw at their knee. Past saturation the background gets no more bandwidth at all, and the latency explodes into microseconds: more requests in flight than Little’s law asks for only lengthen the queue.
Bandwidth: what each level delivers
Latency is the cost of one load; bandwidth is how many bytes per second the level delivers when the whole GPU is asking. A full-occupancy grid re-reads one buffer of a given size, over and over, and every block walks the whole of it. While the buffer fits in an SM’s L1, the reads never leave the SM; up to the L2 size they come from L2; beyond it, from DRAM.
with experiment("hierarchy_bandwidth"):
nth = FULL_BLOCKS * 256
sink = cp.zeros(1, dtype=cp.float32)
reps = int(max(64, 2e9 / (nth * 16)) // 8 * 8) # ~2 GB of loads per launch at every size
rows = []
for k in range(12, 31):
buf = cp.random.random(2**k // 4, dtype=cp.float32)
t = gpu_time(lambda: launch("reread", FULL_BLOCKS, 256, buf, np.int64(2**k // 16 - 1),
np.int32(reps), sink), iters=5)
rows.append(dict(bytes=2**k, gb_s=nth * reps * 16 / t / 1e9))
del buf; free_gpu()
hbw = pd.DataFrame(rows)
RESULTS["hierarchy_bandwidth"] = hbw.to_dict("records")
print(hbw.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
fig, ax = plt.subplots(figsize=(6.4, 3.2))
ax.plot(hbw["bytes"], hbw.gb_s, "-o", color=FOREST, ms=4)
ax.axvline(spec["l2_bytes"], color=EMBER, ls="--", lw=0.8)
ax.axhline(spec["peak_dram_bw_theory"] / 1e9, color=GRAY, ls=":", lw=0.8)
ax.text(hbw["bytes"].iloc[0], spec["peak_dram_bw_theory"] / 1e9, "DRAM datasheet", color=GRAY,
fontsize=8, va="bottom")
ax.set_xscale("log", base=2); ax.set_yscale("log")
ax.set_xlabel("working set (bytes)"); ax.set_ylabel("read bandwidth (GB/s)")
plt.tight_layout(); plt.show() bytes gb_s
4096 3870
8192 3882
16384 3870
32768 3870
65536 3820
131072 1561
262144 1336
524288 1114
1048576 1107
2097152 1094
4194304 1088
8388608 642.5
16777216 430.6
33554432 323.8
67108864 293.8
134217728 271
268435456 268.8
536870912 267.8
1073741824 267.1
[hierarchy_bandwidth: 0.9 s]How big is L1, exactly?
On Turing, L1 and shared memory are one 96 KB array per SM, and the driver decides per kernel how much of it becomes L1. A kernel can state a preference, the carveout: 0 asks for as much L1 as possible, 100 for as much shared memory as possible. The chase is repeated here in fine steps under each setting, with two chain orders: random (as above) and sequential, which is immune to the unlucky set conflicts a random chain can hit just below capacity. The bandwidth kernel is repeated at the power-of-two sizes it supports.
with experiment("l1_capacity"):
sink, cyc = cp.zeros(1, dtype=cp.uint32), cp.zeros(1, dtype=cp.int64)
bsink = cp.zeros(1, dtype=cp.float32)
nth, reps = FULL_BLOCKS * 256, 3048
rows, notes = [], {}
for label, carveout in [("default", None), ("max L1 (carveout 0)", 0), ("max shared (carveout 100)", 100)]:
if carveout is not None:
try:
for name in ("chase", "reread"):
KERN[name].preferred_shared_memory_carveout = carveout
except Exception as e:
notes[label] = f"could not set carveout: {e}"
continue
for kb in [8, 16, 24, 32, 40, 48, 56, 64, 72, 80, 96, 112, 128, 160, 192]:
for order in ("random", "sequential"):
chain, start = build_chain(kb * 1024, order=order)
t = gpu_time(lambda: launch("chase", 1, 1, chain, np.uint32(start), np.int32(20000), sink, cyc),
iters=3, warmup=1)
rows.append(dict(carveout=label, kb=kb, kind=f"chase, {order}", ns=t / 20000 * 1e9,
cycles=int(cyc.get()[0]) / 20000))
del chain
if kb & (kb - 1) == 0: # the bandwidth kernel needs a power of two
buf = cp.random.random(kb * 1024 // 4, dtype=cp.float32)
t = gpu_time(lambda: launch("reread", FULL_BLOCKS, 256, buf, np.int64(kb * 1024 // 16 - 1),
np.int32(reps), bsink), iters=5)
rows.append(dict(carveout=label, kb=kb, kind="bandwidth", gb_s=nth * reps * 16 / t / 1e9))
del buf
free_gpu()
try:
for name in ("chase", "reread"):
KERN[name].preferred_shared_memory_carveout = -1 # back to the driver's default
except Exception as e:
notes["reset"] = str(e)
l1 = pd.DataFrame(rows)
seq = l1[l1.kind == "chase, sequential"]
caps = {c: int(d[d.cycles <= 1.3 * d.cycles.min()].kb.max()) for c, d in seq.groupby("carveout")}
RESULTS["l1_capacity"] = {"rows": l1.to_dict("records"), "capacity_kb": caps, "notes": notes}
print("L1 capacity seen by a sequential chase (KB):", caps, notes or "")
print(l1.pivot_table(index="kb", columns=["carveout", "kind"], values="cycles").round(0).to_string())
bw = l1[l1.kind == "bandwidth"]
if len(bw):
print("\nbandwidth (GB/s):")
print(bw.pivot_table(index="kb", columns="carveout", values="gb_s").round(0).to_string())L1 capacity seen by a sequential chase (KB): {'default': 56, 'max L1 (carveout 0)': 56, 'max shared (carveout 100)': 24}
carveout default max L1 (carveout 0) max shared (carveout 100)
kind chase, random chase, sequential chase, random chase, sequential chase, random chase, sequential
kb
8 41.0 41.0 41.0 41.0 41.0 41.0
16 41.0 41.0 41.0 41.0 41.0 41.0
24 42.0 42.0 42.0 42.0 42.0 42.0
32 42.0 42.0 42.0 42.0 138.0 138.0
40 43.0 43.0 43.0 43.0 195.0 195.0
48 44.0 44.0 44.0 44.0 227.0 227.0
56 44.0 44.0 44.0 44.0 227.0 232.0
64 115.0 118.0 115.0 128.0 232.0 232.0
72 191.0 198.0 191.0 203.0 233.0 232.0
80 231.0 232.0 231.0 232.0 232.0 232.0
96 232.0 232.0 232.0 232.0 232.0 232.0
112 233.0 232.0 232.0 232.0 232.0 232.0
128 232.0 233.0 232.0 232.0 232.0 232.0
160 232.0 232.0 231.0 232.0 232.0 232.0
192 232.0 232.0 232.0 232.0 232.0 232.0
bandwidth (GB/s):
carveout default max L1 (carveout 0) max shared (carveout 100)
kb
8 3903.0 3928.0 3879.0
16 3907.0 3917.0 3870.0
32 3901.0 3916.0 2713.0
64 3840.0 3841.0 1623.0
128 1775.0 1767.0 1432.0
[l1_capacity: 1.2 s]A kernel with no shared memory already gets the large L1: “default” and “max L1” coincide, with 56 KB running at full L1 speed out of the 64 KB configuration. Asking for maximum shared memory shrinks L1 to the 32 KB configuration, 24 KB usable, and the bandwidth at 32 KB drops with it. (The coarse chase above, stepping from 32 KB straight to 64 KB, could only say “between the two”.)
The level below DRAM: the host link
The CPU’s memory is one more level of the hierarchy, the slowest and furthest away, reached over PCIe. Ordinary (“pageable”) host memory has to be staged through a pinned buffer by the driver; memory allocated as pinned can be read by the GPU’s copy engine directly.
with experiment("pcie"):
rows = []
for nbytes in [4, 4096, 2**20, 16 * 2**20, 256 * 2**20, 1024 * 2**20]:
h = torch.empty(nbytes, dtype=torch.uint8)
hp = h.pin_memory()
d = torch.empty(nbytes, dtype=torch.uint8, device=DEV)
rows.append(dict(
bytes=nbytes,
h2d_pageable_s=wall_time(lambda: d.copy_(h), sync=True),
h2d_pinned_s=wall_time(lambda: d.copy_(hp, non_blocking=True), sync=True),
d2h_pinned_s=wall_time(lambda: hp.copy_(d, non_blocking=True), sync=True),
))
del h, hp, d
pcie = pd.DataFrame(rows)
for c in ["h2d_pageable", "h2d_pinned", "d2h_pinned"]:
pcie[c + "_gb_s"] = pcie["bytes"] / pcie[c + "_s"] / 1e9
RESULTS["pcie"] = pcie.to_dict("records")
free_gpu()
print(pcie.to_string(index=False, float_format=lambda v: f"{v:.3g}")) bytes h2d_pageable_s h2d_pinned_s d2h_pinned_s h2d_pageable_gb_s h2d_pinned_gb_s d2h_pinned_gb_s
4 2.13e-05 2.12e-05 1.97e-05 0.000188 0.000188 0.000203
4096 2.44e-05 2.11e-05 1.98e-05 0.168 0.194 0.207
1048576 0.00035 0.000107 9.64e-05 2.99 9.79 10.9
16777216 0.00366 0.00138 0.0013 4.59 12.2 12.9
268435456 0.0401 0.0218 0.0206 6.7 12.3 13
1073741824 0.161 0.087 0.0823 6.69 12.3 13
[pcie: 7.3 s]Part 4 — FLOPs
The CUDA-core roof
fma_peak never touches memory: every thread runs eight independent chains of fused multiply-adds in registers. That is the most arithmetic the FP32 lanes can do. It runs for a couple of seconds under load while nvidia-smi is polled — the T4 is a 70-watt card, and the question is whether it holds its boost clock.
with experiment("compute_roof"):
out = cp.empty(FULL_BLOCKS * 8 * 256, dtype=cp.float32)
nth = out.size
iters = int(0.02 * spec["peak_fp32_flops_theory"] / (nth * 16)) # ~20 ms per launch
t, clocks = sustained(lambda: launch("fma_peak", FULL_BLOCKS * 8, 256, out, np.int32(iters)))
fma = {"flops_per_s": nth * iters * 16 / t, "clocks": clocks}
fma["theory_at_measured_clock"] = spec["fp32_lanes_total"] * 2 * clocks.get("clocks.sm", np.nan) * 1e6
gemm = {}
for label, n, dtype in [("fp32 4096", 4096, torch.float32), ("fp16 4096", 4096, torch.float16),
("fp16 8192", 8192, torch.float16)]:
A = torch.randn(n, n, device=DEV, dtype=dtype); B = torch.randn(n, n, device=DEV, dtype=dtype)
t, clocks = sustained(lambda: A @ B)
gemm[label] = {"flops_per_s": 2 * n**3 / t, "clocks": clocks}
del A, B
RESULTS["compute_roof"] = {"fma_peak": fma, "gemm": gemm}
del out; free_gpu()
rows = [("FMA kernel, CUDA cores, fp32", fma)] + [(f"cuBLAS matmul {k}", v) for k, v in gemm.items()]
print(f"{'':32s} {'TFLOP/s':>8s} {'SM MHz':>7s} {'W':>6s} {'°C':>4s}")
for name, r in rows:
c = r["clocks"]
print(f"{name:32s} {r['flops_per_s'] / 1e12:8.2f} {c.get('clocks.sm', 0):7.0f} "
f"{c.get('power.draw', 0):6.1f} {c.get('temperature.gpu', 0):4.0f}")
print(f"\ntheoretical fp32 at max clock {spec['peak_fp32_flops_theory'] / 1e12:.2f} TFLOP/s; "
f"at the clock measured under load {fma['theory_at_measured_clock'] / 1e12:.2f} TFLOP/s")
if spec["peak_fp16_tensor_flops_spec"]:
print(f"datasheet fp16 tensor-core peak {spec['peak_fp16_tensor_flops_spec'] / 1e12:.0f} TFLOP/s") TFLOP/s SM MHz W °C
FMA kernel, CUDA cores, fp32 6.81 1365 67.1 68
cuBLAS matmul fp32 4096 3.52 705 66.2 69
cuBLAS matmul fp16 4096 20.26 555 66.8 70
cuBLAS matmul fp16 8192 20.81 570 66.5 71
theoretical fp32 at max clock 8.14 TFLOP/s; at the clock measured under load 6.99 TFLOP/s
datasheet fp16 tensor-core peak 65 TFLOP/s
[compute_roof: 8.5 s]Does the roof sink? A 70-watt card under load
The sustained numbers above were taken after a warm-up. The obvious expectation is that the card boosts at first and sinks as the power limit catches up, so this checks it: the card rests for ten seconds, then multiplies the same pair of 4096×4096 FP16 matrices back to back for five seconds. Every multiplication is timed, and the SM clock and power draw are polled alongside.
with experiment("power_cap"):
n = 4096
A = torch.randn(n, n, device=DEV, dtype=torch.float16)
B = torch.randn(n, n, device=DEV, dtype=torch.float16)
A @ B; torch.cuda.synchronize()
time.sleep(10) # let the card cool and clock back up
polled, stop = [], threading.Event()
t_start = time.perf_counter()
def poll():
while not stop.is_set():
s = smi("clocks.sm,power.draw,temperature.gpu")
s["t"] = time.perf_counter() - t_start
polled.append(s)
stop.wait(0.05)
poller = threading.Thread(target=poll, daemon=True); poller.start()
ev_s, ev_e = torch.cuda.Event(enable_timing=True), torch.cuda.Event(enable_timing=True)
runs = []
while time.perf_counter() - t_start < 5.0:
ev_s.record(); A @ B; ev_e.record(); ev_e.synchronize()
runs.append(dict(t=time.perf_counter() - t_start,
tflops=2 * n**3 / (ev_s.elapsed_time(ev_e) / 1e3) / 1e12))
stop.set(); poller.join()
runs, polled = pd.DataFrame(runs), pd.DataFrame(polled)
RESULTS["power_cap"] = {"runs": runs.to_dict("records"), "polled": polled.to_dict("records")}
del A, B; free_gpu()
first = runs[runs.t < 0.1].tflops.median()
last = runs[runs.t > 4.0].tflops.median()
print(f"first 100 ms: {first:.1f} TFLOP/s · after 4 s: {last:.1f} TFLOP/s "
f"({len(runs)} multiplications)")
print(polled.iloc[[0, len(polled) // 4, len(polled) // 2, -1]].to_string(index=False))
fig, ax = plt.subplots(figsize=(7, 3.2))
ax.plot(runs.t, runs.tflops, color=BLUE, lw=1, label="fp16 matmul, TFLOP/s")
ax.set_xlabel("seconds under load"); ax.set_ylabel("TFLOP/s", color=BLUE)
ax2 = ax.twinx()
ax2.plot(polled.t, polled["clocks.sm"], color=EMBER, lw=1, label="SM clock, MHz")
ax2.set_ylabel("SM clock (MHz)", color=EMBER); ax2.spines["right"].set_visible(True)
plt.tight_layout(); plt.show()first 100 ms: 20.8 TFLOP/s · after 4 s: 19.7 TFLOP/s (715 multiplications)
clocks.sm power.draw temperature.gpu t
585.0 31.49 68.0 0.047024
570.0 66.56 70.0 1.225221
555.0 66.85 71.0 2.505841
540.0 66.42 72.0 4.966242
[power_cap: 15.2 s]The expectation was wrong: there was no boost to lose. The first multiplication already ran at about 20 TFLOP/s, and within 150 ms the power was at the limit and the clock around 600 MHz. Yet in the operation sweep of Part 5, short FP16 multiplications reach 40–42 TFLOP/s — which at this card’s tensor-core rate needs the SMs above 1 GHz while they run. So something about what ran before a burst decides its clock.
Where the burst clock comes from
nvidia-smi samples the clock ten times a second at best, far too coarsely for a burst that lasts milliseconds. So the GPU measures itself: after every multiplication in a burst, a one-thread clock_probe kernel spins for 20 µs of wall time (read from the GPU’s own nanosecond timer) and counts the SM cycles that passed. Cycles over nanoseconds is the SM clock at that moment. The same burst — 2048³ FP16 multiplications back to back — is run after two different preludes: ten seconds of idling, and three seconds of memory-bound streaming, which keeps the card busy at a high clock while drawing far less power than the tensor cores do.
with experiment("burst_clock"):
n, N = 2048, 1500
A = torch.randn(n, n, device=DEV, dtype=torch.float16)
B = torch.randn(n, n, device=DEV, dtype=torch.float16)
probe = cp.zeros(3 * N, dtype=cp.int64)
stream_buf = cp.random.random(256 * 2**20 // 4, dtype=cp.float32)
stream_sink = cp.zeros(1, dtype=cp.float32)
A @ B; torch.cuda.synchronize()
def burst(prelude):
if prelude == "after 10 s idle":
time.sleep(10)
else:
t0 = time.perf_counter()
while time.perf_counter() - t0 < 3.0:
launch("stream_read_4", FULL_BLOCKS, 256, stream_buf, np.int64(stream_buf.size // 4), stream_sink)
cp.cuda.Device().synchronize()
before = smi("clocks.sm,power.draw")
evs = []
for i in range(N):
s, e = torch.cuda.Event(enable_timing=True), torch.cuda.Event(enable_timing=True)
s.record(); A @ B; e.record()
launch("clock_probe", 1, 1, probe, np.int32(i), np.int32(20_000))
evs.append((s, e))
torch.cuda.synchronize()
pr = probe.get().reshape(N, 3)
df = pd.DataFrame({"t_ms": (pr[:, 0] - pr[0, 0]) / 1e6, "sm_mhz": pr[:, 2] / pr[:, 1] * 1e3,
"tflops": [2 * n**3 / (s.elapsed_time(e) / 1e3) / 1e12 for s, e in evs]})
return df, before
runs = {}
for prelude in ["after 10 s idle", "after 3 s of memory-bound work"]:
df, before = burst(prelude)
runs[prelude] = {"trace": df.to_dict("records"), "smi_before": before}
first, last = df[df.t_ms < 5], df[df.t_ms > df.t_ms.max() - 100]
f_mhz, l_mhz = first.sm_mhz.median(), last.sm_mhz.median()
mid = (f_mhz + l_mhz) / 2
crossed = df[(df.sm_mhz < mid) if f_mhz > l_mhz else (df.sm_mhz > mid)]
runs[prelude]["halfway_ms"] = float(crossed.t_ms.iloc[0]) if len(crossed) else None
print(f"{prelude}: nvidia-smi just before — {before.get('clocks.sm', 0):.0f} MHz, "
f"{before.get('power.draw', 0):.0f} W")
print(f" first 5 ms: {f_mhz:6.0f} MHz, {first.tflops.median():5.1f} TFLOP/s")
print(f" last 100 ms: {l_mhz:6.0f} MHz, {last.tflops.median():5.1f} TFLOP/s")
if len(crossed) and abs(f_mhz - l_mhz) > 100:
print(f" the clock is halfway {'down' if f_mhz > l_mhz else 'up'} after {crossed.t_ms.iloc[0]:.1f} ms")
RESULTS["burst_clock"] = runs
del A, B, stream_buf; free_gpu()
fig, axes = plt.subplots(1, 2, figsize=(9, 3.4), sharey=True)
for ax, (prelude, r) in zip(axes, runs.items()):
d = pd.DataFrame(r["trace"])
ax.plot(d.t_ms, d.sm_mhz, color=EMBER, lw=1)
ax.set_title(prelude, fontsize=9); ax.set_xlabel("ms into the burst")
ax2 = ax.twinx(); ax2.plot(d.t_ms, d.tflops, color=BLUE, lw=0.8, alpha=0.7)
ax2.set_ylim(0, 50); ax2.spines["right"].set_visible(True)
axes[0].set_ylabel("SM clock (MHz)", color=EMBER)
plt.tight_layout(); plt.show()after 10 s idle: nvidia-smi just before — 585 MHz, 32 W
first 5 ms: 580 MHz, 18.1 TFLOP/s
last 100 ms: 600 MHz, 18.5 TFLOP/s
after 3 s of memory-bound work: nvidia-smi just before — 1245 MHz, 67 W
first 5 ms: 1409 MHz, 41.4 TFLOP/s
last 100 ms: 585 MHz, 18.1 TFLOP/s
the clock is halfway down after 37.3 ms
[burst_clock: 16.4 s]That settles it. After idling, the card starts at its idle clock and stays near it — tensor-core work there already draws close to the limit. After memory-bound work, it starts at the high clock that work left behind, delivers over 40 TFLOP/s, and the power controller pulls the clock down with a response of a few tens of milliseconds. Short multiplications that fit in that window are the 40 TFLOP/s bursts of Part 5.
Part 5 — the roofline
One kernel, one knob
intensity reads a float4, runs k FMAs on each of its four numbers, and writes it back: 8k FLOPs for 32 bytes, an arithmetic intensity of k/4 FLOP per byte. Nothing else changes. At small k the kernel waits on memory and the extra arithmetic is free; past some k the arithmetic becomes the bottleneck and every extra FMA costs time. Where the two regimes meet is the ridge point, and it should sit at peak FLOP/s ÷ peak bandwidth.
with experiment("roofline_knob"):
n4 = 2**24 # 256 MB in, 256 MB out
xi = cp.random.random(n4 * 4, dtype=cp.float32)
yi = cp.empty_like(xi)
rows = []
for k in [0, 1, 2, 4, 8, 16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768, 1024]:
t = gpu_time(lambda: launch("intensity", FULL_BLOCKS * 4, 256, xi, yi, np.int64(n4), np.int32(k)),
iters=7)
flops, nbytes = 8 * k * n4, 32 * n4
rows.append(dict(k=k, ai=flops / nbytes, ms=t * 1e3, flops_per_s=flops / t, bytes_per_s=nbytes / t))
knob = pd.DataFrame(rows)
del xi, yi; free_gpu()
bw_copy = knob[knob.ai <= 1].bytes_per_s.max() # this kernel reads *and* writes
bw_read = RESULTS.get("latency_hiding", {}).get("bw_max", bw_copy) # best read-only stream
fp32_roof = RESULTS["compute_roof"]["fma_peak"]["flops_per_s"]
fp16_sustained = max(v["flops_per_s"] for k, v in RESULTS["compute_roof"]["gemm"].items() if "fp16" in k)
roof = {"bw": max(bw_read, bw_copy), "bw_copy": bw_copy, "bw_read": bw_read, "fp32": fp32_roof,
"fp16_tensor_sustained": fp16_sustained,
"ridge_fp32_knob": fp32_roof / bw_copy}
RESULTS["roofline"] = {"roofs": roof, "knob": knob.to_dict("records")}
print(f"memory roof: {bw_read / 1e9:.0f} GB/s read-only, {bw_copy / 1e9:.0f} GB/s read+write · "
f"fp32 roof {fp32_roof / 1e12:.2f} TFLOP/s · fp16 tensor roof (sustained) {fp16_sustained / 1e12:.1f} TFLOP/s")
print(f"ridge point of this kernel, predicted: {roof['ridge_fp32_knob']:.1f} FLOP/B (fp32 roof ÷ read+write bandwidth)")
print(knob.to_string(index=False, float_format=lambda v: f"{v:.4g}"))memory roof: 275 GB/s read-only, 229 GB/s read+write · fp32 roof 6.81 TFLOP/s · fp16 tensor roof (sustained) 20.8 TFLOP/s
ridge point of this kernel, predicted: 29.7 FLOP/B (fp32 roof ÷ read+write bandwidth)
k ai ms flops_per_s bytes_per_s
0 0 2.405 0 2.232e+11
1 0.25 2.345 5.723e+10 2.289e+11
2 0.5 2.349 1.143e+11 2.286e+11
4 1 2.341 2.293e+11 2.293e+11
8 2 2.343 4.583e+11 2.292e+11
16 4 2.334 9.202e+11 2.3e+11
24 6 2.324 1.386e+12 2.31e+11
32 8 2.307 1.862e+12 2.327e+11
48 12 2.312 2.787e+12 2.322e+11
64 16 2.282 3.765e+12 2.353e+11
96 24 3.596 3.583e+12 1.493e+11
128 32 4.483 3.832e+12 1.198e+11
192 48 5.684 4.534e+12 9.445e+10
256 64 6.998 4.91e+12 7.672e+10
384 96 9.844 5.236e+12 5.454e+10
512 128 12.55 5.474e+12 4.276e+10
768 192 18.29 5.636e+12 2.935e+10
1024 256 23.74 5.788e+12 2.261e+10
[roofline_knob: 1.2 s]Real operations on the measured roofline
The same roofs, and on them the operations a model is actually made of. For each, the FLOPs and the compulsory bytes — each input read once, each output written once — are counted by hand, so the point’s x-coordinate is the best arithmetic intensity the operation could possibly have. (For softmax and layer norm the FLOP count is a convention — about 5 and 8 per element — since exponentials and square roots are not multiply-adds.)
with experiment("op_zoo"):
torch.backends.cudnn.benchmark = True
F = torch.nn.functional
f16, f32 = torch.float16, torch.float32
T = lambda *shape, dtype=f16: torch.randn(*shape, device=DEV, dtype=dtype)
zoo = []
def op(name, family, fn, flops, nbytes):
t = gpu_time(fn)
zoo.append(dict(name=name, family=family, seconds=t, flops=flops, bytes=nbytes,
ai=flops / nbytes, flops_per_s=flops / t, bytes_per_s=nbytes / t))
N = 2**26
a32, b32 = T(N, dtype=f32), T(N, dtype=f32)
op("x + y fp32", "element-wise", lambda: a32 + b32, N, 12 * N)
op("sum(x) fp32", "reduction", lambda: a32.sum(), N, 4 * N)
del a32, b32
a16, b16, c16 = T(N), T(N), T(N)
op("x + y fp16", "element-wise", lambda: a16 + b16, N, 6 * N)
op("x·y + z fp16", "element-wise", lambda: torch.addcmul(c16, a16, b16), 2 * N, 8 * N)
op("relu(x) fp16", "element-wise", lambda: torch.relu(a16), N, 4 * N)
del a16, b16, c16
R, C = 8192, 4096
xs = T(R, C)
op("softmax fp16", "normalization", lambda: torch.softmax(xs, -1), 5 * R * C, 4 * R * C)
op("layer norm fp16", "normalization", lambda: F.layer_norm(xs, (C,)), 8 * R * C, 4 * R * C)
del xs
for n in [256, 512, 1024, 2048, 4096, 8192]:
A, B = T(n, n), T(n, n)
op(f"matmul {n}³ fp16", "matmul", lambda: A @ B, 2 * n**3, 3 * 2 * n * n)
del A, B
for n in [1024, 4096]:
A, B = T(n, n, dtype=f32), T(n, n, dtype=f32)
op(f"matmul {n}³ fp32", "matmul fp32", lambda: A @ B, 2 * n**3, 3 * 4 * n * n)
del A, B
d = 4096
W = T(d, d)
for b in [1, 4, 16, 64, 256, 1024]:
X = T(b, d)
op(f"linear 4096→4096, batch {b} fp16", "linear, batch sweep", lambda: X @ W,
2 * b * d * d, 2 * (d * d + 2 * b * d))
del W, X
xc = T(32, 64, 56, 56).to(memory_format=torch.channels_last)
wc = T(64, 64, 3, 3).to(memory_format=torch.channels_last)
op("conv 3×3, 64→64, 56², batch 32 fp16", "convolution", lambda: F.conv2d(xc, wc, padding=1),
2 * 32 * 64 * 56 * 56 * 64 * 9, 2 * (2 * xc.numel() + wc.numel()))
del xc, wc
Bq, H, L, D = 8, 16, 1024, 64
q, k_, v = T(Bq, H, L, D), T(Bq, H, L, D), T(Bq, H, L, D)
op("attention, L=1024, d=64 fp16", "attention", lambda: F.scaled_dot_product_attention(q, k_, v),
4 * Bq * H * L * L * D, 2 * 4 * Bq * H * L * D)
del q, k_, v
free_gpu()
zoo = pd.DataFrame(zoo)
r = RESULTS["roofline"]["roofs"]
# Short fp16 GEMMs finish before the power limit bites, so they show a second, higher roof:
# the best tensor-core throughput any operation reached in these short bursts.
r["fp16_tensor_burst"] = float(zoo[~zoo.name.str.contains("fp32")].flops_per_s.max())
r["ridge_fp16_burst"] = r["fp16_tensor_burst"] / r["bw"]
r["ridge_fp16_sustained"] = r["fp16_tensor_sustained"] / r["bw"]
compute_roof = [r["fp32"] if "fp32" in nm else r["fp16_tensor_burst"] for nm in zoo.name]
zoo["roof_flops_per_s"] = np.minimum(compute_roof, r["bw"] * zoo.ai)
zoo["pct_of_roof"] = 100 * zoo.flops_per_s / zoo.roof_flops_per_s
zoo["bound"] = np.where(r["bw"] * zoo.ai < np.array(compute_roof), "memory", "compute")
RESULTS["op_zoo"] = zoo.to_dict("records")
print(f"fp16 tensor roof: {r['fp16_tensor_burst'] / 1e12:.1f} TFLOP/s in short bursts, "
f"{r['fp16_tensor_sustained'] / 1e12:.1f} sustained · memory roof {r['bw'] / 1e9:.0f} GB/s")
print(f"ridge points (fp16): {r['ridge_fp16_burst']:.0f} FLOP/B burst, {r['ridge_fp16_sustained']:.0f} sustained\n")
print(zoo[["name", "ai", "seconds", "flops_per_s", "bytes_per_s", "pct_of_roof", "bound"]]
.to_string(index=False, float_format=lambda v: f"{v:.3g}"))fp16 tensor roof: 40.3 TFLOP/s in short bursts, 20.8 sustained · memory roof 275 GB/s
ridge points (fp16): 146 FLOP/B burst, 76 sustained
name ai seconds flops_per_s bytes_per_s pct_of_roof bound
x + y fp32 0.0833 0.00337 1.99e+10 2.39e+11 86.8 memory
sum(x) fp32 0.25 0.000992 6.76e+10 2.71e+11 98.3 memory
x + y fp16 0.167 0.00166 4.04e+10 2.42e+11 88 memory
x·y + z fp16 0.25 0.00215 6.25e+10 2.5e+11 90.8 memory
relu(x) fp16 0.25 0.00116 5.81e+10 2.32e+11 84.4 memory
softmax fp16 1.25 0.00075 2.24e+11 1.79e+11 65 memory
layer norm fp16 2 0.000778 3.45e+11 1.72e+11 62.7 memory
matmul 256³ fp16 85.3 3.39e-05 9.89e+11 1.16e+10 4.21 memory
matmul 512³ fp16 171 3.72e-05 7.21e+12 4.23e+10 17.9 compute
matmul 1024³ fp16 341 7.58e-05 2.83e+13 8.3e+10 70.3 compute
matmul 2048³ fp16 683 0.000426 4.03e+13 5.9e+10 100 compute
matmul 4096³ fp16 1.37e+03 0.00813 1.69e+13 1.24e+10 41.9 compute
matmul 8192³ fp16 2.73e+03 0.0551 2e+13 7.31e+09 49.6 compute
matmul 1024³ fp32 171 0.000922 2.33e+12 1.36e+10 34.2 compute
matmul 4096³ fp32 683 0.0418 3.29e+12 4.82e+09 48.3 compute
linear 4096→4096, batch 1 fp16 1 0.000188 1.78e+11 1.78e+11 64.7 memory
linear 4096→4096, batch 4 fp16 3.99 0.000207 6.49e+11 1.63e+11 59.1 memory
linear 4096→4096, batch 16 fp16 15.9 0.000229 2.34e+12 1.47e+11 53.6 memory
linear 4096→4096, batch 64 fp16 62.1 0.000215 1e+13 1.61e+11 58.6 memory
linear 4096→4096, batch 256 fp16 228 0.000478 1.8e+13 7.9e+10 44.6 compute
linear 4096→4096, batch 1024 fp16 683 0.00164 2.1e+13 3.07e+10 52 compute
conv 3×3, 64→64, 56², batch 32 fp16 287 0.000307 2.41e+13 8.39e+10 59.8 compute
attention, L=1024, d=64 fp16 512 0.00214 1.61e+13 3.14e+10 39.9 compute
[op_zoo: 3.2 s]with experiment("roofline_plot"):
r = RESULTS["roofline"]["roofs"]
knob = pd.DataFrame(RESULTS["roofline"]["knob"])
zoo = pd.DataFrame(RESULTS["op_zoo"])
ai = np.logspace(-2, 4, 300)
fig, ax = plt.subplots(figsize=(7.5, 4.4))
ax.plot(ai, np.minimum(r["fp16_tensor_burst"], r["bw"] * ai) / 1e12, color=INK, lw=1.2,
label="fp16 tensor-core roof, short bursts")
ax.plot(ai, np.minimum(r["fp16_tensor_sustained"], r["bw"] * ai) / 1e12, color=INK, lw=1.0, ls=":",
label="fp16 tensor-core roof, sustained (power-capped)")
ax.plot(ai, np.minimum(r["fp32"], r["bw"] * ai) / 1e12, color=GRAY, lw=1.2, ls="--", label="fp32 CUDA-core roof")
k = knob[knob.k > 0]
ax.plot(k.ai, k.flops_per_s / 1e12, "-o", color=FOREST, ms=3, label="one-knob kernel (fp32)")
fam_colors = {"element-wise": EMBER, "reduction": EMBER, "normalization": PLUM, "matmul": BLUE,
"matmul fp32": GRAY, "linear, batch sweep": BLUE, "convolution": FOREST, "attention": PLUM}
for fam, d in zoo.groupby("family"):
marker = "s" if "linear" in fam else "o"
ax.scatter(d.ai, d.flops_per_s / 1e12, s=18, color=fam_colors.get(fam, INK), marker=marker,
label=fam, zorder=3)
ax.set_xscale("log"); ax.set_yscale("log")
ax.set_xlabel("arithmetic intensity (FLOP / byte)"); ax.set_ylabel("TFLOP/s")
ax.legend(frameon=False, fontsize=7, loc="lower right", ncol=2)
plt.tight_layout(); plt.show()
[roofline_plot: 0.7 s]Why attention sits so far below its roof
scaled_dot_product_attention is not one kernel but a dispatcher over several implementations — FlashAttention, memory-efficient attention, cuDNN, and a plain “math” fallback that materializes the full L × L score matrix. Each is forced in turn here; the ones this GPU can’t run say why.
with experiment("attention_backends"):
import warnings
from torch.nn.attention import sdpa_kernel, SDPBackend
F = torch.nn.functional
Bq, H, L, D = 8, 16, 1024, 64
q, k_, v = (torch.randn(Bq, H, L, D, device=DEV, dtype=torch.float16) for _ in range(3))
flops = 4 * Bq * H * L * L * D
rows = [dict(backend="default (PyTorch's choice)",
seconds=gpu_time(lambda: F.scaled_dot_product_attention(q, k_, v)), why="")]
for be in [SDPBackend.FLASH_ATTENTION, SDPBackend.EFFICIENT_ATTENTION, SDPBackend.CUDNN_ATTENTION,
SDPBackend.MATH]:
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter("always")
try:
with sdpa_kernel(be):
t = gpu_time(lambda: F.scaled_dot_product_attention(q, k_, v))
rows.append(dict(backend=be.name, seconds=t, why=""))
except Exception as e:
first_line = lambda m: (str(m).splitlines() or [""])[0][:140]
reasons = [first_line(w.message) for w in caught]
rows.append(dict(backend=be.name, seconds=None,
why=" | ".join(reasons) or first_line(e)))
att = pd.DataFrame(rows)
att["tflops"] = flops / att.seconds / 1e12
RESULTS["attention_backends"] = att.to_dict("records")
del q, k_, v; free_gpu()
with pd.option_context("display.max_colwidth", 140):
print(att.to_string(index=False, float_format=lambda x: f"{x:.3g}")) backend seconds why tflops
default (PyTorch's choice) 0.00191 18
FLASH_ATTENTION NaN Memory efficient kernel not used because: (Triggered internally at /pytorch/aten/src/ATen/native/transformers/cuda/sdp_utils.cpp:960.) | Memory Efficient attention has been runtime disabled. (Triggered internally at /pytorch/aten/src/ATen/native/transformers/sdp_utils_cpp.h:55 | Flash attention kernel not used because: (Triggered internally at /pytorch/aten/src/ATen/native/transformers/cuda/sdp_utils.cpp:962.) | Flash attention only supports gpu architectures in the range [sm80, sm121]. Attempting to run on a sm 7.5 gpu. (Triggered internally at /pyt | cuDNN attention kernel not used because: (Triggered internally at /pytorch/aten/src/ATen/native/transformers/cuda/sdp_utils.cpp:964.) | cuDNN attention has been runtime disabled. (Triggered internally at /pytorch/aten/src/ATen/native/transformers/cuda/sdp_utils.cpp:680.) NaN
EFFICIENT_ATTENTION 0.00241 14.3
CUDNN_ATTENTION NaN Memory efficient kernel not used because: (Triggered internally at /pytorch/aten/src/ATen/native/transformers/cuda/sdp_utils.cpp:960.) | Memory Efficient attention has been runtime disabled. (Triggered internally at /pytorch/aten/src/ATen/native/transformers/sdp_utils_cpp.h:55 | Flash attention kernel not used because: (Triggered internally at /pytorch/aten/src/ATen/native/transformers/cuda/sdp_utils.cpp:962.) | Flash attention has been runtime disabled. (Triggered internally at /pytorch/aten/src/ATen/native/transformers/sdp_utils_cpp.h:540.) | cuDNN attention kernel not used because: (Triggered internally at /pytorch/aten/src/ATen/native/transformers/cuda/sdp_utils.cpp:964.) | cuDNN MHA only supports gpu architectures in the range [sm80, sm121]. Attempting to run on a sm 7.5 gpu. (Triggered internally at /pytorch/a NaN
MATH 0.0305 1.13
[attention_backends: 0.8 s]FlashAttention and cuDNN’s attention both refuse: they need sm80 or newer, and the T4 is sm75. PyTorch’s default is the memory-efficient kernel — the two times match — and the math fallback, which writes the whole score matrix to memory, is an order of magnitude slower.
The cheapest byte is the one never moved
Three element-wise operations in a row — multiply, add, ReLU — are three kernels, and each reads its input from DRAM and writes its output back. Fused into one kernel, the intermediate values stay in registers: 4 bytes per element instead of 12 (FP16), for the same FLOPs. If the chain is memory-bound, the roofline predicts the fused version is three times faster. torch.compile does the fusing.
with experiment("fusion"):
x = torch.randn(2**26, device=DEV, dtype=torch.float16)
chain = lambda t: torch.relu(t * 1.5 + 0.5)
fused = torch.compile(chain)
fused(x); torch.cuda.synchronize() # compile outside the timing
t_eager, t_fused = gpu_time(lambda: chain(x)), gpu_time(lambda: fused(x))
RESULTS["fusion"] = {"eager_s": t_eager, "fused_s": t_fused, "speedup": t_eager / t_fused,
"bytes_ratio_predicted": 12 / 4, "elements": x.numel()}
print(f"eager (3 kernels) {t_eager * 1e3:.3f} ms · fused {t_fused * 1e3:.3f} ms · "
f"{t_eager / t_fused:.2f}× faster (bytes moved: 3× fewer)")
del x; free_gpu()eager (3 kernels) 3.322 ms · fused 1.184 ms · 2.81× faster (bytes moved: 3× fewer)
[fusion: 17.0 s]Part 6 — the paradox, predicted
The fourth number, taken apart
Part 1’s slowest result was addition “including the copies”: the GPU lost by 10× at every size. Here the same naive line — torch.add(x.to("cuda"), y.to("cuda")).cpu() — is split into its stages and each is timed on its own: copying x in, copying y in, the addition, and copying the result back into a brand-new CPU tensor. Two references sit next to it. One is the cost of what .cpu() does besides copying: allocating a fresh host tensor and touching its memory for the first time. The other is the best the link allows: pinned buffers allocated once and reused, so nothing but the bytes themselves crosses PCIe.
with experiment("copy_path"):
rows = []
for N in [2**10, 2**16, 2**20, 2**24, 2**28]:
x_c, y_c, o_c = torch.randn(N), torch.randn(N), torch.empty(N)
cpu_s = wall_time(lambda: torch.add(x_c, y_c, out=o_c))
def naive_stages():
torch.cuda.synchronize()
t0 = time.perf_counter(); xg = x_c.to(DEV); torch.cuda.synchronize()
t1 = time.perf_counter(); yg = y_c.to(DEV); torch.cuda.synchronize()
t2 = time.perf_counter(); zg = xg + yg; torch.cuda.synchronize()
t3 = time.perf_counter(); zc = zg.cpu()
t4 = time.perf_counter()
return np.array([t1 - t0, t2 - t1, t3 - t2, t4 - t3])
naive_stages()
stages = np.median([naive_stages() for _ in range(3 if N >= 2**24 else 15)], axis=0)
alloc_touch_s = wall_time(lambda: torch.empty(N).zero_())
xp, yp, zp = x_c.pin_memory(), y_c.pin_memory(), torch.empty(N).pin_memory()
xg, yg, zg = (torch.empty(N, device=DEV) for _ in range(3))
def best():
xg.copy_(xp, non_blocking=True); yg.copy_(yp, non_blocking=True)
torch.add(xg, yg, out=zg); zp.copy_(zg, non_blocking=True)
best_s = wall_time(best, sync=True)
rows.append(dict(elements=N, cpu_s=cpu_s, h2d_x_s=stages[0], h2d_y_s=stages[1], add_s=stages[2],
d2h_s=stages[3], naive_s=float(stages.sum()), alloc_touch_s=alloc_touch_s,
best_s=best_s))
del x_c, y_c, o_c, xp, yp, zp, xg, yg, zg
free_gpu()
cpath = pd.DataFrame(rows)
cpath["h2d_gb_s"] = 4 * cpath.elements / cpath.h2d_x_s / 1e9
cpath["d2h_gb_s"] = 4 * cpath.elements / cpath.d2h_s / 1e9
cpath["naive_vs_cpu"] = cpath.cpu_s / cpath.naive_s
cpath["best_vs_cpu"] = cpath.cpu_s / cpath.best_s
RESULTS["copy_path"] = cpath.to_dict("records")
print(cpath.to_string(index=False, float_format=lambda v: f"{v:.3g}")) elements cpu_s h2d_x_s h2d_y_s add_s d2h_s naive_s alloc_touch_s best_s h2d_gb_s d2h_gb_s naive_vs_cpu best_vs_cpu
1024 2.66e-06 3.65e-05 3.29e-05 3.19e-05 2.48e-05 0.000126 3.42e-06 4.87e-05 0.112 0.165 0.0211 0.0546
65536 1.09e-05 0.000117 0.000111 3.92e-05 9.92e-05 0.000366 9.53e-06 9.92e-05 2.24 2.64 0.0298 0.11
1048576 0.000291 0.0011 0.0011 9.86e-05 0.00424 0.00654 0.000111 0.00114 3.8 0.989 0.0445 0.255
16777216 0.00995 0.0164 0.0154 0.000907 0.0682 0.101 0.0284 0.0169 4.08 0.984 0.0986 0.588
268435456 0.143 0.234 0.237 0.0131 0.963 1.45 0.452 0.27 4.59 1.12 0.0987 0.529
[copy_path: 20.8 s]For the largest size, the naive line spends almost all its time on the road: the addition is a few milliseconds, the inputs take about half a second to go in through pageable memory, and the result takes about a second to come back, half of which is the cost of allocating and first touching a fresh host tensor. With pinned buffers allocated once, the round trip shrinks to what the link allows — 12 bytes per element at ~12 GB/s — and even then it is about twice as slow as the CPU adding the vectors in its own memory.
The model
Each machine is now reduced to three measured numbers per operation: a compute roof P, a memory roof B, and a floor t0 — the fixed cost of getting the operation started. The model is t = t0 + max (FLOPs/P, bytes/B). The GPU’s roofs come from Parts 4 and 5 (the FP32 FMA roof, and the read+write bandwidth of the one-knob kernel, since addition writes as much as it reads); the CPU’s from the largest sizes of Part 1. The floor is the one number taken from the Part 1 curves themselves: each operation’s time at its smallest size. Nothing else is fitted. The question is how much of the paradox this explains.
with experiment("paradox_model"):
d = pd.DataFrame(RESULTS["paradox"])
mm, ad = d[d.op == "matmul"], d[d.op == "add"]
roofs = RESULTS["roofline"]["roofs"]
cpu = {"P": float((mm.flops / mm.cpu_s).max()), "B": float((ad["bytes"] / ad.cpu_s).iloc[-1])}
gpu = {"P": roofs["fp32"], "B": roofs["bw_copy"]}
floor = {(op, dev): float(d[d.op == op].sort_values("size")[col].iloc[0])
for op in ("matmul", "add") for dev, col in (("cpu", "cpu_s"), ("gpu", "gpu_s"))}
model = lambda m, t0, F_, B_: t0 + max(F_ / m["P"], B_ / m["B"])
d["cpu_model_s"] = [model(cpu, floor[(o, "cpu")], f, b) for o, f, b in zip(d.op, d.flops, d["bytes"])]
d["gpu_model_s"] = [model(gpu, floor[(o, "gpu")], f, b) for o, f, b in zip(d.op, d.flops, d["bytes"])]
d["speedup_model"] = d.cpu_model_s / d.gpu_model_s
RESULTS["paradox_model"] = {"cpu": cpu, "gpu": gpu, "floors": {f"{o}/{m}": v for (o, m), v in floor.items()},
"rows": d[["op", "size", "speedup", "speedup_model"]].to_dict("records")}
print(f"CPU: P = {cpu['P'] / 1e9:.0f} GFLOP/s, B = {cpu['B'] / 1e9:.1f} GB/s")
print(f"GPU: P = {gpu['P'] / 1e12:.2f} TFLOP/s, B = {gpu['B'] / 1e9:.0f} GB/s")
print("floors: " + ", ".join(f"{o} on {m} {v * 1e6:.1f} µs" for (o, m), v in floor.items()))
print(f"=> asymptotic speed-up: compute-bound {gpu['P'] / cpu['P']:.0f}×, memory-bound {gpu['B'] / cpu['B']:.0f}×")
print(d[["op", "size", "speedup", "speedup_model"]].to_string(index=False, float_format=lambda v: f"{v:.3g}"))
fig, axes = plt.subplots(1, 2, figsize=(9, 3.4), sharey=True)
for ax, op in zip(axes, ["matmul", "add"]):
s = d[d.op == op]
ax.plot(s["size"], s.speedup, "o", color=BLUE, ms=4, label="measured")
ax.plot(s["size"], s.speedup_model, "-", color=INK, lw=1, label="three-number model")
ax.axhline(1, color=GRAY, lw=0.8, ls="--")
ax.set_xscale("log", base=2); ax.set_yscale("log"); ax.set_title(op)
axes[0].set_ylabel("GPU speed-up over the CPU (×)"); axes[0].legend(frameon=False, fontsize=8)
plt.tight_layout(); plt.show()CPU: P = 215 GFLOP/s, B = 22.5 GB/s
GPU: P = 6.81 TFLOP/s, B = 229 GB/s
floors: matmul on cpu 6.2 µs, matmul on gpu 34.8 µs, add on cpu 2.7 µs, add on gpu 20.5 µs
=> asymptotic speed-up: compute-bound 32×, memory-bound 10×
op size speedup speedup_model
matmul 16 0.179 0.183
matmul 24 0.185 0.188
matmul 32 0.189 0.195
matmul 48 0.228 0.214
matmul 64 0.405 0.248
matmul 96 0.531 0.41
matmul 128 0.866 0.721
matmul 192 1.85 1.95
matmul 256 2.54 4.08
matmul 384 6.98 10.3
matmul 512 8.6 16.9
matmul 768 11.9 25.1
matmul 1024 12.4 28.5
matmul 1536 12.8 30.6
matmul 2048 13 31.2
matmul 3072 18.4 31.5
matmul 4096 18.9 31.6
add 1024 0.132 0.158
add 2048 0.132 0.184
add 4096 0.147 0.236
add 8192 0.173 0.338
add 16384 0.231 0.536
add 32768 0.375 0.909
add 65536 0.494 1.57
add 131072 0.847 2.66
add 262144 3.3 4.17
add 524288 5.84 5.89
add 1048576 4.6 7.46
add 2097152 5.55 8.61
add 4194304 8.98 9.33
add 8388608 9.9 9.74
add 16777216 10.6 9.96
add 33554432 10.9 10.1
add 67108864 10.9 10.1
add 134217728 10.9 10.2
add 268435456 11 10.2
[paradox_model: 0.7 s]The machine in one table
with experiment("summary"):
lat = pd.DataFrame(RESULTS["latency"]); hbw = pd.DataFrame(RESULTS["hierarchy_bandwidth"])
pc = pd.DataFrame(RESULTS["pcie"])
L2 = spec["l2_bytes"]
pick = lambda df, lo, hi, col: float(df[(df["bytes"] >= lo) & (df["bytes"] <= hi)][col].median())
# capacities as the pointer chase sees them: the largest chain that still runs at a level's speed
l2_ns, dram_ns = pick(lat, L2 / 8, L2 / 2, "ns"), pick(lat, 4 * L2, 16 * L2, "ns")
l1_cap = int(lat[lat.ns <= 1.3 * lat.ns.min()]["bytes"].max())
l2_cap = int(lat[lat.ns <= (l2_ns + dram_ns) / 2]["bytes"].max())
caps_kb = RESULTS.get("l1_capacity", {}).get("capacity_kb", {})
RESULTS["capacities"] = {"l1_bytes": l1_cap, "l2_bytes": l2_cap, "l1_by_carveout_kb": caps_kb}
l1_size = (" · ".join(f"{v} KB {k}" for k, v in caps_kb.items()) if caps_kb
else f"{l1_cap // 1024} KB / SM measured")
table = [
("registers", f"{spec['regs_per_sm'] * 4 // 1024} KB / SM", "≈ 1 cycle", "—"),
("L1 (data cache)", l1_size,
f"{pick(lat, 2**12, l1_cap, 'ns'):.0f} ns", f"{pick(hbw, 2**12, l1_cap // 2, 'gb_s'):.0f} GB/s"),
("L2", f"{L2 / 2**20:.0f} MB spec, {l2_cap / 2**20:.0f} MB measured", f"{l2_ns:.0f} ns",
f"{pick(hbw, L2 / 8, L2 / 2, 'gb_s'):.0f} GB/s"),
("DRAM", f"{spec['dram_bytes'] / 2**30:.0f} GB", f"{pick(lat, 4 * L2, 16 * L2, 'ns'):.0f} ns",
f"{pick(hbw, 64 * L2, 2**30, 'gb_s'):.0f} GB/s"),
("host RAM over PCIe", "—", f"{pc.h2d_pinned_s.iloc[0] * 1e6:.0f} µs (4-byte copy)",
f"{pc.h2d_pinned_gb_s.max():.1f} GB/s pinned"),
]
RESULTS["summary"] = table
print(pd.DataFrame(table, columns=["level", "size", "latency", "bandwidth"]).to_string(index=False))
print("\nerrors:", list(RESULTS["errors"]) or "none") level size latency bandwidth
registers 256 KB / SM ≈ 1 cycle —
L1 (data cache) 56 KB default · 56 KB max L1 (carveout 0) · 24 KB max shared (carveout 100) 27 ns 3870 GB/s
L2 4 MB spec, 4 MB measured 150 ns 1107 GB/s
DRAM 15 GB 314 ns 268 GB/s
host RAM over PCIe — 21 µs (4-byte copy) 12.3 GB/s pinned
errors: none
[summary: 0.0 s]RESULTS["meta"].update(
executed=datetime.datetime.now(datetime.timezone.utc).isoformat(timespec="seconds"),
torch=torch.__version__, cuda=torch.version.cuda, cupy=cp.__version__,
python=platform.python_version(), kaggle=bool(os.environ.get("KAGGLE_KERNEL_RUN_TYPE")),
)
out_dir = "/kaggle/working" if os.path.isdir("/kaggle/working") else "."
with open(os.path.join(out_dir, "gpu_results.json"), "w") as f:
json.dump(RESULTS, f, default=jsonable, indent=1)
print("saved", os.path.join(out_dir, "gpu_results.json"))saved /kaggle/working/gpu_results.json