The Roofline Model on a Laptop: Measuring the Bandwidth Slope, the FMA Ceiling and the Ridge

Before optimising a loop it helps to know which wall it is pressed against: the rate at which the core can do arithmetic, or the rate at which memory can feed it. The roofline model of Williams, Waterman and Patterson answers that with two machine numbers and one number per kernel. This post, the first in a short series on single-node performance, measures both machine numbers on one core of an Apple-silicon laptop (sysctl -n machdep.cpu.brand_string prints Apple M1, hw.ncpu prints 8) and then places real kernels under the result. Everything below was compiled with Apple clang 16.0.0 (clang --version: Apple clang version 16.0.0 (clang-1600.0.26.6)) at -O3 -mcpu=native, and every number comes from the programs shown.

Arithmetic intensity

The kernel's number is its arithmetic intensity: the floating-point operations it performs per byte it moves between the core and the memory level under study,

Two kernels in double precision make the range concrete. The STREAM triad a[i] = b[i] + s*c[i] loads two 8-byte values, stores one and performs one fused multiply-add: $I=2/24\approx 0.083$. A degree-64 polynomial evaluated by Horner's rule loads one value, stores one and performs 64 multiply-adds in registers: $I=128/16=8$. A cache that allocates on write would also read each destination line before overwriting it, so the hardware may move 32 bytes per triad element rather than 24. STREAM counts 24, and so do I, for the roof and for the kernels alike; what matters is that the convention is applied the same way to both.

The bound

Let $\pi$ be the peak arithmetic rate and $\beta$ the peak bandwidth of the level that feeds the kernel. A kernel with intensity $I$ cannot run faster than either limit allows:

On log-log axes the right-hand side is a line of slope one that meets a horizontal line at the ridge point $I^{\ast }=\pi /\beta$. Left of the ridge the bandwidth roof is the lower bound, so the best a kernel there can do is run at memory speed, and only more bandwidth or more reuse raises that ceiling; right of the ridge the arithmetic roof is the lower one. The position says which roof caps a kernel, not what holds a particular implementation below it. A scalar dot product sits left of the ridge, but if every addition waits for the previous one it can run no faster than a single dependent chain, 1.5 GFLOP/s on this core as measured below, far under its bandwidth roof; splitting the sum into independent chains would leave its intensity unchanged and still raise its speed. The paper also draws ceilings under the flat roof, lower limits that apply when a kernel cannot use some feature of the core, such as SIMD lanes or enough independent operations. Two of them are measured below.

Measuring β with a triad

The bandwidth roof comes from a STREAM-style triad on one core:

__attribute__((noinline)) void triad(double *restrict a, const double *restrict b,
                                     const double *restrict c, double s, size_t n) {
    for (size_t i = 0; i < n; i++) a[i] = b[i] + s * c[i];
}

clang turns the loop body into one fmla.2d per two elements, so the 2-flops-per-element count holds. The harness touches every page before timing, repeats the kernel until each trial has moved about 1.5 GB (or one pass, for the large arrays), discards three warm-up trials, keeps fifteen, and runs the whole sweep three times: each size gets 45 trials. Array sizes run from 4 KiB to 256 MiB per array in half-octave steps.

Three plateaus are visible. About 196 GB/s holds up to 64 KiB per array, which is consistent with the two arrays that are read (128 KiB together) fitting in the 128 KiB L1 data cache of a performance core while the stores drain further out. About 124 GB/s holds while the data lives in the 12 MiB L2, and the curve settles between 8 and 24 MiB of working set. STREAM's run rules ask for each array to be at least four times the last-level caches used in the run, here the 12 MiB L2, so at least 48 MiB, so the DRAM roof comes from the three sizes that qualify:

I take $\beta =46$ GB/s, the mean of the three best values. The spread between best and median is real: this was a laptop in normal use, and DRAM-bound timings were the noisiest of everything I measured.

Measuring π with independent FMA chains

A fused multiply-add depends on its accumulator, so a loop that updates one accumulator runs at the latency of the instruction, not its throughput. The compute roof therefore comes from a kernel that updates $K$ independent accumulators held in registers, with no memory traffic at all:

#define VEC_KERNEL(K)                                                       \
__attribute__((noinline)) double vec_##K(double seed) {                     \
    float64x2_t acc[K];                                                     \
    const float64x2_t x = vdupq_n_f64(1e-5), y = vdupq_n_f64(1e-4);         \
    _Pragma("clang loop unroll(full)")                                      \
    for (int k = 0; k < K; k++) acc[k] = vdupq_n_f64(seed + k);             \
    for (long it = 0; it < ITERS; it++) {                                   \
        _Pragma("clang loop unroll(full)")                                  \
        for (int k = 0; k < K; k++) acc[k] = vfmaq_f64(acc[k], x, y);       \
    }                                                                       \
    float64x2_t s = acc[0];                                                 \
    _Pragma("clang loop unroll(full)")                                      \
    for (int k = 1; k < K; k++) s = vaddq_f64(s, acc[k]);                   \
    return vaddvq_f64(s);                                                   \
}

The inner loop compiles to exactly $K$ fmla.2d instructions plus the loop counter. A scalar twin uses __builtin_fma on double accumulators; that file is built with -fno-vectorize -fno-slp-vectorize added, so the scalar version stays scalar. Each cell is the median of three per-pass medians, each over 21 timings:

Each chain adds about 3 GFLOP/s until sixteen, then nothing. The single chain measures latency; the plateau measures throughput; their ratio, sixteen, is how many multiply-adds must be in flight to keep the units busy, which is Little's law applied to a pipeline. If the FMA latency is four cycles, the numbers imply a clock of about 3.05 GHz and four 128-bit FMAs issued per cycle; I did not measure the clock directly. The best single timing gives $\pi =49.8$ GFLOP/s in double precision. The scalar plateau is half of that: the same number of FMAs per second, one lane instead of two.

