从零复现脉动阵列：FPGA上实现矩阵乘法加速器

1. 为什么值得动手复现一个脉动阵列矩阵乘法这件事放在CPU上跑大家习以为常——三层循环编译器自动向量化顶多再挂个OpenMP。但当你真正把矩阵规模拉到几百乘几百还要在功耗和延迟双重约束下跑推理任务时CPU的通用架构就开始露怯了。谷歌TPU第一代之所以能在推理场景打出漂亮的能效比核心武器之一就是脉动阵列Systolic Array——一种让数据像心跳一样有节奏地在PEProcessing Element阵列中流动的计算结构。我第一次在FPGA上复现脉动阵列是为了做一个定点矩阵乘法加速器目标规模是32×32的INT8矩阵。当时踩了不少坑PE之间的握手时序对不上、边界数据注入错位、综合后时序跑不到目标频率。后来一点点啃下来发现这东西的原理其实不复杂难的是把数据流动的节奏在硬件里精确对齐。这篇文章就是把我从零复现脉动阵列的完整过程拆开讲。适合谁看有Verilog基础、做过至少一个FPGA小项目比如UART、计数器、状态机想往AI加速器方向深入的人。能学到什么脉动阵列的核心原理、PE的微架构设计、数据调度策略、Verilog实现细节、仿真验证方法以及综合布线阶段的实操经验。前置知识组合逻辑与时序逻辑的区别、阻塞赋值与非阻塞赋值的语义、基本的定点数概念。如果这些还不太熟建议先把Verilog的计数器、移位寄存器、简单状态机写一遍再回来。我不会只贴代码而是把每个设计决策背后的为什么讲清楚——为什么PE要这样设计、为什么数据要斜着注入、为什么累加器位宽要留够。这些才是你真正能迁移到其他项目里的东西。2. 脉动阵列到底在算什么原理拆解与方案选型2.1 从一个人搬砖到流水线协作先打个比方。普通矩阵乘法像是一个人搬砖你有一个计算单元它要反复从内存里取A的元素、取B的元素、乘一下、加到结果里然后再取下一对。取数据的时间远大于计算时间计算单元大部分时间在等数据。脉动阵列则是一队人排成方阵传砖每个PE只负责一次乘加算完就把数据传给右边的邻居同时从左边接收新数据。数据像水流一样穿过整个阵列每个PE在每个时钟周期都在干活。这就是脉动这个名字的由来——数据有节奏地脉动式流动。具体到矩阵乘法C A × B假设A是M×KB是K×N那么A的每一行从阵列左侧逐周期注入向右流动B的每一列从阵列顶部逐周期注入向下流动每个PE(i,j)在时刻t接收到A[i][t]和B[t][j]做乘加累加到本地寄存器经过K个周期后PE(i,j)的累加器里就是C[i][j]关键在于A的第i行第k个元素和B的第k列第j个元素必须在同一时刻到达PE(i,j)。这就是数据调度要解决的核心问题。2.2 为什么选输出固定型Output-Stationary脉动阵列有好几种数据流变体常见的有数据流类型特点适用场景输出固定Output-Stationary结果留在PE本地累加A和B流过通用矩阵乘法最直观权重固定Weight-StationaryB矩阵预加载到PEA流过推理场景权重可复用行固定Row-StationaryA的行固定B和部分和流动卷积神经网络我选的是输出固定型理由很实际第一它最容易理解PE的结构最简单——一个乘法器加一个累加器第二对于教学和首次复现来说调试难度最低因为每个PE的结果就在本地不需要追踪流动的部分和第三它天然支持任意K值只要控制好注入周期数就行。权重固定型虽然在实际推理芯片里更常见因为权重可以预加载、反复使用但它需要额外的权重加载阶段和更复杂的控制逻辑。第一次做没必要给自己上难度。2.3 阵列规模怎么定32×32还是8×8阵列规模直接决定了资源消耗和设计复杂度。我一开始想直接上32×32结果综合报告一出来DSP用了1024个目标器件根本放不下。后来退回到8×8先跑通再逐步扩展到16×16。这里有个经验公式每个PE至少消耗1个DSP乘法器和1个累加器若干LUTFF。8×8阵列就是64个DSP16×16是256个32×32是1024个。你先查一下自己板子上FPGA的DSP数量留出30%余量给其他逻辑就知道能上多大阵列了。另外阵列不一定是正方形。如果你的矩阵是M×N且M远大于N可以考虑M×N的非方阵列但控制逻辑会复杂一些。我建议第一次做就选正方形8×8起步。3. PE的微架构设计一个乘法器加一个累加器就够了吗3.1 PE的基本结构一个最简PE需要以下端口输入a_in来自左侧、b_in来自上方、clk、rst_n、en使能输出a_out传给右侧、b_out传给下方内部acc累加器寄存器每个时钟周期当en有效时a_out a_inb_out b_in数据继续流动acc acc a_in * b_in乘加看起来很简单但有几个细节必须处理好。3.2 数据位宽与定点数格式我用的是INT8输入INT32累加。为什么输入用INT8这是量化推理的常见格式8位乘法器在FPGA里映射到DSP48 slice非常高效累加用INT32K个INT8乘积相加最坏情况下每个乘积是127×127≈16129K256时总和约4.1M需要23位。留到32位有充足余量也方便后续做移位和饱和处理如果你要做浮点那DSP消耗会翻好几倍而且时序很难收敛。第一次做老老实实用定点。注意INT8乘法在Verilog里要显式声明为有符号数signed否则综合工具可能按无符号处理结果全错。我在这上面浪费过一整个下午。3.3 累加器的复位与清零策略累加器不能一直累加必须在每轮矩阵乘法开始前清零。有两种做法全局复位所有PE的acc在start信号来时统一清零。简单但需要额外的复位布线首周期清零利用第一个有效周期的en信号在状态机里产生一个clear脉冲我选的是第二种因为它不依赖全局复位网络时序更好收敛。具体做法是在控制状态机里当K计数器为0时让PE执行acc a_in * b_in而不是累加从第二个周期开始才累加。3.4 关键代码片段module pe #( parameter DATA_WIDTH 8, parameter ACC_WIDTH 32 )( input wire clk, input wire rst_n, input wire en, input wire clear, input wire signed [DATA_WIDTH-1:0] a_in, input wire signed [DATA_WIDTH-1:0] b_in, output reg signed [DATA_WIDTH-1:0] a_out, output reg signed [DATA_WIDTH-1:0] b_out, output reg signed [ACC_WIDTH-1:0] acc ); wire signed [2*DATA_WIDTH-1:0] mult_result; assign mult_result a_in * b_in; always (posedge clk or negedge rst_n) begin if (!rst_n) begin a_out 0; b_out 0; acc 0; end else if (en) begin a_out a_in; b_out b_in; if (clear) acc {{(ACC_WIDTH-2*DATA_WIDTH){mult_result[2*DATA_WIDTH-1]}}, mult_result}; else acc acc {{(ACC_WIDTH-2*DATA_WIDTH){mult_result[2*DATA_WIDTH-1]}}, mult_result}; end end endmodule这段代码里有个细节mult_result需要符号扩展到ACC_WIDTH位再累加。如果直接加Verilog会按最大位宽对齐但符号位可能处理不对。我显式做了符号扩展确保负数乘法结果正确。4. 数据调度让A和B在正确的时间相遇4.1 斜向注入的数学推导这是脉动阵列最烧脑的部分。假设阵列是N×N我们要算C A×BA和B都是N×N先考虑方阵。PE(i,j)需要计算C[i][j] Σ(k0 to N-1) A[i][k] * B[k][j]在时刻tPE(i,j)应该收到A[i][t]和B[t][j]。也就是说A的第i行第k个元素应该在时刻k到达PE(i,0)然后向右流动在时刻kj到达PE(i,j)B的第k列第j个元素应该在时刻k到达PE(0,j)然后向下流动在时刻ki到达PE(i,j)所以注入策略是A矩阵第i行的数据从左侧第i个端口注入但第i行要比第0行延迟i个周期开始注入。这样A[i][0]在时刻i到达PE(i,0)A[i][k]在时刻ik到达PE(i,0)向右流动j个周期后在时刻ikj到达PE(i,j)。B矩阵第j列的数据从顶部第j个端口注入第j列比第0列延迟j个周期开始注入。B[0][j]在时刻j到达PE(0,j)B[k][j]在时刻kj到达PE(0,j)向下流动i个周期后在时刻kji到达PE(i,j)。两者在PE(i,j)的到达时刻都是ikj完美对齐。4.2 用移位寄存器实现延迟注入在硬件里这个延迟i个周期怎么实现最简单的方法是用移位寄存器链。对于A矩阵的注入第0行直接注入第1行经过一级寄存器延迟第2行经过两级以此类推。B矩阵同理。// A矩阵行延迟链示例简化 reg signed [7:0] a_delay [0:N-1][0:N-1]; integer i, j; always (posedge clk) begin for (i 0; i N; i i 1) begin a_delay[i][0] a_data[i]; for (j 1; j N; j j 1) a_delay[i][j] a_delay[i][j-1]; end end // a_delay[i][i] 就是延迟i个周期后的A第i行数据这段代码综合出来是一堆寄存器N8时是64个8位寄存器资源消耗可以接受。如果N更大可以考虑用BRAM做延迟线但控制会复杂一些。4.3 控制状态机的设计整个加速器需要一个主状态机来协调IDLE等待start信号LOAD从外部存储器读取A和B矩阵数据存入输入缓冲COMPUTE启动注入运行2N-1个周期因为最晚到达的数据需要ikj的最大值iN-1, kN-1, jN-1总共3N-3个周期但实际有效计算是N个周期加上注入和排空总共约3N个周期DRAIN等待最后一个PE完成累加OUTPUT将PE阵列的结果读出到输出缓冲DONE拉高done信号这里有个容易搞错的地方COMPUTE阶段到底要跑多少个周期。我一开始以为是N个周期结果发现最后一个PE的结果要等到第3N-3个周期才稳定。正确的做法是从第一个数据注入开始数到最后一个数据到达最后一个PE并完成累加总共需要(2N-1) (N-1) 3N-2个周期。保险起见我跑了3N个周期。5. 完整Verilog实现与仿真验证5.1 顶层模块结构顶层模块把PE阵列、延迟链、控制状态机、输入输出缓冲都例化在一起module systolic_top #( parameter N 8, parameter DATA_WIDTH 8, parameter ACC_WIDTH 32 )( input wire clk, input wire rst_n, input wire start, input wire signed [DATA_WIDTH-1:0] a_data [0:N-1], input wire signed [DATA_WIDTH-1:0] b_data [0:N-1], output wire signed [ACC_WIDTH-1:0] c_data [0:N-1][0:N-1], output wire done ); // PE阵列例化 // 延迟链 // 控制状态机 // ... endmodule5.2 仿真测试平台验证是重中之重。我写了一个testbench用随机数生成A和B矩阵然后在Verilog里算一遍参考结果跟硬件输出对比。// 参考模型 integer i, j, k; reg signed [ACC_WIDTH-1:0] ref_c [0:N-1][0:N-1]; initial begin for (i 0; i N; i i 1) for (j 0; j N; j j 1) begin ref_c[i][j] 0; for (k 0; k N; k k 1) ref_c[i][j] ref_c[i][j] a_mat[i][k] * b_mat[k][j]; end end仿真时我用了Vivado自带的仿真器跑1000个随机测试用例全部通过后才上板。实操心得仿真通过不代表上板能跑。我遇到过仿真全对、上板结果错位的情况原因是综合工具对数组的推断跟仿真器不一致。解决办法是避免在可综合代码里使用多维数组作为端口改用扁平化的向量或者用generate展开。5.3 综合与实现的关键约束综合时我加了以下约束create_clock -period 10.000 -name sys_clk [get_ports clk] set_input_delay -clock sys_clk 2.000 [get_ports a_data*] set_output_delay -clock sys_clk 2.000 [get_ports c_data*]目标频率100MHz。8×8阵列在Artix-7上综合下来DSP用了64个LUT约2000个FF约1500个时序余量0.8ns能稳定跑100MHz。如果扩展到16×16DSP变成256个时序余量会降到0.2ns左右可能需要降到80MHz。这时候可以考虑流水化PE内部的乘法器在乘法器和累加器之间插一级寄存器代价是增加一个周期的延迟但时序能大幅改善。6. 常见问题与排查技巧实录6.1 结果全错或部分错位症状仿真结果跟参考模型对不上或者只有对角线元素对。排查思路先检查延迟链的深度是否正确。A的第i行应该延迟i个周期B的第j列延迟j个周期。我犯过的错是把行列搞反了。检查PE的clear信号时序。如果clear晚了一个周期第一个乘积会被丢掉。检查符号扩展。用负数测试用例跑一遍如果负数结果错就是符号位没处理好。6.2 时序不收敛症状综合后时序报告显示setup violation。排查思路看关键路径在哪里。通常是PE内部的乘法器到累加器这条路径。解决办法在乘法器输出和累加器之间插寄存器流水化或者降低目标频率。另一个常见问题是复位网络。如果用了全局复位复位信号的扇出很大时序会差。改用局部同步复位。6.3 资源不够用症状综合报错说DSP或LUT不够。排查思路先算一下理论消耗N×N个DSP每个PE约30个LUTN×N×30个LUT。如果DSP不够可以考虑时分复用用N/2×N/2的阵列跑两轮但控制逻辑会复杂很多。如果LUT不够检查延迟链是不是用了太多寄存器。N16时延迟链是256个8位寄存器可以改用BRAM。6.4 常见问题速查表问题现象可能原因解决方法结果全为0clear信号一直有效检查clear脉冲宽度结果偏大累加器位宽不够溢出增加ACC_WIDTH结果偏小乘法结果被截断检查符号扩展只有部分PE有输出使能信号没连到所有PE检查en信号布线上板后结果随机跨时钟域没处理加同步器时序违例关键路径太长流水化或降频独家避坑技巧在PE阵列的每个PE输出上加一个ILA集成逻辑分析仪抓几个关键PE的acc值跟仿真波形对比。这比盲猜快得多。Vivado的ILA用起来很方便但注意别加太多否则资源不够。7. 从8×8到更大规模扩展思路与优化方向8×8跑通之后下一步自然是扩展。但直接放大阵列规模不是唯一的路有几个方向可以考虑。方向一增加阵列规模。16×16、32×32资源线性增长时序压力也线性增长。适合对吞吐量要求极高的场景。方向二时分复用。用8×8的阵列跑16×16的矩阵分4轮计算。资源不变但计算时间变成4倍。适合资源受限但延迟不敏感的场景。方向三权重固定型改造。把B矩阵预加载到PE的本地寄存器A矩阵流过。这样B不需要每轮重新注入适合推理场景中权重固定的情况。改造工作量中等主要是控制逻辑要改。方向四加入量化与反量化。在PE的输入和输出加移位和饱和逻辑支持INT8输入、INT8输出中间用INT32累加。这是实际部署中最常见的配置。我个人建议先把8×8的输出固定型吃透再根据实际需求选一个方向深入。脉动阵列的精髓在于数据流动的节奏控制这个思维模式一旦建立起来后面做卷积加速器、Transformer加速器都会轻松很多。最后分享一个我在调试时的小习惯每次改完代码先跑一个2×2的小矩阵用纸笔算一遍期望结果跟仿真波形逐周期对比。2×2的时序关系最容易看清确认无误后再放大到8×8。这个笨办法帮我省下了大量盲目调试的时间。