P14 — Hardware-Aware ML System

Run it first. There is a companion page that builds this project's machinery as numbered, independently runnable blocks and then assembles them into one measured system: P14 hands-on — block by block (handson/h*.py). Every number on it was produced by running the code. Read it alongside the milestones below.

Medium · 77 hours · Weeks 111–117 · Stage 5 · C, with CUDA/Metal where available

Table of Contents

• The Loop, Instantiated
• Why This Project Matters
• Prerequisites
• Duration and Size
• Central Technical Questions
• The Measurements That Define This Project
• The Systolic Array, From First Principles
• Showcase — Place a Kernel Before You Write It
• Implementation Milestones
• Concepts To Study
• Primary-Source Readings
• Experiments
• Benchmarks and Metrics
• Correctness Tests
• Failure Tests
• Expected Difficulties
• Scope Boundaries
• Deliverables
• Exit Criteria
• Extension Ideas
• Connections
• References

The Loop, Instantiated

Why This Project Matters

You have spent thirteen projects learning that performance is data movement. This is where that becomes the whole subject.

The specific realisation this project is built around: a naive matrix multiply and an optimised one execute exactly the same number of arithmetic operations. Every difference in runtime — and you will measure a factor of 18× within your own C code, and nearly 900× against a hardware matrix unit — comes from where the operands were when they were needed. Once that is measured rather than believed, you understand the motivation for GPUs, TPUs, and every ML accelerator, and you can reason about new ones from their memory hierarchy rather than their FLOP/s number.

For your trajectory it is also the most direct route to reasoning about inference cost. tools/roofline.py already encodes the key result — decode arithmetic intensity equals batch size, so a batch-1 chatbot uses ~0.3% of an H100's multipliers — and this project is where that stops being a script and becomes understanding.

Prerequisites

• P13-II complete — its blocked kernel and measured GFLOP/s are the baseline here
• P12 helpful — TLB, page, and cache behaviour from the OS side
• C; willingness to read compiler output. From math.md: nothing new

Duration and Size

Medium, 77 hours, 7 weeks.

Central Technical Questions

1. Why is naive matmul so far below peak when it does exactly the right arithmetic?
2. What does blocking do, precisely, to DRAM traffic? Derive the factor.
3. What is arithmetic intensity, and how do you move a kernel along the roofline?
4. Why does a systolic array beat a general-purpose core at this one job? The answer is about operand reuse per fetch, not about clock speed.
5. What does reduced precision actually buy? Not faster arithmetic — smaller operands.
6. Why is inference batch size the dominant performance variable?

The Measurements That Define This Project

All produced on the machine used to build this track (12-core arm64 laptop). Reproduce each on yours in milestone 1; the numbers will differ, the shape will not.

CPU matmul, N=512, fp32, single-threaded C

Four results, each worth a paragraph in your report:

1. Loop reordering alone is 14.3× (1.91 → 27.33 at -O2), with no blocking and no intrinsics. The i,j,k order strides B by N — a cache miss per inner iteration — while i,k,j walks B contiguously. Same arithmetic, same instruction count, one-line change.
2. At -O2, blocking is worse than plain reordering (26.15 vs 27.33 at the best block size). Blocking only pays once the inner loop vectorises: at -O3 with -mcpu=native, BS=32 reaches 58.30 GFLOP/s, 2.1× over plain reordering. The lesson: an optimisation's value depends on what other optimisations are present, and measuring one in isolation can tell you it is useless when it is not.
3. Block size has a clear optimum at 32 and degrades either side. BS=16 wastes vector width and loop overhead; BS=128 overflows L1 (128²×4 bytes × 3 arrays = 192 KB). Compute your own L1 and L2 sizes and predict your optimum before sweeping.
4. Peak vs achieved. Apple's Accelerate BLAS on the same machine, same problem: 1,679 GFLOP/s at N=512 (verified across N=256–2048 at 1,008–1,705 GFLOP/s, with fp32-consistent relative error ~10⁻⁷). Your best hand-written kernel reaches 58.3 — 29× off. That gap is not sloppy coding: Accelerate dispatches to Apple's AMX matrix coprocessor, a dedicated hardware matrix unit. The 29× is this project's thesis, measured on your own laptop before you write the simulator.

The roofline placement

