Dataflow is a schedule for calculating the values of a matrix multiply

This animation shows one concrete, testable case: output-stationary systolic matrix multiplication, computing C = A × B. Each blue activation and amber weight is injected once at an edge, travels through neighboring cells, and is reused at each one it visits. The green partial sum stays local until its output is complete.

The key insight is that reuse is latent in the algebra but invisible to hardware by default. To compute C[i,j] = Σ_k A[i,k]·B[k,j], a conventional processor fetches each A[i,k] and B[k,j] from memory on demand — even though A[i,k] is needed for every column j of C, and B[k,j] for every row i. Nothing forces a cache to hold it for the next output; without explicit scheduling, N³ multiplications incur N³ memory reads and operand reuse is zero by default. The systolic schedule makes reuse structural: inject A[i,k] once at the left edge and let it travel east through all N cells in row i; inject B[k,j] once at the top and let it flow south through all N cells in column j. One memory read funds N multiplications — exactly the ratio a bandwidth-limited system needs.

Correct the mental image: the grid is not a river carrying a generic stream. It is a spatially scheduled computation. Each datum has a named role, a finite route, a known arrival cycle, and is discarded after its last consumer. "Dataflow" here means deciding which operand stays local and which values move, so one slow memory transfer buys many operations.

Output-stationary systolic array

Blue A[i,k] values enter from the left; amber B[k,j] values enter from the top. Their injection times are staggered so matching k values arrive in a cell on the same cycle. Each PE first multiplies the two arriving operands (M), then accumulates the product into its local register (AC), keeping C[i,j] in place. Tokens appear at the bottom-left (A) and top-right (B) corners of the PE they are currently visiting.

Each token has a finite route: one MAC per cell it visits, then it exits the array.

Glossary

PE: Processing Element — one cell in the array; holds a local accumulator register and performs one MAC per clock cycle.
MAC: Multiply-Accumulate — two sub-steps: M multiplies the two arriving operands; AC adds (accumulates) the product into the local register. Together: acc += a × b.
DMA: Direct Memory Access — hardware engine that moves tensor tiles from DRAM to the array's edge buffers without CPU involvement.
TPU: Tensor Processing Unit — Google's application-specific chip for neural-network workloads; its matrix multiply unit is a systolic array.
VLSI: Very Large Scale Integration — the circuit-density era that motivated Kung's systolic work: arithmetic gates became cheap while off-chip I/O remained expensive.

The exact computation

The N×N grid calculates all N² entries of C, each a dot product of length N. For the 3×3 base case:

C[i,j] = A[i,0]·B[0,j] + A[i,1]·B[1,j] + A[i,2]·B[2,j]

Cell (i,j) owns only C[i,j]. On cycle t, it receives A[i,k] and B[k,j] when t = i + j + k. Both operands arrive simultaneously, so the multiply and accumulate happen in one clock. The staggered injection times are what guarantee this alignment.

The original problem

Kung's systolic-architecture argument was driven by VLSI economics: arithmetic became cheap and plentiful faster than off-chip I/O and long, global interconnect. The goal was to obtain many computations per expensive memory access using simple, regular, locally connected processing elements.

Modern ML accelerators inherited that idea because dense linear algebra has exactly the regular reuse pattern the scheme needs. A TPU matrix unit is not a general replacement for memory; it is a specialized device for turning a tiled tensor transfer into a large amount of nearby MAC work.

Naïve	Systolic
Fetch both operands from memory for every MAC.	Inject at edge, forward to neighbors, accumulate locally.
3×3: 54 reads for 27 MACs — zero reuse.	18 injections; each value used by 3 cells — 3× reuse.
Arbitrary reads, global arbitration.	Fixed routes, fixed timing, short wires.

What stays vs. moves

"Dataflow" does not mean all data moves. Three common schedules differ in which operand stays local:

Output-stationary — partial sums stay; inputs and weights move. This animation.
Weight-stationary — weights stay; activations and partial sums move.
Input-stationary — activations stay; weights and partial sums move.

The best choice depends on tensor shape, buffer size, energy cost per movement, and which operand benefits most from reuse.

Array size ≠ matrix size. The physical PE grid does not need to match the dimensions of the matrices being multiplied. A 256×256 array can compute a 16384×16384 matrix multiply by tiling: break each matrix into 256×256 blocks, stream matching pairs of tiles through the array one pair at a time, and accumulate partial results across tiles. The array sees an infinite stream of small problems; the matrices can be arbitrarily large.

How the M and AC steps differ

Each cycle in an active PE has two distinct sub-operations. M (Multiply): the two arriving tokens — one from the left (A), one from above (B) — are multiplied together. The product is a transient value, never stored to memory.

AC (Accumulate): the product is immediately added into the PE's local register C[i,j]. This register persists across cycles and grows until all K products have been summed. Hot cells label both steps so you can see which operands feed the multiply and what they accumulate into.

The role of software

The hardware array provides a fast, fixed route for operands; it does not know what problem it is solving. Software — specifically the compiler and runtime — provides everything else:

Tiling. The compiler breaks large matrices into tiles that fit the array. For a 4096×4096 multiply on a 256×256 array, it emits a loop over 16×16 = 256 tile pairs.
Injection scheduling. The staggered arrival times seen in this animation are not accidental — they are computed by the compiler. It calculates the exact cycle each operand must enter the edge so that matching A[i,k] and B[k,j] meet in cell (i,j). The hardware cannot do this itself.
DMA orchestration. A DMA engine prefetches the next tile from DRAM into an on-chip buffer while the array is busy with the current tile, hiding memory latency. The compiler schedules these transfers to keep the array fed.
Partial-sum accumulation. Each tile contributes a partial result. The runtime accumulates partial sums across tiles — either in the PE registers (if the reduction dimension fits) or in a separate on-chip buffer.
Layout decisions. The compiler chooses whether to store matrices in row-major, column-major, or a tiled layout so that tile extraction requires minimal data movement.

