2D Systolic Array Accelerator

A parameterizable Weight-Stationary 2D Systolic Array implemented in Verilog for hardware-accelerated matrix multiplication — the same architecture at the heart of Google's TPU.

Inspired by Systolic Arrays for VLSI — Kung & Leiserson (1978).

---

What It Does

Computes C = A × B for N×N matrices using an array of N×N Processing Elements (PEs), each performing a single Multiply-Accumulate (MAC) operation per clock cycle.

      B[0]  B[1]  B[2]
       ↓     ↓     ↓
A[0]→[PE]→[PE]→[PE]→
       ↓     ↓     ↓
A[1]→[PE]→[PE]→[PE]→
       ↓     ↓     ↓
A[2]→[PE]→[PE]→[PE]→
       ↓     ↓     ↓

• Activations (rows of A) stream in from the left, flow right
• Weights (matrix B) are preloaded into each PE before compute (weight-stationary)
• Partial sums accumulate downward, column by column

---

Architecture

Dataflow: Weight-Stationary

This array implements the weight-stationary dataflow — the same strategy used by Google's TPU v1. In this scheme the weight matrix B is loaded into the PE array once and remains stationary for the entire multiplication. Activations (matrix A) and partial sums stream through the array:

Data	Loaded / Streamed	Direction	Lifetime
Weights	Preloaded (stationary)	—	Held in PE registers for entire compute phase
Activations	Streamed each cycle	Left → Right	Each value is forwarded to the right neighbour after one cycle
Partial sums	Accumulated	Top → Bottom	Grows as it passes through each row of PEs

The computation proceeds in three phases:

1. Weight Load — Assert weight_load for N consecutive cycles. On cycle t, set weight_row_sel = t and drive weight_data with row t of matrix B. Each PE latches its weight into a local register.
2. Compute — Deassert weight_load and begin streaming the columns of A through act_col_in, one column per cycle. The built-in skew registers (when SKEW_EN=1) delay row r by r cycles so that the diagonal systolic wavefront aligns correctly.
3. Output Drain — Valid partial sums emerge at psum_bottom_out after the pipeline is fully filled. Column c of the result matrix C exits the bottom of column c of the array.

The annotated waveform below shows all four phases for the N=4 Incremental test case — weight loading (orange), skew chain delay (cyan), activation wavefront propagation (green), and partial sum accumulation cascade (pink):

Processing Element (PE)

Each PE is a single-cycle Multiply-Accumulate (MAC) unit with a stationary weight register:

if weight_load:
    weight  ← weight_in          // latch new weight

psum_out ← psum_in + (act_in × weight)   // MAC
act_out  ← act_in                         // forward activation to right neighbour

All operations are fully pipelined — one MAC per PE per clock cycle.

Input Skew

Row r of activations is internally delayed by r flip-flop stages, aligning the systolic wavefront without requiring external skewing logic (SKEW_EN=1 by default). This creates the characteristic diagonal wavefront:

Cycle   Row 0 input   Row 1 input   Row 2 input   Row 3 input
  0       A[0][0]          0             0             0
  1       A[0][1]       A[1][0]          0             0
  2       A[0][2]       A[1][1]       A[2][0]          0
  3       A[0][3]       A[1][2]       A[2][1]       A[3][0]
  4          0          A[1][3]       A[2][2]       A[3][1]
  5          0             0          A[2][3]       A[3][2]
  6          0             0             0          A[3][3]

Array Size & Precision Support

The design is fully parameterizable at elaboration time:

Parameter	Default	Description
N	4	Array dimension (N×N PEs, computes N×N matmul)
DATA_W	8	Bit-width of activations and weights (unsigned)
PSUM_W	20	Bit-width of partial sums (2 × DATA_W + 4)
SKEW_EN	1	Enable built-in input skew registers

Tested configurations:

N	PEs	act_col_in width	psum_bottom_out width	Max product per PE	Max accumulated psum	Status
1	1	8 bits	20 bits	65 025	65 025	✅ Pass
2	4	16 bits	40 bits	65 025	130 050	✅ Pass
4	16	32 bits	80 bits	65 025	260 100	✅ Pass

The PSUM_W default of 2 × DATA_W + 4 provides 4 extra bits of headroom beyond the minimum product width, supporting accumulation of up to 2⁴ = 16 maximum-value products without overflow (i.e. N ≤ 16 with 8-bit data).