From tools/roofline.py, a 4096³ bf16 GEMM on an H100:

Same arithmetic. Same hardware. A 295× range decided entirely by operand reuse. This is the number to quote when someone asks why kernel engineering is a job.

Hardware ridge points

Check the asterisk. Vendor headline FLOP/s are routinely quoted with 2:4 structured sparsity, which does not apply to a dense GEMM and doubles the number. Using the sparse figure as a denominator halves every efficiency you report.

The Systolic Array, From First Principles

The problem. A CPU core computing \(C = AB\) fetches operands from the register file for every multiply-accumulate. Each fetch costs energy and bandwidth, and the register file has few ports, so the core is limited by operand delivery, not by multiplier count. Adding multipliers does not help; they starve.

The idea. Arrange \(k \times k\) multiply-accumulate cells in a grid. Each cell holds one weight. Activations flow horizontally, partial sums flow vertically. Each value entering the array is used by every cell in its row or column before leaving.

The consequence, quantified. For a \(k \times k\) array:

• Values fetched from memory per cycle: \(O(k)\) — one column of activations
• Multiply-accumulates performed per cycle: \(k^2\)
• Operand reuse: \(O(k)\)

At \(k = 256\) (the TPUv1 dimension), each fetched value participates in 256 multiply-accumulates. That is a 256× reduction in operand-delivery pressure, achieved by wiring rather than by caching — no tag comparisons, no misses, no replacement policy.

The array does \(k^2 = 65{,}536\) MACs per cycle. At 700 MHz that is \(2 \times 65{,}536 \times 700\times 10^6 = 91.75\) TOPS, which is the TPUv1's quoted 92 TOPS. You can derive the headline number of a real accelerator from two integers, and doing so in your report is worth more than any amount of description.

The costs, which your simulator must expose:

• Fill and drain. The pipeline takes \(\approx 2k\) cycles to fill and drain, so a small matrix wastes most of its time. Utilisation for an \(M \times N \times K\) GEMM on a \(k \times k\) array is roughly \(\frac{MNK}{k^2(\ldots + 2k)}\) — and measuring that curve is milestone 8.
• Shape mismatch. A matrix whose dimensions are not multiples of \(k\) leaves cells idle. Your simulator must report utilisation, and the utilisation drop on awkward shapes is the most instructive output it produces.
• It does one thing. No branches, no gather, no control flow. Total inflexibility is what buys the efficiency, and that trade — the same restriction-buys-automation trade as P06, now in silicon — is the report's closing argument.

Showcase — Place a Kernel Before You Write It

Fifteen minutes. Two lines of arithmetic decide whether a kernel is worth optimising and in which direction, and they work before any code exists.

# P14 -- place a kernel on the roofline before writing it.
HW = {"h100":(989.4e12, 3.35e12), "a100":(312e12, 2.039e12), "laptop":(1.7e12, 57.5e9)}
def place(name, W, Q):
    p,b = HW[name]; I=W/Q; ridge=p/b
    bound = "COMPUTE" if I>ridge else "MEMORY"
    t = max(W/p, Q/b)
    return I, ridge, bound, t
N=4096
for hw in HW:
    I,ridge,bound,t = place(hw, 2*N**3, (3*N*N)*2)      # perfect reuse, bf16
    print(f"{hw:<8} GEMM {N}^3: I={I:>7.1f}  ridge={ridge:>6.1f}  {bound:<8} {t*1e3:>7.3f} ms")
print()
for hw in HW:
    I,ridge,bound,t = place(hw, 2*7e9*1, 7e9*2)         # decode, batch 1
    print(f"{hw:<8} decode b=1: I={I:>7.1f}  ridge={ridge:>6.1f}  {bound:<8} "
          f"{1/t:>7.0f} tok/s")
print("\\nSame kernel, three machines, two different regimes. Optimising for the")
print("wrong one is why 'we upgraded the GPU and nothing improved' happens.")

h100     GEMM 4096^3: I= 1365.3  ridge= 295.3  COMPUTE    0.139 ms
a100     GEMM 4096^3: I= 1365.3  ridge= 153.0  COMPUTE    0.441 ms
laptop   GEMM 4096^3: I= 1365.3  ridge=  29.6  COMPUTE   80.846 ms