TIP

One dependent chain uses a sixteenth of this core's floating-point capacity. Splitting an accumulator into several partial sums is the cheapest optimisation in this series, and the compiler will not do it for floating-point reductions on its own, because reassociating the additions changes the rounding.

With both numbers the ridge sits at $I^{\ast }=49.8/46\approx 1.1$ flop/byte.

Kernels under the roof

To walk a kernel across the ridge I used the polynomial from the figure, with every coefficient equal to c so that the coefficients live in one register, and the degree D fixed at compile time:

for (size_t i = 0; i < n; i++) {                 /* n = 16 Mi doubles: 128 MiB per array */
    double xi = x[i], acc = c;
    for (int k = 0; k < D; k++) acc = acc * xi + c;   /* fully unrolled */
    y[i] = acc;
}

Each element costs 16 bytes and $2D$ flops, so $I=D/8$ and $D=1,2,\dots ,128$ spans 1/8 to 16 flop/byte. A second version, written with NEON intrinsics, evaluates twelve vectors (24 elements) at a time and interleaves their steps, so that twelve independent chains are always available. The dot product (two loads per multiply-add, $I=1/8$) and the triad complete the set. All kernels were timed round-robin within each trial, so slow periods on the machine hit them equally.

GFLOP/s, medians of three passes of eleven timings. The dot product ran at 4.61 GFLOP/s (roof 5.7) and the triad at 3.40 (roof 3.8).

Left of the ridge these kernels reach the slope. Degrees 1 to 4, the dot product and the triad all land at 80 to 92 % of $\beta I$, as close to the bandwidth roof as the noise on this machine allows, so cheaper arithmetic has little left to gain here.

Right of the ridge the compiled loop gets slower as the work grows. The vectoriser's remarks explain it: for $D\le 8$ clang reports vectorized loop (vectorization width: 2, interleaved count: 4), but for $D\ge 16$ it reports interleaved count: 1. Each loop iteration is then one vector carrying two elements through $D$ dependent multiply-adds, and the only other parallelism is whatever the out-of-order core can overlap between consecutive iterations. The longer the chain, the fewer iterations fit in that window: degree 64 runs at about five chains' worth of the single-chain rate (15.4 / 3.05) and degree 128 at about three. The kernel is compute-bound in the roofline sense and latency-bound in practice.

Exposing independent chains lifts it to a plateau of about 30 GFLOP/s, 61 % of π. The remaining gap has a visible cause in the assembly: NEON's fmla overwrites its addend, and in Horner's rule the addend is the constant c, so every step needs a register copy before the multiply-add:

	mov.16b	v31, v1              ; copy c
	fmla.2d	v31, v17, v8         ; c + x * acc
	mov.16b	v8, v1
	fmla.2d	v8, v16, v9

That doubles the instruction count per flop. Without hardware counters I cannot prove that the copies are what caps the plateau, but the plateau sits where a loop with two instructions per FMA would be expected to sit, well above the scalar ceiling and well below π. The same intrinsics kernel is slower than clang's loop at low degree, where it streams memory less efficiently; I tuned it for the compute-bound end only.

WARNING

The roof is a measurement of one access pattern, not a law. Individual best trials of the triad inside this round-robin harness reached 54 GB/s and of the dot product 49 GB/s, above the 46 GB/s roof from the separate sweep, and a one-read-one-write stream like the low-degree polynomials does not have the same bandwidth as the triad's two-reads-one-write. On this machine I would not draw conclusions from differences smaller than about 15 % in the memory-bound region.

What the ridge means for matrix multiplication

A ridge near 1 flop/byte is low: one flop per byte of DRAM traffic is enough for the arithmetic roof, not DRAM, to be a kernel's ceiling. Dense matrix multiplication clears it by a wide margin on paper. $C=AB$ for $N\times N$ matrices performs $2N^3$ flops; if each matrix crossed the memory bus exactly once (read $A$ and $B$, read and write $C$), the traffic would be $32N^2$ bytes and the intensity

which is 64 flop/byte at $N=1024$, sixty times the ridge. Every extra trip a matrix makes over the bus lowers it, so this is the best case, not a typical value. The textbook loop is nowhere near that. Its inner statement s += A[i][k] * B[k][j] loads two doubles per multiply-add, an intensity of 1/8, and if those bytes came from DRAM the roof would be 5.7 GFLOP/s. The whole difference between $1/8$ and $N/16$ is reuse: whether a value brought into cache is used once or many times before it is evicted. Once the data does live in cache, the slope that matters is the cache's, 124 GB/s for L2 and 196 GB/s for L1, not DRAM's.

The next post measures the classic loop orders and tiling against exactly these lines, and the SIMD post looks at the vector instructions that the flat part of the roof assumes.

Further reading

• Samuel Williams, Andrew Waterman and David Patterson, Roofline: An Insightful Visual Performance Model for Multicore Architectures, Communications of the ACM, 2009, and the longer Berkeley technical report.
• John McCalpin, STREAM: Sustainable Memory Bandwidth in High Performance Computers, with the run rules on array sizes.
• Ulrich Drepper, What Every Programmer Should Know About Memory, for the cache and DRAM mechanics behind the bandwidth plateaus.
• The Empirical Roofline Toolkit from Berkeley Lab, which measures roofs on any machine, and NERSC's notes on roofline performance analysis.
• Sergey Slotin, Algorithms for Modern Hardware, especially the chapters on instruction-level parallelism and memory.