简介SystolicArrayDemo是一份Python示例代码展示脉动阵列Systolic Array的实际运行逻辑。项目面向对高性能计算、硬件加速器以及TPU设计原理感兴趣的学习者通过软件模拟呈现二维计算单元在矩阵乘法中的流水线式数据流动。资源压缩包仅3KB共2个文件一个.py程序负责模拟单元间数据传递、乘法累加与并行调度一份md说明文档介绍设计原理与运行说明结构简洁适合快速阅读和二次修改。目前已有442人学习浏览。演示内容覆盖脉动阵列的矩阵乘法算法、单元间数据流管理、流水线并行机制等关键点并延伸到多线程/多进程模拟并行、硬件调度规则以及TPU这类专用集成电路的加速背景。读者既能借此理解计算密集型任务的加速思路也可将其作为软硬件协同设计的起点通过修改参数观察不同数据流对吞吐量和延迟的影响为后续深入FPGA或ASIC实现打下基础。1. 脉动阵列为什么它的秘密不是算得快而是搬得少第一次在书里看到脉动阵列那个菱形的格子图时很容易把它当成一种“硬件玄学”数据从一个格子流到另一个格子每过一个格子做一次乘加看起来只是把很多乘法器连起来而已。真正反直觉的点在另一个地方——它解决的瓶颈不是算不出来而是喂不进去。一个最内层循环每次都从内存取两个数的矩阵乘法访存量是计算量的两倍而处理器的访存能力远比算力增长得慢。脉动阵列把“从内存搬数据”改成了“在芯片内传数据”同一个数被一排计算单元接力复用从外部只读一次。主流AI加速芯片的矩阵乘单元大多沿用这条思路。打开这份用Python逐时钟周期模拟的 Demo你会看到每个计算单元每一拍手里拿着哪个数、累加器变了多少——黑匣子第一次有了透明窗口。2. 从内存墙到数据复用脉动阵列到底解决了什么问题2.1 内存墙为什么一次乘累加读两个数会要命矩阵乘法的最小操作是乘累加MAC算 C[i][j] 时要对 k 从 0 到 K-1 执行acc A[i][k] * B[k][j]。一次 MAC 需要两个操作数一个来自 A一个来自 B。如果这两个数都在片外内存里那这次计算就欠了内存一笔“搬运债”。把账算清楚M×N×K 次 MAC对应 2×M×N×K 次数据读取。取 MNK1024总访存次数大约 21 亿次。处理器每个时钟周期能做的乘加次数远超它能从内存取数的次数于是出现一个尴尬局面乘法器在空转等数据内存总线在满负荷跑。这就是内存墙——算力提升的速度远高于访存带宽提升的速度再快的乘法器也会被喂数据的速度卡住喉咙。我一般这么理解内存墙它不是一条物理的墙而是一个比值——每做一次算术操作平均要从外面拿几个数进来。普通三重循环的这个比值是 2意味着计算单元每干一份活内存要搬两份数据。搬数据消耗的能量比做乘法高一个数量级所以这个比值直接决定了芯片的上限。2.2 数据复用让同一个数多干几份活矩阵乘法里的数据复用规律非常整齐A[i][k] 要被用来计算 C[i][0]、C[i][1]……一直到 C[i][N-1]一共 N 个输出B[k][j] 要被用来计算 C[0][j]、C[1][j]……一直到 C[M-1][j]一共 M 个输出。也就是说每个输入数据天生就该被复用几十上百次而不是读一次用一次。普通处理器靠缓存来抓这个规律A[i][k] 被读进缓存后后续 N 次使用都可能命中。但缓存命中靠的是“局部性碰运气”——数据在内存里的排布、分块策略、多核竞争都会影响命中率。一旦 miss就要付出几十上百个周期的代价回到内存去取。脉动阵列换了一个思路把“靠缓存碰运气”变成“靠结构保证”。它把 B 的每一列预先“驻留”在对应列的计算单元里A 的行数据从左侧进入后沿着一整行计算单元向右流动。A[i][k] 每经过一个计算单元就和本地权重做一次乘加然后传给下一个邻居。同一份 A 的数据在芯片内部被 N 个计算单元依次使用全程不回内存。复用不再是缓存命中的概率事件而是物理结构上的必然。这种设计把“数据复用关系”直接映射成了“芯片上的邻居关系”。你在纸上把矩阵乘法的依赖图画出来那些连线就是脉动阵列的布线。理解了这一层后面看时序推演就会顺很多。2.3 数量级对比访存/计算比从 2 降到约 0.003做一个可验证的估算。普通实现忽略缓存假设每次 MAC 都从内存取数的访存/计算比是 2 words/MAC。脉动阵列的一版典型账本A 的 M×K 个元素从片外读一次进入阵列后被 N 个 PE 复用B 的 K×N 个权重在计算开始前加载进 PE 本地寄存器之后不再访存输出的 M×N 个结果写回。总访存约 MK KN MN总计算量是 MNK于是实现方式每 MAC 访存wordsMNK1024 时总访存普通三重循环无缓存2约 21 亿脉动阵列权重驻留1/N 1/M 1/K约 314 万当阵列规模到 1024 时脉动阵列把访存/计算比压低了近 700 倍。这 700 倍不是靠提高时钟频率换来的而是靠把数据搬运转化为数据流动换来的。后面第 4 章的 Python 仿真会直接验证这个过程——你能在打印出的网格状态里看到同一个输入值确实是在一拍一拍地向右走而不是反复被重新读取。3. 时钟、权重与驻留逐拍拆解脉动阵列的流水3.1 三个角色权重驻留、输入流动、累加器本地化先把角色分配清楚。这个演示用的是“权重驻留weight stationary”方案它是各种矩阵加速单元里最常见的一种。阵列里有 M×N 个处理单元PE每个 PE 承担一个输出元素 C[i][j] 的全部累加工作。三个角色各管一件事。权重驻留B 矩阵的第 j 列整体装载进第 j 列的每个 PE 的本地寄存器组计算开始前一次性完成计算过程中不流动、不重读。输入流动A 矩阵的第 i 行从左边界逐拍注入每拍向右移动一格沿路在每个 PE 里和本地权重相乘。累加器本地化C[i][j] 的部分和一直待在 PE 内部的累加寄存器里就地累加不参与流动等所有 K 个输入都走完直接从 PE 里读走。为什么选这个方案而不是“部分和流动”方案因为教学场景下它把数据复用讲得最直接——你能亲眼看到同一个输入值沿着一排 PE 走完全程而不是被不同的数据源反复喂进来。部分和流动的变体会把输出从底部排出去多一个排空drain阶段时序更绕留到以后进阶再看。3.2 逐拍推演3×3 矩阵乘法在 9 个 PE 里发生了什么用一组具体数来推演。设 A 是 3×3 的 [[1,2,3],[4,5,6],[7,8,9]]B 是 [[1,0,2],[0,1,0],[2,0,1]]要算 C A×B。直接手算 C 的第一行是 [7, 2, 5]7 1×1 2×0 3×22 1×0 2×1 3×05 1×2 2×0 3×1。只看第 0 行 PE它们分别负责算 C[0][0]、C[0][1]、C[0][2]。每个 PE 本地驻留的权重是 B 的对应列第 0 列 PE 持有 [1, 0, 2]第 1 列 PE 持有 [0, 1, 0]第 2 列 PE 持有 [2, 0, 1]。A 的第 0 行 [1, 2, 3] 从左边逐拍注入。周期 tPE(0,0) 收到PE(0,1) 收到PE(0,2) 收到第 0 行累加器01×1 → 1空空[1, 0, 0]12×0 → 11×0 → 0空[1, 0, 0]23×2 → 72×1 → 21×2 → 2[7, 2, 2]3空3×0 → 22×0 → 2[7, 2, 2]4空空3×1 → 5[7, 2, 5]这个表是整个机制的核心。注意看数据到达 PE 的时间输入 1 在 t0 到达第 0 列t1 到达第 1 列t2 到达第 2 列。输入 3 在 t2 才刚进入第 0 列。所以同一时刻不同列的 PE 其实在处理 A 的不同元素——这正是“脉动”二字的来源数据像队列一样一拍一拍地挪而不是所有 PE 同时面对同一个新数据。另一个容易忽略的细节t0 到 t2 是“填充”阶段右侧的 PE 在等数据流入t4 时所有输入已经走完但最后一个 PE 还在做最后一次累加。填充和排空不是错误而是流水线的正常组成部分。后面章节的参数和避坑都会反复提到它。3.3 关键参数阵列尺寸、累加深度、总周期数怎么定这个演示里几个参数的含义和影响值得单独列出来。理解它们才能把代码改成任意尺寸的矩阵乘法而不翻车。参数含义本例取值对流水的影响M输出行数等于阵列行数3行与行完全并行M 变化不改变总周期数N输出列数等于阵列列数3决定数据要流经多少个 PE直接加长填充阶段K累加深度即内积维度3决定每个 PE 要连续做多少次 MACtotal_cycles从第一个输入注入到最后一个结果完成K N - 1 5少了算不完多了只是空转为什么总周期数是 K N - 1 而不是更自然的 K N推导一下第 0 个输入在 t0 进入第 0 列走到第 N-1 列需要 N-1 拍所以它到达最后一列是 tN-1。第 K-1 个输入在 tK-1 进入第 0 列到达最后一列就是 tK-1N-1KN-2。周期编号从 0 到 KN-2一共 KN-1 个周期。另一个值得注意的点total_cycles 与 M 无关因为 M 行是平行的各行同步推进、互不干扰。这也是为什么这个演示只打印第 0 行就能代表整个阵列的行为。4. Python 逐周期仿真让每个 PE 的状态可见4.1 用二维数组模拟 PE 网格与寄存器代码层面不需要真的定义“硬件单元”用三个二维数组就能把状态装下。第一个是权重表pe_weights[j]存 B 的第 j 列第二个是累加器acc[i][j]对应 C[i][j]第三个是输入寄存器reg[i][j]表示 PE(i,j) 当前某一拍手里握着的 A 的元素。这三个数组就是全部状态。import numpy as np M, N, K 3, 3, 3 A np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtypeint) B np.array([[1, 0, 2], [0, 1, 0], [2, 0, 1]], dtypeint) # 权重驻留把 B 的第 j 列整体预装载给第 j 列的每个 PE # pe_weights[j][k] 对应原矩阵 B[k][j] pe_weights [[B[k][j] for k in range(K)] for j in range(N)] # 累加器C[i][j] 就地累加 acc [[0] * N for _ in range(M)] # 输入寄存器reg[i][j] 是 PE(i, j) 当前持有的 A 的元素 reg [[0] * N for _ in range(M)] total_cycles K N - 1这三个数组和硬件的对应关系是pe_weights对应 PE 的本地寄存器组acc对应累加寄存器reg对应输入锁存器。程序的每一轮for t循环就是硬件里的一个时钟周期。4.2 先注入后计算同一个周期内数据不能穿透两级 PE时钟循环是最容易写错的地方核心难点是“在同一拍里数据不能同时跨越多个 PE”。硬件里相邻 PE 之间是寄存器锁存数据只能在时钟边沿更新一次软件里如果用朴素的双层循环从左到右覆盖右边的 PE 会读到左边这一拍刚更新的值等于一个周期内数据穿越了好几个 PE结果完全错位。for t in range(total_cycles): # 阶段1数据注入与右移模拟时钟边沿的寄存器更新 for i in range(M): feed A[i][t] if t K else 0 # 片外读入新元素K 之后补零 for j in reversed(range(N)): # 从右往左更新避免覆盖本拍旧值 if j 0: reg[i][j] feed else: reg[i][j] reg[i][j - 1] # 阶段2每个 PE 用当前寄存器里的输入做一次乘累加 for i in range(M): for j in range(N): k t - j # 到达 PE(i, j) 的输入序号 if 0 k K: acc[i][j] reg[i][j] * pe_weights[j][k] # 直接做矩阵乘法验证结果 C A B assert acc C.tolist(), 结果不一致 print(acc)关键在reversed(range(N))这一步更新顺序从右往左第 j 列读到的reg[i][j - 1]还是上一拍的值不会被本拍刚写入左邻的新值污染。这模拟了硬件里“所有 PE 在同一时钟沿同时锁存”的语义——要么一起更新要么都保持旧值绝不能出现左边刚算完右边马上就用的情况。第二层关键在k t - j这一行。它表达的本质是一个输入从第 0 列出发向右每走一列就多花一拍所以到达第 j 列时它对应的 A 元素序号是“全局时间减去走过的列数”。if 0 k K的越界判断处理的是填充阶段和排空阶段——左侧还没有数据进来k 0或者 K 个输入已经全部走完k K。4.3 打开黑匣子打印每一拍的网格状态结果正确不等于理解正确。把每个周期的reg和acc打出来才能看到“脉动”到底是什么样子。在循环里加一段调试输出if t 5: print(f--- t{t} ---) print(reg row0:, reg[0]) print(acc row0:, acc[0])运行后第 0 行的轨迹如下--- t0 --- reg row0: [1, 0, 0] acc row0: [1, 0, 0] --- t1 --- reg row0: [2, 1, 0] acc row0: [1, 0, 0] --- t2 --- reg row0: [3, 2, 1] acc row0: [7, 2, 2] --- t3 --- reg row0: [0, 3, 2] acc row0: [7, 2, 2] --- t4 --- reg row0: [0, 0, 3] acc row0: [7, 2, 5]看reg row0的变化方式t0 时 1 在第 0 列t1 时 1 移到第 1 列并且 2 进入第 0 列t2 时 1 到第 2 列、2 在第 1 列、3 在第 0 列。整个序列像一条斜线向右推进这就是“波前”的形状。硬件里实际发生的事就是这条斜线在二维网格上扫过所有 PE。acc row0的变化也值得盯一下t2 时第 0 列已经算完7 到齐但第 2 列还只有第一个输入贡献的 2t4 时第 2 列才补上最后一个 3×1变成 5。这说明不同 PE 的完成时间天然错开不是所有结果都在同一个周期ready。这正是后面要讨论的流水线填充/排空问题。这份代码的验证方式依赖assert——每次跑完都强制和直接矩阵乘法比对。改参数时只要 M、N、K 保持一致pe_weights的转置布局不出错断言就能兜底。5. 仿真的坑与排查时序、边界和那些“看起来对但错了”的结果5.1 数据在一个时钟周期内穿透多个 PE现象打印reg时发现t1 时右侧的 PE 已经拿到了本来要 t2 才该到达的数据整个矩阵的结果列全部错位而且和直接矩阵乘法比对不上。原因更新数组时从左往右遍历reg[i][j] reg[i][j-1]读到的左边值已经被本周期改过了数据在同一拍里多跳了一格。这在硬件里不可能发生——寄存器锁存保证所有 PE 在同一时钟沿拿到的都是上一拍的值。解决更新方向改成reversed(range(N))从右往左覆盖或者干脆用两个数组做双缓冲上一拍的值算完本拍再更新。核心原则只有一条每个周期内PE 只能“看”到自己寄存器里原有的数据新注入的数据要等下一拍才能被邻居使用。用断言比对结果是最快的兜底手段。5.2 周期数算错结果总是缺最后一列现象把total_cycles设成K N - 2跑完发现最后一列的结果全缺项有时改回K N打印时最后几拍全是空转看起来又像代码有问题。原因对总周期数的推导没搞清楚。最后一个输入第 K-1 个在 tK-1 进入第 0 列到达第 N-1 列还要走 N-1 拍最终完成时刻是 tKN-2。周期编号从 0 开始所以需要的周期数就是 KN-1。少一个周期最后一列 PE 的最后一次累加来不及做多一个周期只是浪费几拍空转倒不会算错。解决别背公式直接按“最后一个输入到达最后一列要走多少拍”现场推导一遍。改参数时先对小的 M、N、K 跑断言确认边界无误再上大矩阵。这个坑基本是每个写仿真的人都要踩一次的。5.3 累加器“看起来慢半拍”其实是流水还在填充现象在 t2 时打印acc发现右侧 PE 的累加器还是 0以为权重装载错了或者代码漏了步骤。原因这不是 bug是流水线还没填满。第 2 列 PE 的第一个输入要到 t2 才到达它天然比第 0 列晚两拍开始工作。同理排空阶段左侧 PE 已经算完右侧 PE 还在累加。上一节打印的轨迹里acc row0在 t3 时第 0 列已经是 7 而第 2 列还是 2就是这个现象。解决判断“某 PE 是否完成”不能只看某个固定时刻要看它的累加器是否已经累计了 K 次。更稳妥的做法是在循环结束后统一比对AB而不是在中间某个周期就急着下结论。把填充和排空当成流程的一部分来看而不是当成错误信号。5.4 把行为级仿真当成周期精确模型现象看完代码后觉得硬件也是“先移位、再乘加”两个串行阶段讨论时序时用软件的执行顺序去推硬件的周期。原因这份 Python 仿真是行为级加周期级混合模型——它保证了“每拍数据只移一格”的周期语义但把计算和移位写成了两个 for 循环这是为了代码可读性不代表硬件里真的分两步走。真实电路里组合逻辑在同一个周期内完成乘法加法数据在时钟边沿才更新寄存器两个事件在物理上是重叠的。解决用这份仿真理解数据流和复用关系即可不要拿它推导能量消耗、关键路径或线延迟。想学微架构级细节得去看带门级时序的模型软件仿真在“每周期状态”这一层是对的再往下管到纳秒级就超出它的职责范围了。这个边界意识能省掉不少后续学习里的困惑。6. 设计取舍为什么脉动阵列没有一统天下6.1 面积换带宽利用率没有免费午餐脉动阵列把访存降了一个数量级代价是把大量面积花在了本地寄存器和互联线上。每个 PE 要为 K 个权重准备存储整个阵列存下了 K×M×N 个权重值虽然逻辑上可以广播共享物理上仍要付出可观的寄存器成本。更麻烦的是它对数据供给的形状极其敏感矩阵尺寸和阵列尺寸不匹配时要么填充零元素浪费 PE要么需要额外的数据重排逻辑。软件侧要把数据切成阵列大小合适的块再按节拍精确喂入这比普通的 SIMD 向量化难调度得多。所以脉动阵列不是万能解它是在“片上面积充裕、片外带宽紧张”这个特定约束下做的最优权衡。你在自己的项目里要不要跟进这条思路取决于你的瓶颈到底是乘法器不够用还是内存带宽确实喂不动计算单元。6.2 我的调试验证习惯做这个 Demo 的那段时间有个习惯救了我好几次无论改了什么参数哪怕只是换了个矩阵数值我都先跑一遍断言再把前三个周期的reg状态打印出来人工核一遍。断言管正确性打印管直觉两个都过才敢说这个改动没问题。后来做其它时序类仿真我也强制自己走这套流程——先看波前有没有正常推进再看结果对不对。不是信不过代码是这种流水结构一旦错位错误会被扩散到后面所有周期靠肉眼盯一行输出根本盯不出来。希望这个习惯也能帮到你少踩几个坑。本文还有配套的精品资源点击获取