h100     decode b=1: I=    1.0  ridge= 295.3  MEMORY       239 tok/s
a100     decode b=1: I=    1.0  ridge= 153.0  MEMORY       146 tok/s
laptop   decode b=1: I=    1.0  ridge=  29.6  MEMORY         4 tok/s
\nSame kernel, three machines, two different regimes. Optimising for the
wrong one is why 'we upgraded the GPU and nothing improved' happens.

Identical arithmetic, three machines, two regimes. The GEMM is compute-bound everywhere; batch-1 decode is memory-bound everywhere, at 1 FLOP/byte against ridges of 30 to 295. Optimising the wrong one is why "we upgraded the GPU and nothing improved" is such a common sentence.

Implementation Milestones

If no GPU is available, milestone 9 is replaced by extending the simulator to model a multi-core vector machine and comparing all three architectures in simulation. The learning objective — understanding why the architectures differ — is preserved.

Concepts To Study

• The memory hierarchy: registers, L1/L2/L3, DRAM, HBM — capacity, latency, and bandwidth at each level, for your machine
• Cache mechanics: lines, associativity, replacement, prefetching; conflict misses and why power-of-two strides are pathological
• Arithmetic intensity and the roofline model
• Loop transformations: interchange, tiling, unrolling, and how each changes the access pattern
• SIMD: vector width, alignment, why the compiler often fails to vectorise, and how to check
• Systolic arrays: dataflow (weight-stationary, output-stationary, row-stationary), fill/drain, utilisation
• Reduced precision: fp16, bf16, fp8, int8; the exponent/mantissa trade; why bf16 won for training
• Quantization: symmetric/asymmetric, per-tensor/per-channel, calibration, accumulate-in-higher-precision
• Batching: why intensity in decode equals batch size
• Kernel fusion and operator scheduling
• Amdahl and the memory wall

Primary-Source Readings

Budget: 13 hours.

Experiments

E9 and E14 together are the project. E9 shows how array dimension trades against utilisation. E14 puts your hand-written CPU kernel (58 GFLOP/s measured), your platform's BLAS/AMX path (1,679 GFLOP/s measured), and your simulated systolic array on one chart, and answers the project's title question with three of your own numbers.

E7 has a trap worth falling into. int8 will show a large throughput gain and, on a real model, a possibly catastrophic accuracy loss — because transformer activations have outlier channels whose dynamic range destroys per-tensor scaling. That is the LLM.int8() result, and reproducing it yourself is far more instructive than reading it. Measure per-channel scaling as the fix and quantify the recovery.

Benchmarks and Metrics

Correctness Tests

1. Every variant matches a reference matmul. For fp32, relative error < 1e-5; state and justify the tolerance, since reassociation changes results legitimately.
2. The simulator matches a reference matmul exactly in integer arithmetic. A cycle-accurate simulator that computes the wrong answer is worthless.
3. Edge shapes: non-square, non-power-of-two, dimensions smaller than the block or array size, and 1×N and N×1.
4. Quantization round-trip error within the derived bound.
5. Accumulator overflow: int8 inputs must accumulate in int32; test the boundary deliberately.
6. Numerical stability at large N — sum order matters; compare against a Kahan-summed float64 reference.
7. The simulator's cycle count matches a hand-computed value for a small case. Do a 4×4 array by hand on paper.

Failure Tests

Expected Difficulties

1. Compiler flags dominate your results. The -O2/-O3 table above shows a conclusion reversing with flags. Report the exact flags with every number, and run the whole sweep at both settings.
2. Auto-vectorisation is opaque. Use -Rpass=loop-vectorize (clang) or read the assembly. "I added #pragma omp simd and it got faster" is not an explanation.
3. You will not beat BLAS. Expect 20–30× off. That is the correct outcome and it is the project's thesis; do not spend two weeks trying to close it.
4. Hardware counters vary by platform. On Apple Silicon, perf-equivalents are limited. Derive DRAM traffic analytically where you cannot measure it, and say which you did.
5. The simulator can absorb unlimited time. Cycle-accurate MACs, fill/drain, and utilisation reporting are the scope. A full memory-hierarchy model is the extension.
6. Thermal throttling on a laptop will corrupt long sweeps. Interleave variants rather than running each to completion, and log temperature or at least run order.

Scope Boundaries

