第7卷第4期 智能系统学报 Vol.7 No.4 2012年8月 CAAI Transactions on Intelligent Systems Aug.2012 D0I:10.3969/j.issn.1673-4785.201202002 网络出版地址：http://ww.cnki.net/kems/detail/23.1538.TP.20120712.1121.010.html 基于FPGA的全流水双精度浮点矩阵乘法器设计 刘沛华，鲁华祥，龚国良，刘文鹏 (中国科学院半导体研究所神经网络实验宣，北京100083) 摘要：在数字通信、图像处理等应用领域中需要用到大量的矩阵乘法运算，并且它的计算性能是影响系统性能的 关键因素.设计了一个全流水结构的并行双精度浮点矩阵乘法器以提高计算性能，并在Xilinx Virtex-5LXI55现场可 编程门阵列(PGA)上完成了方案的实现.乘法器中处理单元(PE)按阵列形式排列，在一个FPGA芯片上可集成10 个PE单元实现并行计算.为了提高工作频率，PE单元采用流水线结构，并运用C-sow时序重排技术解决了环路流 水线上“数据相关冲突”的问题.仿真结果表明，该乘法器的峰值计算性能可达到5000MFL0PS.此外，对不同维数的 矩阵乘法进行了实验，其结果也证实了该设计达到了较高的计算性能。 关键词：矩阵乘法；现场可编程门阵列(FPGA);环路流水线；C-slow时序重排技术；乘法器设计 中图分类号：TP332.2文献标志码：A文章编号：16734785(2012)04030205 Design of an FPGA-based double-precision floating-point matrix multiplier with pipeline architecture LIU Peihua,LU Huaxiang,GONG Guoliang,LIU Wenpeng (Lab of Artificial Neural Networks,Institute of Semiconductors,Chinese Academy of Science,Beijing 100083,China) Abstract:Many application areas,such as digital communication and image processing,make extensive use of ma- trix multiplication operations,and the computational performance of this operation is critical for the whole system. A parallel double-precision floating-point matrix multiplier with pipeline architecture was designed to improve the computational performance.The design was implemented in a Xilinx Virtex-5 LX155 field programmable gate array (FPGA).Up to 10 processing elements were integrated in a single FPGA device,and they were arranged as an ar- ray to achieve parallel computation.The processing elements employed pipelined architecture to increase the speed, and C-slow retiming was applied to solve the data-related conflicts issues on the loop pipeline.The post-Route sim- ulation results show that the peak performance of the matrix multiplier can achieve 5000 MFLOPS.In addition,the matrix multiplication experiments with different dimensions were carried out,and the results confirm that the design achieved high computational performance. Keywords:matrix multiplication;FPGA;loop pipeline;C-slow retiming;multiplia design 矩阵乘法是数字信号处理领域中的基本操作， 用处理器串行计算来实现，严重制约了计算速度.要 广泛应用于各种电路计算中，例如数字通信领域的 提高矩阵乘法的计算性能，可以通过提升工作频率 DCM变换、快速FFT变换以及图像处理中的3-D变 和算法并行度来实现.现场可编程门阵列(field pro- 换等都用到了大规模的矩阵乘法运算.由于矩阵乘 grammable gate array,FPCA)具有强大的计算性能 法计算复杂性较高（通常为0(n)),其计算性能直 和逻辑分析能力，特别是它具有并发式的硬件结构 接影响到系统的整体性能.然而传统的矩阵乘法多 和出色的浮点计算性能，适合对矩阵乘法进行硬件 加速，是当前的研究热点。 收稿日期：2012-02-06.网络出版日期：201207-12. 目前，采用PGA实现矩阵乘法计算的研究已 基金项目：国家自然科学基金资助项目(61076014)；江苏省高校自然 科学基金资助项目(10KJA510042);先导项目 经取得一些成果.在定点矩阵乘法方面，Amira等在 (XDA06020700). 通信作者：刘沛华.E-mail:pelph123@163.com. FPGA上实现了8位定点的矩阵乘法器，但是该设
第4期 刘沛华，等：基于FPGA的全流水双精度浮点矩阵乘法器设计 ·303· 计所需要的带宽与矩阵规模成比例增加，限制了该 图1中信号分为2类：数据信号和控制信号.其 设计的可扩展性：Jang等设计的矩阵乘法器只需 中a、b、MULm_RESULT、ADD_RESULT、Q_FB和 要固定的带宽，但是所需要的存储单元大小与矩阵 output都是64bit数据信号，其他都是控制信号，具 规模成正比2]，在浮点矩阵乘法方面，Campell等设 体功能描述见表1. 计了一个并行结构矩阵乘法器，该设计中的各个计 表1PGA内部资源使用情况 算单元之间不需要通讯，具有可扩展性，但其所需的 Table 1 Usage of FPGA internal resources 存储空间随矩阵维数的增加而增大，并且计算效率 信号名称 功能描述 不高3]；田翔等设计了一个实时双精度矩阵乘法 a、b 浮点乘法器的2个输人 器，并在FPGA上完成了方案的实现，但是其计算单 元的工作频率不高，限制了计算性能的提升 乘法器的输出结果，作为加法器的一个 MULTI RESULT 输入 本文设计并在PGA上实现了一个计算性能较 高、可扩展性良好的并行双精度浮点矩阵乘法器.为 ADD_RESULT加法器的输出结果 提高工作频率，乘法器中的计算单元采用流水线结 Shift register的输出，作为浮点加法器的 Q_FB 构，并运用C-slow时序重排技术解决了环路流水线 另一个输入 上“数据相关冲突”的问题，提高了计算效率.此外， 与MULTI_RDY是一对握手信号，置高 本设计所需要的带宽和存储单元大小都是固定的， MULTI OP ND 表示乘法器的输入有效 故可扩展性好 1)与MULTI_OP._ND是一对握手信号， 1 矩阵乘法器设计 MULTI RDY 置高表示乘法器的输出有效；2)、控制 Shift register移位 对于矩阵乘法C=A×B,其中A、B和C分别 与ADD_RDY是一对握手信号，置高表 是M×L、L×N和M×N维矩阵，其计算方法如式 ADD_OP_ND 示加法器的输入有效 (1): 1)与ADD_OP_ND是一对握手信号，置 年=含a×g1≤i≤M,1≤j≤N(① ADD RDY 高表示加法器的输出有效；2)控制Shi register移位 式(1)的计算复杂度为2×M×W×L,即0(n3).为 RD EN Sh计register的读使能信号 降低算法复杂度，本文设计了一个包含P个处理单 NEW_OP 元(processing element,PE)的并行双精度浮点矩阵 Shi近register的清零信号 乘法器，其中PE采用流水线结构，并运用C-slow时 PE工作时，MULTI_OP_ND置高表示乘法器的 序重排技术解决环路流水线上“数据相关冲突”的 输入（α，b)有效，此时启动乘法器，当浮点乘法器计 问题，提高了计算效率.设计中所有操作数均为符合 算完毕时，MULTI._RDY输出高电平表示MULTI EEE754标准的64bit双精度浮点数 RESULT有效；当ADD_OP_ND变为高电平时启动 1.1处理单元设计 浮点加法器，当浮点加法器计算完毕时，ADD_RDY PE是构成浮点矩阵乘法器的基本单元.每个PE 输出高电平表示乘加过程结束.此外，MULTI_RDY 包含一个浮点乘法器、一个浮点加法器和用于存储计 和ADD_RDY还作为Shift register的控制信号，控制 算结果的存储单元(Shift register),其结构如图1所示. Shift register移位.当一组元素计算完毕之后，将输 MULTT OP ND 人信号RD_EN置为高电平，读取Shi进register中存 PE 储的结算结果；当需要计算新的元素时，将NEW MULTI RDY MULTI RESULT ADD OP ND OP置为高电平，Shift register的寄存器全部清零，便 Q FB 可以进行新一轮的乘加运算 ADD RDY TADD RESUL 1.2环路流水线设计 output 1..n- 为了提高计算吞吐率，整个PE采用流水线结 RD E Shift register 构.与一般的结构相比，流水线结构能达到更高的时 NEW OP 钟频率，但是输出结果与输入之间会有时钟延迟，延 图1PE单元结构 迟的时钟周期数等于流水线的级数.流水线的级数 Fig.1 Structure of PE (Latency)越高，乘法器与加法器的工作频率就越高
·304 智能系统学报 第7卷 对于无反馈回路的流水线结构来说，输出结果 个任务的模式，它将多个计算任务交叉编排在一起， 相对输入之间的延迟不会影响整个系统的顺序执 按次序打入流水线.要实现C-slow模式，需要对输 行.但是当流水线结构中存在反馈回路时，若不妥善 入数据进行预处理，将输入数据重新编排，使得a、b 解决延迟问题，流水线上就会出现“数据相关冲 端口的数据分别满足式(2)、(3)： 突”.以PE的数据通道为例，图1中加法器端口2 dinA[k·(v+1)·T+sT+T]=ak,(2) 的输入数据来自于上一次乘积累加操作的结果，这 dinB[k·(v+1)·T+sT+T,]=bk,(3) 便构成了一个反馈回路，只有保证MULTI_RESULT 式中：T。为计算的起始时间，d_inA[t]、d_inB[t]分 到达加法器端口1的时间与上次乘累加的结果(Q_ 别表示t时刻a、b端口的输入数据，1≤k≤L,0≤s≤ FB)到达端口2的时间一致，才能确保PE的正确运 u.由式(2)、(3)可知，上述操作的目的是把求解cg, 行.否则，必然导致流水线时序紊乱，无法完成给定 C+1w,…,C+这v+1个计算任务交叉编排在一条 的计算任务，这就是所谓的“数据相关冲突”，下面 流水线上执行.同时，上述操作能保证对于c,的计算 通过剖析图2来阐述这个问题， 任务而言，前后2组输入数据保持v+1个时钟周期 T+T- 的延迟，而且不同c的输入数据不发生重叠.所以， a 64b0a6450■a,64b0 b 64b0b64b0b64b0 这样一条流水线上能完成：+1计算任务的交叉执 MULTI OP ND 行，即一个PE单元花费L(T2+T)时间能完成c, MULTI RDY ADD RDY T C+1,…,c+的v+1个元素的求解过程.从而，图1 MULTT RESULT 164'b0MX64b0MMX64'b0 所示的shi讥register需要v+1个寄存器来存储计算 ADD RESULT 6404Y64b0A2X64b0 结果，即n=v+1. Q FB 64b0O64b0☐O6460 1,3并行结构矩阵乘法器设计 CLK 几ΠΠΠΠΠΠ几几ΠΠΠ几 根据矩阵乘法的简单并行算法，在FPGA芯片 图2PE的时序波形 上实现P个PE单元，这些PE单元按照图3所示的 Fig.2 Timing diagram of PE 1×P阵列形式排列.PE单元之间不存在信息交互， 如图2所示，MULTIRDY比MULTI_OP_ND延 它们独立地完成各自的计算任务.由式(2)和(3)可 迟T1;ADD_RDY比ADD_OP_ND(即图2中的 知，每个PE单元进行计算时要用到A的n行和B MULTI_RDY)延迟T2.设CLK的周期为T,浮点乘 的某1列数据，整个PE阵列一次计算需要用到A 法器和浮点加法器的Latency值分别为u和v,则 的n行和B的P列数据.将输入矩阵A、B分别按行 T1=uT,T2=vT.设某个计算元素c的前后2组输入 和按列进行分块：A=[ATA…AW]T,B=[B1B2… 数据分别是a1、b,和a2、b2.开始计算时Shift register Bw],其中A:表示A的第i行，B表示B的第j列 的寄存器被全部清零，A1=a1b1+0,A1经过一个寄 将A的n行和B的P列作为图3所示系统的输人， 存器延迟得到Q1;同样，M2=a2b2.由图2可以看 图中预处理模块1和预处理模块2的功能分别对 出，只有当a2、b2和a1、b1之间保持v+1个时钟周 A:和B进行处理，使得它们分别满足式(2)和(3) 期的延迟时(T2+T=(v+1)T),才能保证M2和Q1 中d_inA[t]、d_inB[t]的要求并将其送入各个PE 同时到达浮点加法器2个输入端口，进而得到正确 单元的a、b端口进行计算. 的累加结果：A2=Q1+M2·否则，就会导致M2和 B B Q1到达浮点加法器输入端口的时间不一致，发生 预处理 预处理 预处理 “数据相关冲突”，无法得到正确的计算结果 处 模块2 模块2 模块2 上述分析表明，只有保证计算c:的任意2组输 人(ak、bg和a.k+1、bk+1)之间保持v+1个时钟周期 块 的延迟，PE才能正常工作，然而，这样一条流水线上 E PE 只有一个计算任务时，时序会出现大量的空闲时间 图3并行矩阵乘法器结构 T2 (如图2所示)，空闲的比例为T+为了提高时间 Fig.3 Timing diagram of PE 通过这种PE的阵列结构，可以完成任意维数 利用率，本文采用C-slow模式5]来解决上述“数据 的矩阵乘法运算.假设A和B分别为M×L、L×N 相关冲突”的问题.C-slow是一种解决环路流水线问 维矩阵，对于任意的M、N和L值，可以通过下述算 题的时序重排技术，不同于传统流水线上只执行一
第4期 刘沛华，等：基于FPGA的全流水双精度浮点矩阵乘法器设计 ·305· 法计算C=A×B的结果： 线后仿真的结果，该矩阵乘法器在未做优化的情况 for i=0 to M/n-1 loop 下工作频率能达到250MHz.由此可知该矩阵乘法 for j=0 to N/P-1 loop 器的峰值计算性能可达到5OO0 MFLOPS. all PE,(x=1,2,,P)do in parallel 表2FPGA内部资源使用情况 for r =1 to n loop Table 2 Usage of FPGA internal resources C+r办+g=0 内部资源 资源使用量/个占总资源比例/% end loop Slice Registers 16137 16 for =1 to n loop Slice LUTs 14590 复 for r=1 to n loop BRAM 12 6 Cintrjp+l =aitrk Xbkjpt Cintrp BUFG 2 6 end loop DSP48E 90 end loop 70 end loop 2.2平均计算性能 end loop 并行矩阵乘法器的平均计算性能可以通过计算 从以上算法可以看出，使用并行矩阵乘法器进 2个矩阵相乘所需的总时间来求得，如式(5)所示. 行计算时，循环的次数是传统串行算法的1/P,即计 PERF F/t. (5) 算复杂度降低为O(n/P).同时由于该并行矩阵乘 式中：PRF表示并行矩阵乘法器的平均计算性能， 法器中的各个PE单元是相互独立的，因此可以方 F表示2个矩阵相乘总共需要完成的双精度浮点操 便地扩展到多片PGA上实现并行计算， 作次数，t为计算时间。 本文分别以2个128bit×128bit的矩阵相乘和 2矩阵乘法器性能分析 2个256bit×256bit的矩阵相乘的实例来分析该设 下面以在FPGA上实现的并行矩阵乘法器来对 计的平均计算性能，如表3所示。 上述设计的性能进行分析.本文选用Xilinx Virtex-5 表3平均计算性能对比 Table 3 LX155芯片实现该设计.PE中的浮点乘法器和浮点 Comparison of computation performance 加法器使用Xilinx公司提供的floating-point IP核生 性能 本文设计 设计A4 成.通过对运行速度及该器件中DSP48E单元、CLB 单元等资源进行综合考虑，对并行矩阵乘法器进行 矩阵维数 128×128256×256 100×100200×200 如下设置：1)P核生成浮点乘法器时，DSP48E的使 浮点数操 作次数 22 22 2×1016×109 用等级设置为Medium Usage(即单个浮点乘法器使 用9个DSP48E单元)，Latency的值设定为15；2)P 计算时间/ 986.37672.8 1351 10759 核生成浮点加法器时，DSP48E的使用等级设置为 (t·us1) 所需带宽/ No Usage,Latency的值设定为9；3)设定矩阵乘法器 4 4 2.4 2.4 (GB·s1) 中PE单元的个数P=10.本设计中所使用的FPGA PERF/ 4252 4373 1481 1490 开发环境和仿真环境分别为Xilinx ISE Design Suite MFLOPS 13.1 Mentor Graphics Modelsim SE 6.5a. PERF -/% PERFpo 85.04 87.46 49.4 49.7 2.1峰值计算性能 理想情况下，每个PE单元在一个时钟周期内 PERFp/ 5000 3000 可以完成1次双精度浮点乘法操作和1次双精度浮 MFLOPS 点加法操作，因此整个矩阵乘法器的计算性能可计 由表2、3可以看出，本文设计的并行矩阵乘法 算为 器的峰值计算性能可达到5000 MFLOPS,平均计算 PERFnk =2 x P x f. (4) 性能可以保证在峰值计算性能的85%左右.而田翔 式中：PERF表示矩阵乘法器的峰值计算性能（每 等的设计A未采用流水线结构，工作频率只有60 秒百万次浮点操作)，P为矩阵乘法器中PE单元的 MHz,即使在一片FPGA上集成了25个PE单元，它 个数，f为矩阵乘法器工作的时钟频率。 的峰值计算性能只能达到30O0 MFLOPS.而且， PGA内部资源的使用情况见表2.根据布局布 设计A的平均浮点计算性能只能保持在蜂值计算
·306 智能系统学报 第7卷 性能的50%左右.由此可见，本文设计在计算性能 circuitry by retiming[C]//Proceedings of the 3rd Caltech 上有大幅度提高. Conference On VLSI.Rockville,Maryland,1983:87-116. 在矩阵乘法计算中，若FPGA的I/O带宽小于 [6]SU Ming,ZHOU Lili.Maximizing the throughput-area effi- 一定值，并行矩阵乘法器中的PE单元就会出现等 ciency of fully-parallel low-density parity-check decoding 待状态，此时，带宽便成为制约计算性能的因素.当 with C-slow retiming and asynchronous deep pipelining [C]//The 25th International Conference on Computer De- /0带宽达到或者高于这个值后，每个PE单元的计 sign.Washington,DC,USA,2007:93-100. 算性能则成为制约并行矩阵乘法器计算性能的主要 [7]肖宇，王建业，张伟.基于P核的数选式浮点矩阵相乘 因素.在本文的设计中，PE的计算性能主要由工作 设计[J].集成电路应用，2011,37(6)：52-55. 频率决定，在工作频率为250MHz的情况下，只要 XIAO Yu,WANG Jianye,ZHANG Wei.Floating-point ma- I/O带宽达到4GB/8,便不会对整个系统的计算性 trix multiplication design based on IP core[].Application 能产生影响 of Integrated Circuits,2011,37(6):52-55. [8]许芳，席毅，陈虹.基于FPGA Nios--Ⅱ的矩阵运算硬件 3结束语 加速器设计[J].电子测量与仪器学报，2011,25(4)： 本文设计了一个全流水结构的并行双精度浮点 376-383. 矩阵乘法器，并在Xilinx xc5vxl55FPGA上实现了 XU Fang,XI Yi,CHEN Hong.Design and implementation of matrix hardware acceleration based on FPGA/Nios-II[J], 该方案.矩阵乘法器内部的PE单元采用流水线结 Joumnal of Electronic Measurement and Instrument,2011, 构，并运用C-slow时序重排技术解决了环路流水线 25(4):376-383. 中“数据相关冲突”的问题，提高了计算效率.实验 [9]黎铁军，李秋亮，徐炜遐.一种128位高性能全流水浮 结果表明，对于不同维数的矩阵乘法，本设计都有较 点乘加部件[J].国防科技大学学报，2010,32(2)：56 高的计算性能.同时，本文设计的并行矩阵乘法器结 60. 构，其内部的各个PE单元相互独立，因而具有很好 LI Tiejun,LI Qiuliang,XU Weixia.A high performance 的可扩展性，在后续的研究工作中，需要提出更为合 pipeline architecture of 128 bit floating-point fused multiply 理的并行结构，通过多片PGA的并行计算来进一 add unit[J].Journal of National University of Defense Tech- 步提高矩阵乘法器的计算性能， nolog7,2010,32(2):56-60. 作者简介： 参考文献： 刘沛华，女，1985年生，硕士研究 生，主要研究方向为电路与系统、神经 [1]AMIRA A,BENSAALI F.An FPGA based parameterizable 网络。 system for matrix product implementation[C]//IEEE Work- shop on Signal Processing Systems (SPIS'02).San Diego, 2002:75-79. [2]JANG J,CHOI S,PRASANNA V KK.Area and time effi- cient implementations of matrix multiplication on FPGAs 鲁华祥，男，1965年生，研究员，主 [C]//2002 IEEE International Conference on Field Pro- 要研究方向为智能信息处理、神经网络 grammable Technology.Seoul,Korea,2002:93-100. 技术及其应用.近年来，作为项目负责 [3]CAMPBELL S J,KHATRI S P.Resource and delay effi- 人或骨干研究人员完成国家重大科技 cient matrix multiplication using newer FPGA devices[C]/ 攻关项目3项、国家“863”计划项目3 Proceedings of the 16th ACM Great Lakes Symposium on 项、国家自然科学基金重点项目3项， VLSI.Philadelphia,USA,2006:308-311. 发表学术论文50余篇。 [4]田翔，周凡.基于FPGA的实时双精度浮点矩阵乘法器 设计[J].浙江大学学报：工学版，2008,42(9)：1611- 龚国良，男，1982年生，博士研究 1615. 生，主要研究方向为优化算法、神经网 TIAN Xiang,ZHOU Fan.Design of field programmable gate 络、模式识别等， array based real-time double precision floating-point matrix multiplier[J].Joural of Zhejiang University:Engineering Science,2008,42(9):1611-1615. 5]LEISERSON C,ROSE F,SAXE J.Optimizing synchronous
