Systolic arrays

You know a neural network is mostly GEMM. A systolic array is the most direct hardware for it: an $R \times C$ mesh of MAC PEs with registered nearest-neighbor links, fed from on-chip buffers at two edges. Kung’s 1982 argument still frames the design: a GEMM is compute-bound ($O(n^3)$ operations on $O(n^2)$ data), so a structure that uses each fetched value many times can raise throughput without raising memory bandwidth.²

The details that decide whether a given array is any good are:

• Dataflow: which operand is stationary in the PE (weights, outputs, inputs, or Eyeriss’s row-stationary mix), which sets what each level of memory is read for.⁴
• Mapping: how a GEMM’s $M$, $K$ and $N$, or a convolution lowered to a GEMM, are assigned to rows, columns and time, and how many folds that takes.⁶
• Utilization: fill and drain (about $2R + C - 2$ cycles per fold), fragmentation, weight-load stalls and memory-bound layers. TPU v1 averaged 23% of peak MACs doing useful work across six production workloads.⁷

This chapter puts numbers on each, using published TPU, Tesla FSD, Eyeriss and Gemmini data, and derives the runtime, utilization, traffic and energy models in “Under the hood”.

The PE and its timing

In the weight-stationary form used by TPU matrix units, PE($k, j$) holds $w = B[k][j]$ and on every cycle computes

with $a_{\mathrm{out}}$ registered to the right neighbor and $p_{\mathrm{out}}$ to the one below. Activation $A[m][k]$ enters row $k$ at cycle $m + k$ (the skew) and reaches column $j$ at $m + k + j$. The partial sum for output $(m, j)$ is at row $k$ at the same cycle $m + k + j$, so the operands always coincide, and $C[m][j]$ leaves the bottom of column $j$ at cycle $m + (R - 1) + j$. The outputs leave skewed and are de-skewed into accumulators.⁶⁷

Properties that follow, and that Kung listed as the criteria for a good systolic design:

• Each input is used many times ($C$ times per activation, $M$ times per weight).
• Extensive concurrency: all $R \cdot C$ PEs can work at once in steady state.
• Only a few kinds of simple cells.
• Simple, regular data and control flow: every wire is nearest-neighbor, so wire length doesn’t grow with array size; the only global signal is the clock.²

Kung contrasted these “pure” systolic arrays with semi-systolic ones that broadcast an input to all cells or sum results through a fan-in adder tree. Broadcast and fan-in also achieve reuse, but their wires grow with the array and are hard to scale without slowing the clock.² The original Kung–Leiserson matrix multiplier was a hexagonally connected array in which all three matrices moved; today’s accelerators use the rectangular, one-operand-stationary form.¹

TPU v1 as the reference implementation

TPU v1’s matrix unit is 256 × 256 8-bit MACs. Products (16-bit) go to 4 MiB of 32-bit accumulators below the unit, organized as 4096 vectors of 256; the unit produces one 256-element partial sum per cycle. Activations come from a 24 MiB on-chip Unified Buffer; weights come from off-chip DRAM through a Weight FIFO four tiles deep. The unit holds one 64 KiB weight tile plus a second for double buffering, “to hide the 256 cycles it takes to shift a tile in”.⁷ Mixed 8/16-bit operands run at half speed and 16/16 at a quarter.⁷

Definitions, precisely

• WS: each weight stays in a PE’s register file to maximize convolutional and filter reuse; the PE runs every MAC that uses that weight before it is replaced. Partial sums are reduced across PEs and, if not finished, spilled to the global buffer, so buffer capacity can limit how many filters are in flight.⁴
• OS: each output accumulates in one PE. Eyeriss distinguishes variants by which outputs are computed together: one output channel with many pixels (for convolution layers), many channels with many pixels, or many channels with one pixel (for fully connected layers).⁴
• IS: the mirror of WS for the activation operand. SCALE-Sim models OS, WS and IS on a systolic array and finds WS beats IS when a layer has more output pixels than weights, and the reverse otherwise, because loading the stationary matrix costs cycles with no compute; the less often it must be reloaded, the better.⁵
• NLR: no register file in the PE; all area goes to a larger global buffer. It has the lowest DRAM energy but pays buffer-level energy for nearly every access.³
• RS: each PE performs a 1-D convolution of a filter row with an input row, keeping the row of weights, a window of activations and one partial sum in its 0.5 kB register file; rows of PEs then sum their outputs vertically for the 2-D convolution.⁴

