What Is GPU Systolic Array Device Architecture

A GPU systolic array device architecture is the specialized hybrid semiconductor microarchitecture that integrates 2D planar arrays of rhythmic, neighbor-connected processing elements (PEs)—historically pioneered by H.T. Kung and Charles Leiserson in 1978—directly into the streaming execution pipelines of modern Graphics Processing Units (GPUs). In classical GPU SIMT execution, every arithmetic operation performed by an ALU (CUDA core) requires fetching two source operands from a centralized register file and writing the computed result back through shared crossbars. While flexible for general-purpose code, this register-file-centric dataflow incurs catastrophic energy consumption on dense General Matrix Multiplications (GEMMs): in an $N \times N$ matrix multiplication, memory and register access energy scales as $O(N^3)$, rapidly thermal-throttling the GPU die. A GPU systolic array device architecture circumvents this register-energy bottleneck by embedding dedicated 2D matrix computing engines (such as modern Tensor Cores and Matrix Execution Units) within each Streaming Multiprocessor (SM). In this architecture, once activation and weight operands enter the boundary of the PE grid, they pulse rhythmically across adjacent neighbor PEs over lockstep clock cycles. Each PE performs a Multiply-Accumulate (MAC) operation and passes its input operands directly to its spatial neighbors via sub-micrometer interconnects without returning to the register file or memory hierarchy. Structurally, a GPU systolic array device architecture unites four co-designed silicon components: (1) *A 2D Mesh of Pipelined Processing Elements (PEs)*, containing dual-path multiplier-accumulator units supporting multi-precision FP16, BF16, FP8, and INT8/INT4 math; (2) *Spatial Dataflow Controllers*, implementing Weight-Stationary (WS), Output-Stationary (OS), or Input-Stationary (IS) scheduling to optimize data reuse; (3) *Tile Distribution Registers & Skew FIFOs*, formatting and staging incoming matrix slices from shared memory to synchronize diagonal wavefront propagation across the array; and (4) *Asynchronous Tensor Memory Accelerators (TMA)*, streaming multi-dimensional tensor tiles directly from High Bandwidth Memory (HBM3e/HBM4) through L2 cache into on-chip SRAM, bypassing SM register files entirely. At a high level, the architectural anatomy of a GPU systolic array and its dataflow lifecycle looks like this:

Tensor tile memory dispatch
Asynchronous TMA loads 2D sub-matrix tiles (A and B) directly from HBM into on-chip Shared Memory
Wavefront skewing & boundary staging
Skew FIFOs delay incoming matrix rows/columns to align diagonal systolic arrival wavefronts
2D systolic execution & spatial data forwardingZero-register operand reuse, localized neighbor propagation, and single-cycle MAC accumulation
Weight stationary tile pinning
Stationary weight matrix B loaded into local PE registers; remains fixed during tile compute
Horizontal activation propagation
Activation inputs A stream rhythmically eastward across adjacent PEs every clock cycle
Pipelined MAC dot-product accumulation
Each PE computes P = A · B + C; partial sum accumulates vertically down the column
Spatial neighbor forwarding
Operands forwarded directly between adjacent PEs with zero register file read/write energy
↻ Multi-GHz rhythmic systolic clock cycle
16x energy efficiency leap
Local neighbor forwarding slashes operand movement power by >90% versus SIMD ALUs
Multi-PetaFLOPS GEMM throughput
Densely packed 2D PE grids achieve near-100% MAC arithmetic utilization

A few principles are essential for understanding how GPU systolic array device architectures are structured, how spatial dataflow stationary topologies minimize register bandwidth, and why diagonal wavefront skewing is mathematically required:

Device structural hierarchy and physical classifications. Modern GPU systolic array device architectures are organized across four foundational technological domains:
1. *2D Systolic Processing Element (PE) Grid Microarchitecture:*
- Each systolic array within a GPU SM (e.g., a 4th/5th generation Tensor Core block) consists of an $M \times N$ orthogonal grid of PEs (typically $16 \times 16$ or $32 \times 32$ PEs).
- *Anatomy of a Processing Element:* A single PE contains:
- An input activation register ($A_{\text{in}}$) and horizontal forwarding register ($A_{\text{out}}$).
- A weight storage register ($B$) or vertical input register ($B_{\text{in}}$).
- A high-speed multi-precision multiplier (FP16, BF16, FP8, or INT8).
- An accumulator register ($C$) maintaining 32-bit single-precision (FP32) intermediate sum precision to prevent numerical underflow/overflow:

