Parameterized Systolic Matrix-Multiply Accelerator

SystemVerilog · cocotb · Icarus Verilog · Vivado 2019.2 · Xilinx ZCU104 (xczu7ev)

A matrix multiply accelerator built three ways, each one fixing the previous one's limitation. Same AXI4-Stream interface on all three, same cocotb testbench against a Python golden model under randomized backpressure, then characterized through synthesis and implementation to see what each design actually costs.

Version	Compute units	Compute latency	Operand storage	Dimensions
sequential/	1 MAC	~N³ cycles	registers	NxN
systolic_reg/	N² PEs	~3N cycles	registers	NxN
systolic_bram/	M·K PEs	~M+N+K cycles	banked memory	MxN · NxK

Out-of-context synth and impl, XCZU7EV-2, Vivado 2019.2. Fmax from worst negative slack against a 500 MHz constraint. Square configs to isolate scaling.

N	sequential	systolic_reg	systolic_bram
	LUT / FF / Fmax	LUT / FF / Fmax	LUT / FF / Fmax
4	422 / 316 / 334 MHz	1,982 / 801 / 389 MHz	1,988 / 595 / 350 MHz
8	710 / 1,105 / 328 MHz	7,547 / 3,268 / 351 MHz	7,408 / 2,322 / 398 MHz
16	1,993 / 4,243 / 260 MHz	29,740 / 13,336 / 305 MHz	28,860 / 9,333 / 359 MHz
32	7,348 / 16,709 / 214 MHz	121,739 / 54,176 / 238 MHz	116,305 / 38,199 / 276 MHz
64	26,882 / 66,239 / 178 MHz	exceeds device	exceeds device

• Banking wins on both axes: fewer LUTs and ~30% fewer flops than systolic_reg at every N≥8, and it closes higher.
• Fmax drops monotonically with size on all three. Sequential falls hardest (334 → 178 MHz) because its operand storage and output muxing scale as N² even though its compute doesn't.
• N=64 needs ~464k LUTs against 230,400 available. Both systolic variants fail identically, so the design is compute-bound not storage-bound: the N² PE count dominates, and banking improves efficiency without moving the ceiling.
• No DSP inference at WIDTH=8. Vivado maps 8x8 multiplies to LUT and CARRY8 instead of DSP48E2.

Forcing ram_style = "block" confirms the banking maps cleanly onto real primitives (BRAM count is exactly 2N, one per bank) but loses at these sizes:

N	default (LUTRAM)	forced BRAM	Fmax delta
4	1,988 LUT, 0 BRAM, 350 MHz	2,145 LUT, 8 BRAM, 307 MHz	-43 MHz
8	7,408 LUT, 0 BRAM, 398 MHz	7,607 LUT, 16 BRAM, 309 MHz	-89 MHz
16	28,860 LUT, 0 BRAM, 359 MHz	28,640 LUT, 32 BRAM, 290 MHz	-69 MHz
32	116,305 LUT, 0 BRAM, 276 MHz	116,101 LUT, 64 BRAM, 253 MHz	-24 MHz

A bank holds N words of 8 bits, so 128 bits at N=16 against an 18 kbit RAMB18. Each block sits under 1% used and the LUT savings never show up. Below N=8 area goes up, since the address and control logic costs more than the storage it removes. Banking would only pay off in block RAM with much deeper banks.

Systolic dataflow

For an MxN · NxK product the array is M rows by K columns; each PE does N multiply-accumulates over the shared dimension. Shown 3x3.

                 Matrix B streams DOWN (skewed by column)

           b_edge[0]   b_edge[1]   b_edge[2]
               |           |           |
               v           v           v
a_edge[0] --> +-----+ --> +-----+ --> +-----+
(A row 0)     |PE 00|     |PE 01|     |PE 02|
              +-----+     +-----+     +-----+
               |           |           |
               v           v           v
a_edge[1] --> +-----+ --> +-----+ --> +-----+
(A row 1)     |PE 10|     |PE 11|     |PE 12|
              +-----+     +-----+     +-----+
               |           |           |
               v           v           v
a_edge[2] --> +-----+ --> +-----+ --> +-----+
(A row 2)     |PE 20|     |PE 21|     |PE 22|
              +-----+     +-----+     +-----+

Matrix A streams RIGHT (skewed by row)

• Output-stationary: PE(i,j) owns C[i][j] and does one MAC per cycle.
• Each PE registers its passthroughs, so every hop costs a cycle.
• The skew (A row i enters i cycles late, B column j enters j cycles late) falls out of that per-hop register, not a scheduled signal. Correct global timing from a purely local rule, which is why systolic arrays scale: communication is nearest-neighbour only, so wire length stays constant as the array grows.

Layout

sequential/         1 MAC baseline
systolic_reg/       NxN PE grid, operands in registers
systolic_bram/      same grid, operands in banked memory, MxN * NxK
    rtl/
        chip.sv         top: AXI4-Stream slave and master
        control.sv      FSM, LOAD -> COMPUTE -> LOAD
        datapath.sv     operand storage, edge feeding, result streaming
        array.sv        PE grid and interconnect
        pe.sv           one MAC + passthrough registers
        bram.sv         banked operand memory
    tb/
        golden_model.py reference triple-nested multiply
        testbench.py    cocotb driver
fpga/
    constraints/    clk.xdc
    scripts/        run_ooc.tcl, sweep.sh
    results/        per-config timing, utilization, power
docs/
    log.md          what I tried and why

Quickstart

Requires: Icarus Verilog, cocotb, Vivado 2019.2

1. Simulate (100 random trials with per-cycle backpressure, against the golden model):

    cd systolic_bram      # or sequential/ or systolic_reg/
    make

2. One synthesis config:

    vivado -mode batch -source fpga/scripts/run_ooc.tcl -tclargs systolic_bram 8

3. Full sweep:

    bash fpga/scripts/sweep.sh

Roadmap

Part 1:

• Sequential baseline
• Output-stationary systolic array
• Banked operand memory, generalized to MxN · NxK
• cocotb regression against a golden model under backpressure
• PPA sweep across three architectures and five sizes
• Forced-BRAM ablation

Part 2:

• Deploy to hardware: package as IP, block design with the PS and AXI DMA, measure end-to-end throughput
• Deeper banks (wider elements or tiled operands) where block RAM wins
• Edge cases: identity, all-max (accumulator width stress), all-zeros