Systolic Arrays and Dataflow Scheduling

A matrix multiply performed on a general-purpose core is bound by data movement. Each output element is a dot product, and a straightforward implementation reads both operands from memory for every multiply-add. For a multiply of two $n \times n$ matrices that is $2n^3$ operand fetches for $n^3$ multiplies, so the arithmetic-to-memory ratio is two fetches per operation — far below what any memory hierarchy can sustain. The fix that matters is not a faster multiplier. It is a structure in which each operand, once fetched, is used many times before it leaves the array. A systolic array is the most direct expression of that idea: a grid of cells, each holding one value and passing it to a neighbour on the next cycle.

One cell, and what it costs

The repeating element is a multiply-accumulate cell. It holds an accumulator, receives a value from its left neighbour and a value from above, computes the product, adds it to the accumulator, and forwards both inputs to its right and downward neighbours. In the simplest form the cell is one multiplier, one adder, one accumulator register, and four wires.

The operand that arrives from the left travels along a row; the operand from above travels down a column. Each cell therefore sees a sequence of pairs, and the pairing is determined by when each value arrives. The array as a whole performs a matrix multiply if the timing is arranged so that cell $(i,j)$ receives $A[i,k]$ at the same cycle it receives $B[k,j]$, for every $k$. Getting that alignment right is the entire art of using the structure.

The consequence for wiring is significant. Every cell needs connections to its neighbours, and the wires run the length of the array. A $256 \times 256$ array has 65,536 cells and 512 boundary connections, and the internal wire length grows with the array dimension. Clock distribution across a large grid is the practical limit on how large an array can be built before it must be partitioned into tiles with registered boundaries.

Reuse, and the ratio that follows

The reuse argument is easy to state precisely. In an $n \times n$ multiply on an $n \times n$ array, each element of $A$ enters the array once and is used $n$ times as it travels across a row. Each element of $B$ enters once and is used $n$ times as it travels down a column. Total operand traffic into the array is $2n^2$ values for $n^3$ multiply-accumulates, so the arithmetic-to-input ratio is $n/2$. A $128 \times 128$ array therefore performs 64 multiply-accumulates per value that crosses its boundary.

That ratio is the design’s whole value, and it is also its constraint. The array is efficient only when the matrices are large relative to the array dimension. A multiply of two $32 \times 32$ matrices on a $128 \times 128$ array leaves three quarters of the cells idle for most of the computation, and the boundary traffic is a much larger fraction of the total work. The structure is sized for a problem shape, and problems of other shapes either pad, which wastes cells, or are split into tiles, which requires accumulating partial results across passes.

#include <stdint.h>

/* One cell of a weight-stationary array: the weight stays, activations flow. */
struct cell {
    int32_t weight;      /* loaded once during configuration */
};

/* One reduction step. The activation arrives from the left; the partial sum
   arrives from the row above and leaves below with this product folded in. */
static inline int32_t step(const struct cell *c, int32_t activation, int32_t partial_in)
{
    return partial_in + activation * c->weight;
}

The variant above is weight-stationary: each cell holds one element of one operand permanently, and the other operand streams past. This is the arrangement most often used for inference, because weights can be loaded once at configuration and reused across an unbounded number of inputs. The alternative, output-stationary, keeps the accumulator in place and streams both operands, which requires both to be delivered in a skewed order and yields a different balance of boundary traffic against internal movement.

The three dataflows and what each one moves

The choice of which operand to hold stationary determines what has to move and how often.

Weight-stationary. Weights are loaded into the array once and remain. Activations stream in from one edge and partial sums stream out of another. The array can process a stream of independent inputs without reloading, which makes it the natural fit for a fixed model applied to many samples. Its cost is that the accumulator cannot stay in the cell across a batch boundary: partial sums must be written out and read back when a new input starts, unless the array is deep enough to complete the reduction internally.

Output-stationary. Each cell owns one output element and accumulates it in place for the whole reduction. Both operands stream through, and the partial sums never leave the array until the reduction completes. This minimizes accumulator traffic, which is often the largest single term, because partial sums are wide and written repeatedly. The cost is that both operands must be delivered in a skewed schedule, so the buffers that feed the array must be able to present different rows at different times.

Row-stationary and hybrid schedules. The intermediate arrangement keeps a row of one operand and a column of the other in local storage and iterates over the reduction dimension in an order that maximizes reuse at every level of the hierarchy. This is the formulation that treats the array as one level in a memory hierarchy rather than as a standalone structure, and it is what makes the design scale to arrays whose operands do not fit in registers.

The three differ in which term of the traffic equation they minimize, and none minimizes all of them. The trade is fundamental: an operand held in place is not being moved, and every operand that moves occupies a wire and a buffer slot.

Skew, and why the schedule is not optional

Consider a $2 \times 2$ array computing $C = A \cdot B$. Cell $(0,0)$ needs $A[0,0]$ and $B[0,0]$ at the same time. Cell $(0,1)$ needs $A[0,1]$ and $B[1,0]$. If both operands enter the array at cycle zero, the values arrive at the wrong cells. The fix is to delay the input to each row and each column by a number of cycles proportional to its index, so that the wavefronts meet where they should.

For the $2 \times 2$ case, $A$‘s row 0 enters at cycle 0 and row 1 at cycle 1; $B$‘s column 0 enters at cycle 0 and column 1 at cycle 1. The result for $C[0,0]$ emerges from cell $(0,0)$ after the first product, and the diagonal of results emerges one per cycle. This skewing is what gives the array its name: data pulses through the grid like a wave, and the latency from the first input to the last output is proportional to the array dimension plus the reduction length.

The skew is a property of the schedule, not of the hardware. The hardware provides the wires and the cells; the order in which values are presented, the delays applied at the edges, and the accumulation of partial results across tiles are all decided by whatever generates the stream. This is why a systolic array is inseparable from its compiler: the array is a datapath, and the schedule is a program.

Where the approach stops paying

The structure has four limits, and each one shows up in a different regime.

Small problems. Below the array dimension, utilization falls off as the square of the ratio of problem size to array size. A workload dominated by small matrices cannot keep the array busy, and a general-purpose core with a vector unit handles those cases better.

Non-uniform shapes. The array is efficient for dense, rectangular work. A sparse matrix has zeros that the array will happily multiply, and the fraction of useful work is the density. Structured sparsity — a fixed pattern of non-zeros per row — can be exploited by pruning the datapath, but unstructured sparsity cannot, and the array’s efficiency drops to the density.

Data-dependent control flow. A cell has no program counter and no branch. Work that requires decisions per element — branching, gathering, sorting — does not map onto a grid whose every cell does the same thing every cycle. This is the same constraint that limits any wide, statically scheduled datapath, and it is more severe here because a systolic array has no per-lane predication at all: every cell performs the same operation in the same cycle, or the schedule is wrong.

Utilization under a fixed workload. An array sized for the largest model in a family is idle for the smaller ones. Since the array’s area and static power are paid regardless, a system that runs mixed workloads must either partition the array into independently schedulable tiles or accept an average utilization well below peak. Tiling is the standard answer, and it converts the array into a set of smaller arrays whose schedules are coordinated, which brings back the boundary traffic the array was meant to eliminate — one level up.

The trade the structure makes is therefore precise. It exchanges a large amount of wiring and a rigid schedule for an arithmetic-to-memory ratio that scales with the array dimension. That is a good trade for dense, regular, large problems, and a poor one for anything else. The alternatives that keep more flexibility while retaining operand reuse are the topic of Shared Virtual Addressing and Peer-to-Peer Transfers in the memory direction and of Scheduling and Preemption on Accelerator Queues in the control direction.