$$C_{\text{out}} = C_{\text{in}} + (A_{\text{in}} \times B)$$

2. *Spatial Dataflow Topologies (Stationary Modes):*
- To maximize energy efficiency, systolic architectures select a stationary dataflow that holds one matrix operand stationary inside the PE array while streaming the others:
- *Weight-Stationary (WS) Dataflow:* Matrix weights $B$ are pre-loaded into each PE's local register and held fixed throughout the computation. Activation vectors $A$ stream horizontally from the left, while partial sums $C$ stream vertically downward, accumulating products along the column. WS dataflow is ideally suited for convolutional layers and large LLM prefill GEMMs where weight reuse is maximized.
- *Output-Stationary (OS) Dataflow:* The accumulator $C$ remains stationary inside each PE throughout the calculation, accumulating products over time. Matrix $A$ streams horizontally and Matrix $B$ streams vertically. Once accumulation finishes, the completed result $C$ is drained. OS minimizes accumulator movement energy.
- *Input-Stationary (IS) Dataflow:* Activations $A$ remain pinned in the PEs while weights and partial sums move, optimal for streaming real-time batch-1 inference.
3. *Wavefront Synchronization & Input Skewing FIFOs:*
- Because operands take exactly one clock cycle to propagate to the adjacent PE, feeding matrix rows and columns simultaneously would cause misaligned product accumulation.
- *Skew FIFOs:* Matrix inputs must be skewed diagonally in time. For an $N \times N$ systolic array, row $i$ of Matrix $A$ is delayed by $i$ clock cycles, and column $j$ of Matrix $B$ is delayed by $j$ clock cycles before injection at the array boundary:

$$T_{\text{arrival}}(i, j) = i + j + 1$$

This guarantees that activation $A_{i,k}$ and weight $B_{k,j}$ meet simultaneously at PE $(i, j)$ at clock cycle $k + i + j$.
4. *Memory Hierarchy Integration & Asynchronous TMA Engine:*
- Traditional GPU SIMD pipelines load data from L1/Shared Memory into the register file before feeding ALUs. GPU systolic arrays deploy dedicated *Tensor Memory Accelerators (TMA)*.
- The TMA engine executes asynchronous hardware-managed transfers, copying multi-dimensional tensors directly from HBM3e/HBM4 through L2 into Shared Memory using tensor coordinate descriptors without consuming SM execution threads or register slots, maintaining continuous systolic saturation.

The microarchitectural floorplan of a GPU Streaming Multiprocessor featuring embedded Systolic Tensor Arrays and the internal anatomy of a 2D Processing Element (PE) are illustrated below:

GPU Systolic Array Device Architecture: Tensor Engine & PE Grid Integration of 2D systolic matrix multiply arrays inside GPU Streaming Multiprocessors SM Microarchitecture with Systolic Engines SIMT Warp Schedulers & Scoreboards Dispatches asynchronous WMMA / MMA instructions to Systolic Arrays 256KB SM Register File (Bypassed during systolic dataflow) SIMD Core Array 32× FP32 CUDA Cores 32× INT32 Cores 1 MAC / cycle / core Systolic Tensor Array (16×16) 256 Pipelined MAC PEs FP8, FP16, BF16, INT8 Matrix GEMM 256 MACs / cycle (Zero RF energy) 228KB Shared Memory & Async TMA Unit Streams matrix tiles directly from HBM3e into systolic boundary FIFOs SIMD ALUs (Gray) | Systolic Tensor Array (Gold) | SRAM/TMA (Green) Systolic Dataflow & Energy Scaling Math 1. Peak Systolic Throughput per Array: Throughput = 2 · N_pe · f_clock = 2 · (M · N) · f_clock 16×16 PE array @ 2.0 GHz = 1024 GFLOPS / SM (FP16/BF16) 2. Operand Movement Energy Scaling: E_simd ∝ O(N^3 · E_rf) vs. E_systolic ∝ O(N^2 · E_sram + N^3 · E_pe) Local neighbor forwarding eliminates >90% of register access energy 3. Diagonal Wavefront Execution Latency: Cycles_total = M + N + K - 2 [K = matrix inner dimension] Pipelined execution completes K dot products across full 2D array Array Size: 16×16 | PEs: 256 | Precision: FP8/FP16/INT8 | Accum: FP32 Spatial Data Reuse Neighbor-to-neighbor forwarding Weight-Stationary Dataflow Pins weights; streams activations Asynchronous TMA Bypasses SM registers from HBM GPU systolic array device architecture marries general-purpose SIMT flexibility with dedicated 2D matrix arrays slashing data movement energy and driving multi-PetaFLOPS throughput for dense deep learning workloads