Under the same area and 256 PEs, Eyeriss’s analysis found WS lowest for weight energy and OS lowest for partial-sum energy, but RS lowest in total, because it optimizes all three data types. In fully connected layers there is little reuse and every dataflow spends a significant amount of its energy reading weights.³

Two observations from the original systolic work still guide the choice. When there are more weights than cells, a design where partial results stay generally needs less I/O than one where they move. And moving partial results needs wider paths than moving weights, because results carry more bits for numerical accuracy.² In modern terms: WS ships 32-bit sums of 8-bit products between PEs; OS ships only the 8-bit operands.

SCALE-Sim’s authors conclude that while hyper-parameters and array size affect which dataflow is best, fixing one usually costs little, so the cost of supporting several should be weighed carefully.⁵ Gemmini, an open-source generator, makes it a configuration option: its PEs run either WS or OS.¹¹ At the extreme, analog crossbars that store weights as conductances are “the ultimate form of a weight stationary dataflow”, covered in Sparsity, in-memory and analog compute.³

The Expert view adds the array size $n$ (2–8), the reduction depth $K$ (1–24) and up to 24 rows of $A$; $B$ is $K \times n$. The cycle count follows the SCALE-Sim model exactly: $\text{folds} \times (3n + T - 2)$, with $T = M$ for WS (folds over $K$) and $T = K$ for OS (folds over $M$), and no overlap between folds.⁶ Try:

• $M = 16$, $K = 4$ versus $M = 4$, $K = 16$ in each dataflow. Which dimension does each one amortize fill and drain over?
• Read the bar chart: where are the load, fill and drain ramps, and how wide is the flat top?
• Compare “buffer reads” with “minimum reads”. When does each dataflow re-read an operand, and which one?
• At $n = 8$, how large must $M$ be before WS utilization passes 50%?

Buffer reads count every operand entering the array edge; results written count every vector leaving it, including WS partial sums sent to the accumulators between $K$-folds. Memory never stalls in this model.

Mapping GEMM dimensions

SCALE-Sim describes any dataflow as a projection of the GEMM $C(M \times N) = A(M \times K) \cdot B(K \times N)$ onto two spatial dimensions $S_R$, $S_C$ and one temporal dimension $T$. For OS, $M$ maps to $S_R$, $N$ to $S_C$ and $K$ to $T$: each PE accumulates over $K$ in time. WS and IS keep one operand in space and stream the other, so the dimension that maps to time is the one not held by the stationary operand.⁶ In this chapter’s sim, WS holds $B$ ($K$ on rows, $N$ on columns) and streams the $M$ rows of $A$; OS holds $C$ and streams $K$.

On an $R \times C$ array the mapping needs $F = \lceil S_R/R \rceil \cdot \lceil S_C/C \rceil$ folds.⁶ Which dimension is folded decides which data is re-read: a WS array folding over $K$ sends partial sums to the accumulators after every fold and reads the activations once per column tile; an OS array folding over $M$ re-reads all of $B$ for every row tile.

Convolution dimensions

For a layer with batch $N$, $C$ input channels, $K$ filters of $R \times S$, and $P \times Q$ outputs, lowering gives a filter matrix $K \times CRS$ and a data matrix $CRS \times NPQ$, with each input duplicated up to $RS$ times.¹² In GEMM terms the reduction depth is $C \cdot R \cdot S$, the output columns are the $K$ filters and the streamed dimension is $N \cdot P \cdot Q$. Three consequences:

• Early layers tend to have few channels but large filters; late layers have many channels but small filters. The product $CRS$, the reduction depth, is usually fairly large either way, which is why lowering performs consistently.¹²
• Layers with fewer output channels than the array has columns leave columns idle. TPU v1 measured that on active cycles of one CNN only about half the 65,536 MACs held useful weights “because some layers in CNN1 have shallow feature depths”.⁷
• Depthwise convolution (one filter per channel) has a reduction of only $R \cdot S$, typically 9, so reuse is low; Gemmini reports MobileNetV2 maps poorly to spatial arrays for exactly this reason.¹¹