In scope: dense matmul on CPU, the block/vectorise/tile progression, roofline analysis, quantization, a cycle-accurate systolic simulator, a GPU implementation if hardware allows.

Out of scope: actual FPGA synthesis (extension); a full accelerator compiler; sparse matmul; convolution (matmul is the workload); distributed or multi-GPU; a real ML framework integration beyond P13's kernels; power measurement requiring instrumentation you do not have.

Deliverables

1. hwaware/ — C kernels, the sweep harness, the simulator, GPU code if applicable
2. The roofline chart with every variant plotted — the single most legible artifact in Stage 5
3. REPORT.md answering "why can specialised hardware beat a CPU by 100×?" with your own three numbers (58 / 1,679 / simulated)
4. Notebook entries for E2, E7, E9, E14
5. The systolic simulator as a standalone tool, with the TPUv1 derivation as its validation case

Exit Criteria

• CPU progression complete: naive → reordered → blocked → vectorised, all measured at both -O2 and -O3
• E2 complete: block-size optimum found and compared against your cache-size prediction
• Every variant placed on a roofline with measured or derived intensity
• E7 complete: fp32/fp16/int8 throughput and accuracy, with the outlier problem observed and per-channel scaling measured as the fix
• Systolic simulator is cycle-accurate, produces correct results, and reproduces the TPUv1 TOPS figure from k=256 and 700 MHz
• E10/E11 complete: utilisation vs shape and fill/drain overhead measured
• E14 complete: CPU vs BLAS vs simulated accelerator on one chart
• Percentages of peak quoted against dense figures, with the sparsity asterisk noted
• REPORT.md written with a falsified prediction

Extension Ideas

• Memory-hierarchy simulation: add HBM bandwidth and on-chip SRAM capacity limits to the array simulator, then write a tiler that schedules for it. This turns the simulator into a design-space exploration tool and is the strongest extension here.
• FPGA implementation of a small systolic array. Real, and a genuinely distinctive portfolio piece.
• Triton or CUDA kernel for fused attention, measured against PyTorch's.
• Energy modelling: per-operation energy from published figures, comparing architectures on FLOP/joule rather than FLOP/s — which is the metric that actually drives accelerator design.

Connections

Backward: P13-II supplies the kernels and the fusion framing. P12 supplies cache and TLB understanding from the OS side. P01's cost model supplies the workload shapes.

Forward:

• → P15: "can hardware-aware batching substantially reduce end-to-end embedding latency?" is one of the candidate research questions, and this project is its foundation
• → Retroactively: re-read P02's constant-factor result and P13's dispatch crossover. All three are the same phenomenon at different scales, and saying so explicitly in the P14 report is the synthesis that Stage 5 is for

References

• Jouppi, N. P. et al. In-Datacenter Performance Analysis of a Tensor Processing Unit. ISCA 2017.
• Williams, S., Waterman, A., Patterson, D. Roofline: An Insightful Visual Performance Model for Multicore Architectures. CACM 52(4), 2009.
• Goto, K., van de Geijn, R. A. Anatomy of High-Performance Matrix Multiplication. ACM TOMS 34(3), 2008.
• Chen, Y.-H., Emer, J., Sze, V. Eyeriss: A Spatial Architecture for Energy-Efficient Dataflow for Convolutional Neural Networks. ISCA 2016.
• Kung, H. T., Leiserson, C. E. Systolic Arrays for VLSI. Sparse Matrix Proceedings, 1978. The original.
• Drepper, U. What Every Programmer Should Know About Memory. Red Hat, 2007.
• Micikevicius, P. et al. Mixed Precision Training. ICLR 2018.
• Dettmers, T., Lewis, M., Belkada, Y., Zettlemoyer, L. LLM.int8(): 8-bit Matrix Multiplication for Transformers at Scale. NeurIPS 2022.
• Sze, V., Chen, Y.-H., Yang, T.-J., Emer, J. S. Efficient Processing of Deep Neural Networks. Morgan & Claypool, 2020.
• Hennessy, J. L., Patterson, D. A. Computer Architecture: A Quantitative Approach, 6th ed. Morgan Kaufmann, 2017. Chapter 7 on domain-specific architectures.
• Hennessy, J. L., Patterson, D. A. A New Golden Age for Computer Architecture. CACM 62(2), 2019. The framing argument for this whole project.