Systolic Arrays

A systolic array is a 2D grid of MAC units where operands are pumped between nearest neighbors, so each value fetched from memory is reused across an entire row or column of PEs.

Mechanism

An $R\times C$ grid of processing elements, each holding a register and doing one MAC per cycle. Operands enter at the edges and shift one hop per cycle; no PE addresses memory directly.

• Peak = $2\cdot R\cdot C$ FLOP/cycle. A 256×256 array at 1 GHz = 131 TFLOP/s from one block.
• Memory ports required: $O(R+C)$ per cycle, not $O(RC)$. That $O(n)$ vs $O(n^2)$ gap is the entire architectural argument: bandwidth grows linearly while compute grows quadratically with array size.
• Fill/drain latency $\approx R+C$ cycles; a 256×256 array needs ~512 cycles of pipeline before the first result and after the last.

Dataflow variants (what stays put)

Variant	Stationary in PE	Streams through	Best when
Weight-stationary	one weight	activations in, partial sums accumulate down	weights reused across many activations (large batch/large $M$)
Output-stationary	one accumulator	weights and activations both stream	long reduction dimension $K$; avoids partial-sum movement
Row-stationary	a row of weights + reuse of activation rows	sliding-window overlaps	convolutions, where the same activation feeds many outputs

Choosing between them is choosing which operand’s movement energy you eliminate. Off-chip DRAM access costs on the order of 100–1000x the energy of a MAC, and even a local register read is ~1–2x a MAC, so the dataflow choice is primarily an energy argument, then a bandwidth one.

Weight-stationary is the common choice for large matmul engines; output-stationary appears where $K$ is large and precision of accumulation matters.

The utilization cliff

Utilization for a GEMM $M\times K\times N$ mapped to an $R\times C$ array is roughly

$U=\frac{M}{\lceil M/R\rceil R}\cdot \frac{N}{\lceil N/C\rceil C}\cdot \frac{K}{K+(R+C)}$

three multiplicative penalties:

1. Quantization of M and N to the array dimensions. $M=129$ on a 128-row array uses 2 passes at 50%, a step function, not a gradient.
2. Fill/drain amortization: small $K$ never amortizes the $R+C$ pipeline depth. $K=64$ on a 256×256 array is under 15% utilization before anything else.
3. Batch-1 decode is the worst case. A GEMV has $M=1$: one row of a 256-row array is active → 0.4% utilization. This is the single most important fact about systolic arrays for LLM inference, and it is the same phenomenon as low tensor-core efficiency at small M (Tensor Cores).

Compare to a GPU’s SIMT array, which degrades more gracefully because its granularity is a 32-lane warp, not a 128–256-wide array edge.

Mitigations

• Smaller or partitionable arrays: many independent 32×32 or 64×64 tiles instead of one 256×256, trading peak density for shape flexibility.
• Batching / continuous batching to raise $M$ above the array width.
• Speculative decoding: turns $M=1$ into $M=k$ verification, directly buying utilization.
• Fold other dimensions into M: multiple sequences, multiple heads, beam candidates.
• Split-K to keep the array fed when $K$ is the only large dimension, at the cost of a reduction.
• Padding is correct but wasteful: measure achieved MACs, not padded ones.

Tip

Always ask for utilization at your actual shapes, not peak TOPS. Array-width quantization plus fill/drain means a spec-sheet-fast part can deliver single-digit percent on decode-shaped GEMVs.

Dataflow Architectures · Tensor Cores · GEMM Tiling · Performance Modeling · Accelerator Compilers