Accumulator sizing and precision

TPU v1 sized its accumulators from the roofline: about 1350 operations per weight byte are needed to reach peak, rounded up to 2048 accumulator vectors, then doubled to 4096 so the compiler could double-buffer.⁷ TPU v2 and v3 use a 128 × 128 array with “streaming LHS and results” and a “stationary RHS”, bfloat16 multiplies and float32 accumulation.⁸ Tesla’s FSD accelerator accumulates 8-bit × 8-bit products into 30-bit accumulators.¹⁰ The rule behind all three: the accumulator must be wide enough for $K$ products without overflow or excessive rounding, and that width is what makes moving partial sums expensive. The number formats themselves are in Number formats.

The fixed cost per fold

For one fold on an $R \times C$ array SCALE-Sim gives $\tau = 2R + C + T - 2$ for OS, WS and IS alike. In OS, the far PE receives its first operands at $R + C - 2$, accumulates for $T$ cycles and the results drain in $R$. In WS, loading the weights takes $R$ cycles (no skew needed, since nothing computes during the load), the first input reaches the last column after $C - 1$, streaming takes $T$ and the reduction down the column takes $R - 1$.⁶ The overhead $2R + C - 2$ is independent of $T$, so for a full fold

(Computed from the formula; it assumes no overlap between folds, which double-buffered designs such as TPU v1 partly achieve.)

Fragmentation

With $F = \lceil S_R/R \rceil \lceil S_C/C \rceil$ folds, the useful share of PE slots is $S_R S_C / (F R C)$. For TPU v1’s 600 × 600 example: $360{,}000 / (9 \times 65{,}536) \approx 61\%$ on 256 × 256 versus $360{,}000 / (4 \times 262{,}144) \approx 34\%$ on 512 × 512, consistent with the reported 18 µs versus 32 µs.⁷ Fragmentation is worse in two dimensions than one, which the authors compared to internal fragmentation of large memory pages.⁷

Measured: TPU v1 performance counters

Table 3 of the TPU paper splits the matrix unit’s time across six production workloads (MLP0 and 1, LSTM0 and 1, CNN0 and 1). Mean values: array active 28% of cycles, useful MACs 23% of peak, unused MACs on active cycles 5%, weight stalls 43% and weight shifting 12%. CNN0 reaches 78.2% useful MACs and 86 of 92 peak TOPS; CNN1 reaches only 22.5%, half its active MACs idle from shallow layers and 28.1% of cycles stalled on weights for its fully connected layers.⁷ Four of the six workloads sit on the bandwidth-limited slope of the TPU’s roofline, whose ridge is at about 1350 operations per weight byte.⁷

Mitigations, and what they cost

• Overlap folds. Double-buffer weights so the next tile shifts in during compute (TPU v1). Costs a second weight register per PE.⁷
• Smaller, more numerous arrays. 128 × 128 instead of 256 × 256 (TPU v2), and more of them per core (two in TPU v3, four in TPU v4). Costs reuse per operand and so buffer bandwidth.⁸⁹
• Bigger batches raise $T$ for WS, limited by latency targets and memory capacity.⁷
• Pick the dataflow by shape: WS amortizes overhead over $M$, OS over $K$. The sim makes the mirror symmetry visible.⁵
• Scale out. SCALE-Sim’s follow-up asks whether one big array or several partitions is faster for a given layer, and builds the runtime model above to search that space.⁶

Points worth noting:

• Useful MACs track active cycles closely for MLP0, MLP1, LSTM0 and CNN0, so there the loss is time (stalls), not space (fragmentation). CNN1 loses about half its active MACs to unused slots, and LSTM1 about 40%.
• Weight shifting is 13–17% of cycles in the MLPs and LSTMs even with double buffering: when each weight tile is used for few rows, its 256-cycle shift can’t be hidden behind compute.
• Eyeriss reports its row-stationary dataflow 1.4–2.5× more energy-efficient than WS, OS and NLR in AlexNet’s convolution layers and at least 1.3× in fully connected layers at batch sizes of 16 or more.⁴
• The TPU v4 authors attribute part of their energy advantage to reuse: a 128 × 128 MXU reuses each 128-entry input 128 times, so the on-chip SRAM is read less per operation.⁹