On a TPU, XLA performs all of the above. On a GPU using CUDA, the programmer (or a library like cuBLAS) handles tiling and scheduling explicitly. The systolic array makes the hardware simple and predictable; the complexity shifts entirely into the compiler.

Does this work for tensors?

Yes — with one key observation: almost every tensor operation of practical interest decomposes into one or more matrix multiplications, and a systolic array is purpose-built for exactly that primitive.

Batched matmul. Transformer attention computes Q·Kᵀ and scores·V independently for each head and each batch element. Each is a separate matrix multiply; the array processes them sequentially or in a pipelined schedule. The tensor rank is higher (batch × heads × sequence × d_model), but the innermost operation is always a 2-D matmul.
Convolutions. A spatial convolution over a feature map can be unrolled into a matrix multiply via im2col (each receptive field becomes a row of a matrix). The systolic array then computes the whole convolution as a single tiled matmul.
General tensor contractions. Any tensor contraction C[a,b] = Σ_c A[a,c]·B[c,b] is a matmul after appropriate reshaping. Contractions over multiple indices (like einsum) can be decomposed into a sequence of such operations.
What does not map cleanly. Operations with data-dependent access patterns ★ (sparse attention, gather/scatter, tree reductions) do not fit the fixed-route, regular-reuse model. Systolic arrays are underutilized or idle during these phases — a key reason modern accelerators pair a systolic matrix unit with a separate vector/SIMD unit for element-wise and irregular work.

★ Data-dependent access patterns — why the systolic array can't handle them

Data-dependent access patterns are memory accesses where the address you need to read depends on a value produced at runtime — not on a loop index or a static offset the compiler can compute ahead of time.

Why the systolic array cannot handle them. The array's entire efficiency argument rests on one thing: the compiler knows, before the first clock cycle fires, exactly which operand will be in which cell at which cycle. That foreknowledge is what lets it pre-schedule DMA transfers, stagger injection times, and guarantee that A[i,k] and B[k,j] arrive at PE(i,j) simultaneously. Take away that foreknowledge and the mechanism collapses — not because the hardware breaks, but because you can no longer tell it what to do.

Concrete examples:

Sparse matrix multiply. In a sparse matrix, the column index of each non-zero value is itself data stored in an index array. Before you can decide "this value goes to column j," you have to read the index. The systolic array routes by position, not by value. You would need to either pad the sparse matrix back to dense — restoring full cost and gaining nothing — or build a routing network that inspects values and redirects operands, which is a different, more expensive piece of hardware.
Gather/scatter. A gather is out[i] = in[idx[i]]. The addresses read from in are the contents of idx, which are unknown until idx is fetched. A systolic DMA prefetch requires knowing what to fetch before the array needs it. With gather indices, you cannot issue the prefetch until after you have read idx — which is exactly when the array would have needed the data to already be there. The latency hiding disappears.
Dynamic attention masking. If the set of active query-key pairs is determined by a score threshold computed from the data ("only keep the top-k attention scores"), you do not know the access pattern until after a first pass has already run. That first pass is itself a dense computation, so you end up with: dense matmul (array-friendly) → data-dependent pruning (array-hostile) → sparse matmul over pruned set (array-hostile again). The middle step breaks the pipeline.
KV-cache eviction and retrieval. During long-context LLM inference, systems that evict and page cache entries must decide which entries to keep based on attention scores — a data-dependent selection. Fetching the right K and V entries back requires indexed loads whose addresses are not known until the scores are computed. The array has to wait, or the surrounding memory system must speculate wildly.
Graph neural networks. Each node aggregates values from its neighbors. The neighbor list is graph structure — it is data. Two nodes in the same batch may have 3 neighbors and 300 neighbors respectively. The systolic array cannot be pre-programmed with a fixed schedule because the shape of the reduction is different for every node and every input graph.

The structured sparsity loophole. NVIDIA's "2:4 sparsity" (two non-zeros per four elements, enforced during training) converts an otherwise irregular pattern into a regular one. The positions of the two non-zeros per group can be encoded compactly, and the hardware can decode that encoding and redirect operands in a deterministic, compile-time-schedulable way. It is still data-dependent in the general sense, but the constraint (exactly 2-of-4, no exceptions) makes the worst case bounded and the routing circuit simple enough to embed in a modified tensor core. Fully unstructured sparsity has no such bound and cannot be handled this way.

What handles these patterns instead. Modern chips deal with this by pairing the systolic array with other execution engines: a vector/SIMD unit for element-wise ops and irregular reductions; a general-purpose load/store unit with hardware-assisted gather/scatter; and in some designs (Cerebras, Graphcore) a fully programmable router per PE that handles irregular patterns natively, at the cost of the efficiency the fixed-route model gains from regularity.

Checkable predictions of this picture

The first completed output appears only after a pipeline fill delay equal to 2(N−1)+(K−1) cycles.
Every cell performs exactly K useful MACs — no cell is skipped or double-counted.
Each A[i,k] injection is consumed once per output column it passes through (N times total). Same for B[k,j] and rows.
Switching to the naïve view shows the same MAC count with 2× the memory reads — no reuse by construction.
For a matrix–vector product, A reuse collapses (only one output column), leaving B reuse only. The array becomes partially idle — one reason systolic arrays are not universal for LLM token-by-token decode.