To use different precisions, override DATA_W (and optionally PSUM_W) at elaboration:

cmake -GNinja .. -DN=8 -DDATA_W=16 -DPSUM_W=40

Timing

Phase	Duration
Weight load	N cycles
Pipeline fill	N − 1 cycles
Compute	N cycles
Total latency	2N − 1 cycles from first activation to last valid output

---

Verification

The design is verified using Verilator (C++ testbench) across three array sizes and four test cases:

Test Case	Description
All Ones	A = B = all 1s
Identity	A = B = identity matrix
Incremental	A[i][j] = i·N+j+1, B = 2·A
Near Overflow	A = B = all 255 (max 8-bit value)

All test cases pass for N = 1, 2, and 4.

VCD waveform dumps are generated for each test case and can be viewed in GTKWave. A pre-configured demo.gtkw save file is included for quick waveform inspection of the N=2 Incremental case.

Synthesised Layout

The images below show the synthesised chip layout, illustrating the regular, tiled structure of the PE grid:

---

Design Notes

Overflow Headroom

The default PSUM_W = 2 × DATA_W + 4 gives 4 bits of accumulation headroom. For 8-bit unsigned data, the worst-case single product is 255 × 255 = 65 025 (fits in 16 bits). With 4 extra bits, the accumulator can hold up to N = 16 such products (65 025 × 16 = 1 040 400, which fits in 20 bits). For larger N or wider data, increase PSUM_W accordingly — a safe rule of thumb is PSUM_W ≥ 2 × DATA_W + ⌈log₂(N)⌉.

Skew Trade-offs

The built-in skew registers (SKEW_EN=1) simplify the external interface by handling activation alignment internally, at the cost of N × (N−1) / 2 additional flip-flops. For an N=4 array this is only 6 extra 8-bit registers — negligible. For very large arrays where register budget is tight, set SKEW_EN=0 and implement skewing in the data feeder or memory controller.

Weight Loading Protocol

Weights are loaded one row at a time using a row-select signal (weight_row_sel). This means the weight bus width scales as O(N × DATA_W) rather than O(N² × DATA_W), keeping the I/O pin count manageable. The trade-off is that weight loading takes N cycles instead of 1. For workloads where weights change frequently, consider double-buffering the weight registers (not implemented here).

Scalability Considerations

• Area: Scales as O(N²). Each PE contains one DATA_W-bit multiplier, one PSUM_W-bit adder, and three registers.
• Frequency: The critical path runs through one multiply-add inside a single PE. Since there are no combinational chains across PE boundaries, the design should achieve high clock frequencies.
• Memory bandwidth: The array consumes N activations per cycle (left edge) and produces N partial sums per cycle (bottom edge), giving O(N) bandwidth vs. O(N²) computation — the classic systolic advantage.
• Utilization: For matrices smaller than N×N, the outer PEs will be idle. Tiling larger matrices across multiple passes is left to the host controller.

Verification Philosophy

The testbench computes a golden reference in C++ (mirroring the NumPy golden model) and checks every intermediate partial sum at the exact cycle it becomes valid — not just the final output. This cycle-accurate checking validates the systolic timing, the skew alignment, and the MAC pipeline simultaneously. Any off-by-one error in the dataflow would cause a mismatch.

---

Build & Run

Prerequisites

• Verilator ≥ 4.200
• CMake ≥ 3.16
• A C++17-capable compiler (GCC / Clang)

Run a single size

mkdir build && cd build
cmake -GNinja .. -DN=4   # set N=1, 2, or 4
ninja
./VArray_tb

Run all sizes (N = 1, 2, 4)

bash test_all.sh

This builds and runs all three configurations, producing VCD waveform files in the project root.

---

Project Structure

.
├── pe.v                  # Processing Element module
├── systolic_array.v      # Top-level N×N systolic array
├── systolic_array_tb.cc  # Verilator C++ testbench
├── golden_model.py       # NumPy golden reference
├── CMakeLists.txt        # Verilator CMake build
├── test_all.sh           # Runs all N configurations
└── graphics/             # Architecture diagrams and waveform screenshots

---

References

• DiP: A Scalable, Energy-Efficient Systolic Array for Matrix Multiplication Acceleration :
• Kung, H. T., & Leiserson, C. E. (1978). Systolic Arrays for VLSI. CMU-CS-79-103.
• Verilator Documentation
• Google TPU Paper — Jouppi et al., ISCA 2017