Reuse arithmetic

In an $R \times C$ WS array, each activation read is used $C$ times, each weight read $T$ times per load, and each partial sum leaving the array represents $R$ MACs. Operations per operand read therefore scale with the array edge, which is why quartering a 256 × 256 array into four 128 × 128 arrays halves operations per operand (4× to 2× on Google’s chart, normalized to 64 × 64) while improving utilization 1.6×.⁸ The buffer must supply the difference: four arrays of edge 128 need twice the operand bandwidth of one array of edge 256 at the same FLOPS.

Wide sums on the wires

Kung noted in 1982 that paths carrying partial results must be wider than paths carrying weights, because results carry more bits for accuracy.² With int8 operands and 32-bit accumulation, a WS column carries 32-bit sums PE to PE, four times the width of the activations crossing the row. Tesla’s FSD accelerator takes the other side: it computes 96 output pixels by 96 output channels at a time with the input-channel and kernel loops innermost, accumulating 8-bit × 8-bit products into 30-bit sums held in the MAC array, and lists “in place data reuse vs result movement” among its power choices.¹⁰

Clocking a big array

Nearest-neighbor wiring removes long data wires but not the clock. Kung warned that large 2-D systolic arrays may need a slower global clock to absorb clock skew, while one-dimensional arrays tolerate large skew between their ends.² Today that is the job of clock tree synthesis, but the clock remains the one signal that spans the whole array.

Flexibility tax

Supporting several dataflows needs muxes and extra registers in every PE, and SCALE-Sim found fixing one dataflow rarely costs much; Gemmini makes it a build-time choice so designers can measure.⁵¹¹ For layers whose matrices are small, the shape mismatch matters more than the dataflow: a tiny $M$ on a WS array is mostly fill and drain. More on batch and decode in The memory wall.

The models behind the simulation and the numbers above, step by step, with their assumptions.

1. Dataflow as a loop order

A GEMM is three nested loops. A dataflow chooses which index is spatial (unrolled across PEs) and which is temporal, and the stationary operand is the one whose indices don’t change in the innermost temporal loop. Illustrative pseudo-code for one tile:

# C[M][N] += A[M][K] * B[K][N] on an R x C array
# Weight-stationary: B[k][n] pinned in PE(k, n)
for k0 in 0..K step R:            # fold: load B[k0:k0+R][0:C]
  for m in 0..M:                  # time: stream rows of A
    parallel for k, n in R x C:   # space
      C[m][n] += A[m][k0+k] * B[k0+k][n]

# Output-stationary: C[m][n] pinned in PE(m, n)
for m0 in 0..M step R:            # fold: new tile of outputs
  for k in 0..K:                  # time: stream A columns and B rows
    parallel for m, n in R x C:   # space
      C[m0+m][n] += A[m0+m][k] * B[k][n]
  drain C[m0:m0+R][0:C]

Tesla’s slides show the same idea for convolution: the seven-deep loop nest is reordered so output pixels and output channels step by 96 in the outer loops and input channels and kernel positions run innermost.¹⁰

2. Cycle-exact timing

Weight-stationary, $R = C = n$, one fold, cycles counted from the start of the fold:

1. Load: weight rows shift in from the top, bottom row first, one row per cycle: $n$ cycles. Nothing computes, so no skew is needed.⁶
2. Stream: $A[m][k]$ enters row $k$ at stream cycle $s = m + k$ and meets PE($k, j$) at $s = m + k + j$.
3. Exit: $C[m][j]$ leaves the bottom at $s = m + (n - 1) + j$. The last one, $m = M - 1$, $j = n - 1$, leaves at $s = M + 2n - 3$.
4. Total: $n + (M + 2n - 2) = 3n + M - 2$, which is SCALE-Sim’s $2R + C + T - 2$ with $T = M$.⁶

