Scalable Parameterized Systolic Array Matrix Multiplier (FPGA/RTL)

A fully parameterized N×N systolic array for hardware matrix multiplication, written in Verilog. Built solo for Problem Statement 2 of SSPI, the signal-processing competition at Spectrum (IIT (ISM) Dhanbad's electronics dept fest) — 1st place winner.

The design scales from 2×2 to 16×16 (and beyond, in principle) purely through parameters — no structural rewrite needed — and is verified with a self-checking testbench against a software golden model across 9 test cases per configuration, achieving 100% mathematical accuracy.

Abstract

This project implements the architectural design, Verilog RTL, and functional verification of a scalable, parameterized systolic array matrix multiplier, computing A × B = C in hardware. The design addresses the classical von Neumann memory bottleneck of scalar matrix multiplication by using localized, rhythmic data propagation between processing elements instead of centralized memory access. It was verified systematically from a 2×2 baseline up to a 16×16 grid, with simulation confirming 100% mathematical accuracy and predictable O(N) latency scaling — making the architecture well suited to high-throughput DSP, control systems, and machine learning accelerator front-ends.

Design philosophy: why a systolic array?

Standard matrix multiplication on a scalar CPU is memory-bound — every multiply-accumulate requires a fetch and a store, and the operation is inherently O(N³). This project's architecture sidesteps that limitation through localized data propagation: once data enters the array at its edges, it flows seamlessly between adjacent Processing Elements (PEs) without ever returning to a central memory bank. Each PE performs one MAC per cycle and passes its operands on to its neighbors, so throughput scales with the size of the array rather than being capped by memory bandwidth.

Architectural novelty

The core design decision that shapes this whole project is that the hardware is fully parameterized — the matrix dimension N is abstracted away from the physical wiring entirely. Rather than hardcoding structural connections for a specific size, the top module uses Verilog generate blocks to dynamically unravel the 2D routing mesh at elaboration time. The practical result: this is a genuine, reusable IP core. Changing N from 2 to 16 requires touching exactly one parameter — no rewiring, no manual instantiation, no structural changes to pe.v or systolic_array_top.v at all.

Static timing analysis (STA) optimization

From a physical-design standpoint, timing closure at high clock frequencies was a first-class design constraint, not an afterthought. The combinational logic depth inside each PE is deliberately kept to a single multiply-add, and pipeline registers physically isolate every Processing Element — every input and output of pe.v is registered. This means the critical path length is independent of array size: a 16×16 array has the same per-stage combinational depth as a 2×2 array, so scaling N up doesn't degrade the achievable clock frequency the way a less disciplined design would.

Architecture

                    B[0][j] B[1][j] B[2][j] ... (streamed from north, skewed)
                       │       │       │
   A[i][0] ──────────► PE ──► PE ──► PE ──► ...
                       │       │       │
   A[i][1] ──────────► PE ──► PE ──► PE ──► ...
                       │       │       │
   A[i][2] ──────────► PE ──► PE ──► PE ──► ...
                       │       │       │
                      ...     ...     ...

The Processing Element (pe.v)

The basic computational node is a pipelined MAC unit. On every positive clock edge, it executes two operations concurrently:

• Accumulation — multiplies the current a_in and b_in, and adds the product into an internal 64-bit sum_out register: sum_out <= sum_out + (a_in * b_in).
• Data propagation — passes the unmodified a_in to its east neighbor (a_out) and b_in to its south neighbor (b_out), both one cycle later.

Every signal in and out of the PE is registered — there's no combinational pass-through anywhere in the datapath, which is what keeps the STA story simple regardless of array size.

Top module & data skewing (systolic_array_top.v)

The top module instantiates an N×N grid of PEs using generate blocks, so N, the data width, and the accumulator width are all Verilog parameters. Because data only moves exactly one row or column per clock cycle, the inputs can't be injected in raw row/column order — they have to be structurally skewed before they enter the array:

• Row i of matrix A is delayed by i clock cycles before injection.
• Column j of matrix B is delayed by j clock cycles before injection.

This staggering ensures that A[i][k] and B[k][j] physically arrive at the same PE, PE(i,j), on the same clock cycle — which is the entire mechanism that makes the systolic dataflow correct. In this project the skewing is generated by the testbench's apply_skewed_stream task, standing in for what a real streaming input controller would do upstream of the array in a full system.

Ports are implemented as flattened, packed buses (inp_west, inp_north, each N*DATA_W bits wide) rather than unpacked 2D arrays, since Verilog-2001 module ports don't support unpacked arrays — each lane is sliced out with [(N-1-i)*DATA_W +: DATA_W].

Latency analysis

The latency of the systolic pipeline is mathematically deterministic, not something that has to be measured empirically after the fact:

1. Stream injection — it takes 2N − 1 clock cycles to fully stream the skewed data elements into the edges of the array.
2. Propagation & flush — an additional N + 1 cycles are needed for the final accumulation to complete at the deepest node, PE(N−1, N−1).

That gives a total of 3N + 1 cycles until done goes high — confirmed by simulation for every tested array size:

N	Cycles	Latency @ 10ns clock
2	7	75 ns
4	13	135 ns
6	19	195 ns
8	25	255 ns
16	49	495 ns

This is the practical payoff of the architecture: O(N) time complexity at the array boundary, versus the O(N³) op count of the underlying computation — a real hardware acceleration advantage over scalar processing, and the number that actually matters when justifying "why build this in hardware" in the first place.

Verification

tb/tb_systolic_array.v is a self-checking testbench (built and run in ModelSim) that:

1. Computes a software golden reference (GOLDEN[i][j] = Σ A[i][k]·B[k][j]) independently for each test case.
2. Streams operands into the array with correct skewing via apply_skewed_stream.
3. Waits for done, then reads every PE's sum_out register directly via hierarchical reference and compares it against the golden model.
4. Repeats this across 9 test cases per array size: small hand-picked integers, an all-zero matrix, identity × identity, boundary/max values (0xFF, to stress-test the 64-bit accumulator against overflow), and 5 rounds of randomized stress testing.
5. Prints a per-element PASS/FAIL scoreboard, plus a grand summary with total accuracy.

Result: 100% match against the golden model across all test cases, verified for N = 2, 4, 6, 8, and 16.

A notable implementation constraint worth mentioning: ModelSim (the version used here) doesn't allow variable-index access into generate-block hierarchies at simulation time. That means the output-capture routine that reads every PE's sum_out can't be written as a simple nested loop — it has to be manually unrolled into explicit hierarchical references (uut.pe_row[i].pe_col[j].u_pe.sum_out) for every (i,j) pair, once per array size tested. It's not elegant, but it's a real constraint of the toolchain rather than a design choice, and the full unrolled listings for each size are preserved in the design report's appendix.

Note on hardware status: this was verified through RTL simulation only, not yet synthesized or run on a physical FPGA. Timing, resource utilization (LUTs/DSPs/BRAMs), and post-synthesis achievable frequency are not yet characterized — see Future Work.

Simulation waveforms

ModelSim waveforms for each array size, showing the clock, skewed input streams, the done flag, and accumulated PE outputs matching the expected products.

N = 2 — 7 cycles, 75 ns latency

N = 4 — 13 cycles, 135 ns latency

N = 6 — 19 cycles, 195 ns latency

N = 8 — 25 cycles, 255 ns latency — successfully processes maximum hexadecimal values (0xFF) inside 64-bit accumulators without overflow

N = 16 — 49 cycles, 495 ns latency — full 16×16 grid, 256 active PE nodes

Conclusion

The implemented systolic array matrix multiplier meets and exceeds its original design goals. By parameterizing the top-level architecture, the design scales natively from 2×2 up through 16×16 configurations without any structural rewrite. Localized MAC processing fully decouples throughput from memory access, while per-PE pipeline segmentation keeps the physical timing story simple and predictable regardless of array size.

Repo structure

.
├── rtl/
│   ├── pe.v                    # Processing element (pipelined MAC)
│   └── systolic_array_top.v    # Parameterized N×N array + skew/latch logic
├── tb/
│   └── tb_systolic_array.v     # Self-checking testbench, golden-model comparison
├── sim_results/
│   └── waveform_N*.png         # ModelSim waveform screenshots, N = 2, 4, 6, 8, 16
├── docs/
│   └── Design_Report.pdf       # Full write-up: architecture, timing analysis, waveforms, appendix
└── README.md

Running the simulation

Requires a Verilog simulator (tested on ModelSim/QuestaSim; should also work on Icarus Verilog with minor tweaks since it avoids SystemVerilog-only constructs).

# ModelSim / QuestaSim
vlib work
vlog tb/tb_systolic_array.v          # pulls in systolic_array_top.v and pe.v via `include
vsim -c tb_systolic_array -do "run -all; quit"

# Icarus Verilog (open source alternative)
iverilog -o sim.out tb/tb_systolic_array.v -I rtl
vvp sim.out

To change the matrix size, edit localparam N in tb/tb_systolic_array.v — the RTL itself requires no changes. Note that the testbench's capture_outputs task is unrolled per array size (see Verification above), so changing N to a value other than 2, 4, 6, 8, or 16 will require writing a new unrolled capture block for that size — the full worked examples for each existing size are in the design report's appendix as a reference.

Future work

• Synthesize on a real FPGA target (e.g. Xilinx Artix-7 / Intel Cyclone) and report LUT/DSP/BRAM utilization and max achievable Fmax.
• Add an AXI-Stream or simple handshake wrapper so the array can be driven by a real memory-mapped system instead of a testbench-generated skew.
• Explore weight-stationary / output-stationary dataflow variants to compare against the current input-streamed design, and reduce idle PE cycles during fill/drain.
• Reduce accumulator width dynamically based on DATA_W and N instead of a fixed 64 bits, to save area for smaller configurations.
• Generalize the testbench's output-capture routine to avoid manual unrolling per array size (e.g. via a SystemVerilog interface array or a different simulator without the hierarchical-indexing restriction).

Background

Built and submitted solo for Problem Statement 2 of SSPI (Spectrum, IIT (ISM) Dhanbad's electronics department fest) — awarded 1st place. Full architectural rationale, pipeline/latency derivation, and simulation waveforms for N = 2, 4, 6, 8, 16 are in docs/Design_Report.pdf.

License

MIT — see LICENSE.