In addition to SM-level integration, the internal 2D spatial dataflow mechanics (Weight-Stationary diagonal wavefront propagation and PE circuit design) are essential. The 2D systolic mesh and PE circuit schematics are illustrated below:

2D Systolic Wavefront Dataflow & Processing Element (PE) Weight-stationary matrix multiplier grid with diagonal input skewing and localized accumulation Weight-Stationary 2D Mesh Dataflow Col 0 Col 1 Col 2 A[0,k] → PE(0,0) W[0,0] PE(0,1) W[0,1] PE(0,2) W[0,2] A[1,k] → PE(1,0) W[1,0] PE(1,1) W[1,1] PE(1,2) W[1,2] A[2,k] → PE(2,0) W[2,0] PE(2,1) W[2,1] PE(2,2) W[2,2] ↓ C[0,j] ↓ C[1,j] ↓ C[2,j] Horizontal: Activations (Gold) | Vertical: Partial Sums (Green) Individual Processing Element (PE) Circuit Ain Register From West PE Weight Reg (B) Stationary (Pinned) Multiplier Unit P = Ain × B (FP16/FP8) Cin Register From North PE Adder / Accumulator Cout = Cin + P (FP32) Aout Register To East PE → Cout Register (To South PE ↓) Single-cycle MAC operation with local register-free bypass Systolic Breakthrough: Diagonal wavefront propagation reuses operands across 256 PEs without touching the primary register file 2D systolic arrays replace power-hungry centralized register crossbars with rhythmic spatial neighbor forwarding delivering near-100% MAC arithmetic density while keeping energy consumption within package thermal limits

Boundary fill-drain pipeline overheads, array dimension scaling ($N \times N$), and numerical precision quantization. Designing high-performance GPU systolic array device architectures requires balancing structural hardware trade-offs:
- *Fill and Drain Pipeline Latency Overhead:* A 2D systolic array of dimensions $M \times N$ computing an inner dimension $K$ requires $M + N - 2$ cycles just to fill the array with the first diagonal wavefront, and another $M + N - 2$ cycles to drain the final accumulated partial sums. The effective compute utilization ($\eta$) is governed by:

$$\eta = \frac{K}{K + M + N - 2}$$

If the inner matrix dimension $K$ is small (e.g., $K = 16$ in small batch sizes), $\eta < 50\%$, and the systolic array spends half its time idling during fill/drain phases. Consequently, GPU systolic arrays operate at peak efficiency on large GEMMs where $K \ge 128\text{–}256$, achieving $>95\%$ continuous MAC utilization.
- *Array Dimension Scaling ($16 \times 16$ vs. $128 \times 128$):* While dedicated NPUs (such as Google TPU's $128 \times 128$ or $256 \times 256$ Matrix Multiply Units) favor massive monolithic arrays, GPUs deploy smaller, modular $16 \times 16$ systolic Tensor Core arrays. Smaller arrays minimize fill/drain latency penalties on diverse non-GEMM workloads, reduce silicon dark-silicon penalties, and allow flexible tiling across thousands of independent thread warps executing concurrently.
- *Microscaling & Multi-Precision Quantization (FP8, MXFP6, MXFP4):* As deep learning models scale into hundreds of billions of parameters, 16-bit systolic floating-point units become memory-bandwidth and area constrained. Modern GPU systolic arrays incorporate *Microscaling (MX)* formats. Instead of standard IEEE FP16, groups of 32 elements share a single 8-bit scale factor (E8M0) while individual weights and activations are stored in 4-bit (E2M1) or 6-bit (E3M2) floating-point formats. The systolic multiplier units compute 4-bit multiplications at $4\times$ the throughput and $8\times$ the energy efficiency of FP16, with intermediate FP32 accumulators preserving gradient fidelity.

Read a GPU systolic array device architecture through a *spatial-dataflow-energy-minimization* lens: a specialized computing microarchitecture where rhythmic, neighbor-to-neighbor operand forwarding replaces power-hungry centralized register files to maximize dense matrix arithmetic throughput.

Take what is gpu systolic array device architecture further

Ask the copilot about this term, or have our engineers assess it against your process.

Ask CFSGPT about this → Request a technical assessment → Open the CFS app