Output-stationary: $A[m_0 + r][k]$ enters row $r$ at cycle $k + r$, $B[k][j]$ enters column $j$ at $k + j$, so PE($r, j$) does its $k$-th MAC at cycle $k + r + j$. The last MAC is at $K + 2n - 3$; then $n$ cycles shift the results out the bottom: $K + 3n - 2 = 2R + C + T - 2$ with $T = K$. The busy count at stream cycle $s$ is the number of $(r, j)$ with $0 \le s - r - j < T$, a diagonal band; it reaches $n^2$ only if $T \ge 2n - 1$. That is why the sim’s 4 × 4 default ($T = 4$) peaks at 12 busy PEs.

3. Runtime and utilization with folds

⁶

Assumptions: operands always available (no memory stalls), no overlap between folds, outputs leave without stalling, and every fold pays the full overhead even when partly empty. Real designs overlap weight loading with compute (TPU v1 double-buffers) and may pipeline consecutive folds, so $\tau$ is a pessimistic bound on overhead and an optimistic one on memory.⁵⁷ SCALE-Sim’s authors checked its cycle counts against an in-house RTL model for one case: output-stationary multiplication of matrices the size of the array.⁵

4. Buffer traffic

Counting each operand that enters the array edge as one buffer read (the sim’s model, $N = C$):

• WS: reads = $KN$ (each weight once) + $MK$ (each activation once per column tile); accumulator writes = $MN \lceil K/R \rceil$, of which all but the last fold also need a read-modify-write.
• OS: reads = $MK$ (each activation once) + $KN \lceil M/R \rceil$ ($B$ re-read for every row tile of outputs); result writes = $MN$.
• Lower bound for either: $MK + KN$, every operand read once. Without reuse a MAC needs two operand reads, so $2MKN$.³

Example from the sim: $n = 4$, $M = 16$, $K = 4$. WS reads $16 + 64 = 80$ for 256 MACs (3.2 MACs per read, the minimum possible); OS needs 4 row tiles and reads $64 + 64 = 128$ (2.0 MACs per read). Swap to $M = 4$, $K = 16$ and the roles reverse in cycles, though here both read 128 because $A$ and $B$ are both 64 values.

5. Energy from reuse

Eyeriss models a value’s data-movement energy from how many times it is reused at each level. If a value is read from DRAM $a$ times, each of those copies from the global buffer $b$ times, passed across the array $c$ times and read from the register file $d$ times, then⁴

with $\mathrm{EC} = 200$, 6, 2 and 1 MAC-equivalents. Partial sums are counted the same way for their reads and writes.⁴ A worked example with those costs (illustrative): an activation that must meet 128 weights. Read from the buffer for every MAC it costs $128 \times 6 = 768$. Read once from the buffer and passed PE to PE across a 128-wide WS row it costs $6 + 128 \times 2 = 262$, about 3× less. Read from DRAM every time it would cost $128 \times 200 = 25{,}600$. Sze et al. estimate AlexNet’s 724M MACs would need nearly 3000M DRAM accesses with no reuse and 61M with perfect on-chip reuse.³

6. Lowering convolution

Filter tensor $F$ ($K$ filters × $C$ channels × $R \times S$) reshapes to $F_m$ of $K \times CRS$; the data matrix $D_m$ is $CRS \times NPQ$, built by copying each $R \times S \times C$ input patch into a column; $O_m = F_m \cdot D_m$ is $K \times NPQ$.¹² Each input value appears up to $R \cdot S$ times in $D_m$, so materializing it multiplies activation memory and traffic by up to $R \cdot S$; cuDNN instead forms tiles of $D_m$ lazily in on-chip memory, and Gemmini does the equivalent in hardware at the array edge.¹²¹¹ Sze et al. note that the same redundancy, or a complex access pattern, is the price of treating convolution as plain matrix multiplication, and that its row-stationary dataflow exploits the convolutional reuse directly instead.³

7. What the simulation leaves out

• Memory: buffers deliver every operand on time; there are no DRAM, bandwidth or bank-conflict stalls.
• Overlap: folds run back to back with the full load, fill and drain each time (no double buffering).
• Shape: $B$ always has as many columns as the array ($N = n$), so column fragmentation isn’t modeled.
• Numerics: small positive integers with unlimited precision; no 8-bit saturation or accumulator width.
• Energy and area: reads are counted, not weighted by the cost of each memory level.