劉沛華，魯華祥，龔國良，劉文鵬

(中國科學(xué)院半導(dǎo)體研究所神經(jīng)網(wǎng)絡(luò)實驗室，北京100083)

矩陣乘法是數(shù)字信號處理領(lǐng)域中的基本操作，廣泛應(yīng)用于各種電路計算中，例如數(shù)字通信領(lǐng)域的DCM變換、快速FFT變換以及圖像處理中的3-D變換等都用到了大規(guī)模的矩陣乘法運算.由于矩陣乘法計算復(fù)雜性較高(通常為O(n3))，其計算性能直接影響到系統(tǒng)的整體性能.然而傳統(tǒng)的矩陣乘法多用處理器串行計算來實現(xiàn)，嚴(yán)重制約了計算速度.要提高矩陣乘法的計算性能，可以通過提升工作頻率和算法并行度來實現(xiàn).現(xiàn)場可編程門陣列(field programmable gate array，F(xiàn)PGA)具有強(qiáng)大的計算性能和邏輯分析能力，特別是它具有并發(fā)式的硬件結(jié)構(gòu)和出色的浮點計算性能，適合對矩陣乘法進(jìn)行硬件加速，是當(dāng)前的研究熱點.

目前，采用FPGA實現(xiàn)矩陣乘法計算的研究已經(jīng)取得一些成果.在定點矩陣乘法方面，Amira等在FPGA上實現(xiàn)了8位定點的矩陣乘法器，但是該設(shè)計所需要的帶寬與矩陣規(guī)模成比例增加，限制了該設(shè)計的可擴(kuò)展性［1］;Jang等設(shè)計的矩陣乘法器只需要固定的帶寬，但是所需要的存儲單元大小與矩陣規(guī)模成正比［2］.在浮點矩陣乘法方面，Campell等設(shè)計了一個并行結(jié)構(gòu)矩陣乘法器，該設(shè)計中的各個計算單元之間不需要通訊，具有可擴(kuò)展性，但其所需的存儲空間隨矩陣維數(shù)的增加而增大，并且計算效率不高［3］;田翔等設(shè)計了一個實時雙精度矩陣乘法器，并在FPGA上完成了方案的實現(xiàn)，但是其計算單元的工作頻率不高，限制了計算性能的提升［4］.

本文設(shè)計并在FPGA上實現(xiàn)了一個計算性能較高、可擴(kuò)展性良好的并行雙精度浮點矩陣乘法器.為提高工作頻率，乘法器中的計算單元采用流水線結(jié)構(gòu)，并運用C-slow時序重排技術(shù)解決了環(huán)路流水線上“數(shù)據(jù)相關(guān)沖突”的問題，提高了計算效率.此外，本設(shè)計所需要的帶寬和存儲單元大小都是固定的，故可擴(kuò)展性好.

1 矩陣乘法器設(shè)計

式(1)的計算復(fù)雜度為2×M×N×L，即O(n3).為降低算法復(fù)雜度，本文設(shè)計了一個包含P個處理單元(processing element，PE)的并行雙精度浮點矩陣乘法器，其中PE采用流水線結(jié)構(gòu)，并運用C-slow時序重排技術(shù)解決環(huán)路流水線上“數(shù)據(jù)相關(guān)沖突”的問題，提高了計算效率.設(shè)計中所有操作數(shù)均為符合IEEE 754標(biāo)準(zhǔn)的64bit雙精度浮點數(shù).

1.1 處理單元設(shè)計

PE是構(gòu)成浮點矩陣乘法器的基本單元.每個PE包含一個浮點乘法器、一個浮點加法器和用于存儲計算結(jié)果的存儲單元(Shift register)，其結(jié)構(gòu)如圖1所示.

對于矩陣乘法C=A×B，其中A、B和C分別是M×L、L×N和M×N維矩陣，其計算方法如式(1):

圖1 PE單元結(jié)構(gòu)Fig.1 Structure of PE

圖1中信號分為2類:數(shù)據(jù)信號和控制信號.其中 a、b、MULTI_RESULT、ADD_RESULT、Q_FB 和output都是64 bit數(shù)據(jù)信號，其他都是控制信號，具體功能描述見表1.

表1 FPGA內(nèi)部資源使用情況Table 1 Usage of FPGA internal resources

PE工作時，MULTI_OP_ND置高表示乘法器的輸入(a，b)有效，此時啟動乘法器，當(dāng)浮點乘法器計算完畢時，MULTI_RDY輸出高電平表示MULTI_RESULT有效;當(dāng)ADD_OP_ND變?yōu)楦唠娖綍r啟動浮點加法器，當(dāng)浮點加法器計算完畢時，ADD_RDY輸出高電平表示乘加過程結(jié)束.此外，MULTI_RDY和ADD_RDY還作為Shift register的控制信號，控制Shift register移位.當(dāng)一組元素計算完畢之后，將輸入信號RD_EN置為高電平，讀取Shift register中存儲的結(jié)算結(jié)果;當(dāng)需要計算新的元素時，將NEW_OP置為高電平，Shift register的寄存器全部清零，便可以進(jìn)行新一輪的乘加運算.

1.2 環(huán)路流水線設(shè)計

為了提高計算吞吐率，整個PE采用流水線結(jié)構(gòu).與一般的結(jié)構(gòu)相比，流水線結(jié)構(gòu)能達(dá)到更高的時鐘頻率，但是輸出結(jié)果與輸入之間會有時鐘延遲，延遲的時鐘周期數(shù)等于流水線的級數(shù).流水線的級數(shù)(Latency)越高，乘法器與加法器的工作頻率就越高.

對于無反饋回路的流水線結(jié)構(gòu)來說，輸出結(jié)果相對輸入之間的延遲不會影響整個系統(tǒng)的順序執(zhí)行.但是當(dāng)流水線結(jié)構(gòu)中存在反饋回路時，若不妥善解決延遲問題，流水線上就會出現(xiàn)“數(shù)據(jù)相關(guān)沖突”.以PE的數(shù)據(jù)通道為例，圖1中加法器端口2的輸入數(shù)據(jù)來自于上一次乘積累加操作的結(jié)果，這便構(gòu)成了一個反饋回路，只有保證MULTI_RESULT到達(dá)加法器端口1的時間與上次乘累加的結(jié)果(Q_FB)到達(dá)端口2的時間一致，才能確保PE的正確運行.否則，必然導(dǎo)致流水線時序紊亂，無法完成給定的計算任務(wù)，這就是所謂的“數(shù)據(jù)相關(guān)沖突”.下面通過剖析圖2來闡述這個問題.

圖2 PE的時序波形Fig.2 Timing diagram of PE

如圖2所示，MULTI_RDY比MULTI_OP_ND延遲T1;ADD_RDY比 ADD_OP_ND(即圖2中的MULTI_RDY)延遲T2.設(shè)CLK的周期為T，浮點乘法器和浮點加法器的Latency值分別為u和v，則T1=uT，T2=vT.設(shè)某個計算元素c的前后2組輸入數(shù)據(jù)分別是 a1、b1和 a2、b2.開始計算時 Shift register的寄存器被全部清零，A1=a1b1+0，A1經(jīng)過一個寄存器延遲得到 Q1;同樣，M2=a2b2.由圖2可以看出，只有當(dāng)a2、b2和a1、b1之間保持v+1個時鐘周期的延遲時(T2+T=(v+1)T)，才能保證M2和Q1同時到達(dá)浮點加法器2個輸入端口，進(jìn)而得到正確的累加結(jié)果:A2=Q1+M2.否則，就會導(dǎo)致M2和Q1到達(dá)浮點加法器輸入端口的時間不一致，發(fā)生“數(shù)據(jù)相關(guān)沖突”，無法得到正確的計算結(jié)果.

式中:T0為計算的起始時間，d_inA［t］、d_inB［t］分別表示 t時刻 a、b 端口的輸入數(shù)據(jù)，1≤k≤L，0≤s≤v.由式(2)、(3)可知，上述操作的目的是把求解cij，ci+1，j，…，ci+v，j這 v+1 個計算任務(wù)交叉編排在一條流水線上執(zhí)行.同時，上述操作能保證對于cij的計算任務(wù)而言，前后2組輸入數(shù)據(jù)保持v+1個時鐘周期的延遲，而且不同cij的輸入數(shù)據(jù)不發(fā)生重疊.所以，這樣一條流水線上能完成v+1計算任務(wù)的交叉執(zhí)行，即一個PE單元花費L(T2+T)時間能完成cij，ci+1，j，…，ci+v，j的 v+1 個元素的求解過程.從而，圖 1所示的shift register需要v+1個寄存器來存儲計算結(jié)果，即n=v+1.

1.3 并行結(jié)構(gòu)矩陣乘法器設(shè)計

根據(jù)矩陣乘法的簡單并行算法，在FPGA芯片上實現(xiàn)P個PE單元，這些PE單元按照圖3所示的1×P陣列形式排列.PE單元之間不存在信息交互，它們獨立地完成各自的計算任務(wù).由式(2)和(3)可知，每個PE單元進(jìn)行計算時要用到A的n行和B的某1列數(shù)據(jù)，整個PE陣列一次計算需要用到A的n行和B的P列數(shù)據(jù).將輸入矩陣A、B分別按行和按列進(jìn)行分塊:A=［AT1AT2… ATM］T，B=［B1B2…BN］，其中Ai表示A的第i行，Bj表示B的第j列.將A的n行和B的P列作為圖3所示系統(tǒng)的輸入，圖中預(yù)處理模塊1和預(yù)處理模塊2的功能分別對Ai和Bj進(jìn)行處理，使得它們分別滿足式(2)和(3)中d_inA［t］、d_inB［t］的要求并將其送入各個 PE單元的a、b端口進(jìn)行計算.

圖3 并行矩陣乘法器結(jié)構(gòu)Fig.3 Timing diagram of PE

通過這種PE的陣列結(jié)構(gòu)，可以完成任意維數(shù)的矩陣乘法運算.假設(shè)A和B分別為M×L、L×N維矩陣，對于任意的M、N和L值，可以通過下述算法計算C=A×B的結(jié)果:

從以上算法可以看出，使用并行矩陣乘法器進(jìn)行計算時，循環(huán)的次數(shù)是傳統(tǒng)串行算法的1/P，即計算復(fù)雜度降低為O(n3/P).同時由于該并行矩陣乘法器中的各個PE單元是相互獨立的，因此可以方便地擴(kuò)展到多片F(xiàn)PGA上實現(xiàn)并行計算.

2 矩陣乘法器性能分析

下面以在FPGA上實現(xiàn)的并行矩陣乘法器來對上述設(shè)計的性能進(jìn)行分析.本文選用Xilinx Virtex-5 LX155芯片實現(xiàn)該設(shè)計.PE中的浮點乘法器和浮點加法器使用Xilinx公司提供的floating-point IP核生成.通過對運行速度及該器件中DSP48E單元、CLB單元等資源進(jìn)行綜合考慮，對并行矩陣乘法器進(jìn)行如下設(shè)置:1)IP核生成浮點乘法器時，DSP48E的使用等級設(shè)置為Medium Usage(即單個浮點乘法器使用9個DSP48E單元)，Latency的值設(shè)定為15;2)IP核生成浮點加法器時，DSP48E的使用等級設(shè)置為No Usage，Latency的值設(shè)定為9;3)設(shè)定矩陣乘法器中PE單元的個數(shù)P=10.本設(shè)計中所使用的FPGA開發(fā)環(huán)境和仿真環(huán)境分別為Xilinx ISE Design Suite 13.1 和 Mentor Graphics Modelsim SE 6.5a.

2.1 峰值計算性能

理想情況下，每個PE單元在一個時鐘周期內(nèi)可以完成1次雙精度浮點乘法操作和1次雙精度浮點加法操作，因此整個矩陣乘法器的計算性能可計算為

式中:PERFpeak表示矩陣乘法器的峰值計算性能(每秒百萬次浮點操作)，P為矩陣乘法器中PE單元的個數(shù)，f為矩陣乘法器工作的時鐘頻率.

FPGA內(nèi)部資源的使用情況見表2.根據(jù)布局布線后仿真的結(jié)果，該矩陣乘法器在未做優(yōu)化的情況下工作頻率能達(dá)到250 MHz.由此可知該矩陣乘法器的峰值計算性能可達(dá)到5 000 MFLOPS.

表2 FPGA內(nèi)部資源使用情況Table 2 Usage of FPGA internal resources

2.2 平均計算性能

并行矩陣乘法器的平均計算性能可以通過計算

2個矩陣相乘所需的總時間來求得，如式(5)所示.

式中:PERF表示并行矩陣乘法器的平均計算性能，F(xiàn)表示2個矩陣相乘總共需要完成的雙精度浮點操作次數(shù)，t為計算時間.

本文分別以2個128 bit×128 bit的矩陣相乘和2個256 bit×256 bit的矩陣相乘的實例來分析該設(shè)計的平均計算性能，如表3所示.

表3 平均計算性能對比Table 3 Comparison of computation performance

由表2、3可以看出，本文設(shè)計的并行矩陣乘法器的峰值計算性能可達(dá)到5 000 MFLOPS，平均計算性能可以保證在峰值計算性能的85%左右.而田翔等的設(shè)計A未采用流水線結(jié)構(gòu)，工作頻率只有60 MHz，即使在一片F(xiàn)PGA上集成了25個PE單元，它的峰值計算性能只能達(dá)到3 000 MFLOPS［4］.而且，設(shè)計A的平均浮點計算性能只能保持在峰值計算性能的50%左右.由此可見，本文設(shè)計在計算性能上有大幅度提高.

在矩陣乘法計算中，若FPGA的I/O帶寬小于一定值，并行矩陣乘法器中的PE單元就會出現(xiàn)等待狀態(tài)，此時，帶寬便成為制約計算性能的因素.當(dāng)I/O帶寬達(dá)到或者高于這個值后，每個PE單元的計算性能則成為制約并行矩陣乘法器計算性能的主要因素.在本文的設(shè)計中，PE的計算性能主要由工作頻率決定，在工作頻率為250 MHz的情況下，只要I/O帶寬達(dá)到4 GB/s，便不會對整個系統(tǒng)的計算性能產(chǎn)生影響.

3 結(jié)束語

本文設(shè)計了一個全流水結(jié)構(gòu)的并行雙精度浮點矩陣乘法器，并在Xilinx xc5vlx155 FPGA上實現(xiàn)了該方案.矩陣乘法器內(nèi)部的PE單元采用流水線結(jié)構(gòu)，并運用C-slow時序重排技術(shù)解決了環(huán)路流水線中“數(shù)據(jù)相關(guān)沖突”的問題，提高了計算效率.實驗結(jié)果表明，對于不同維數(shù)的矩陣乘法，本設(shè)計都有較高的計算性能.同時，本文設(shè)計的并行矩陣乘法器結(jié)構(gòu)，其內(nèi)部的各個PE單元相互獨立，因而具有很好的可擴(kuò)展性.在后續(xù)的研究工作中，需要提出更為合理的并行結(jié)構(gòu)，通過多片F(xiàn)PGA的并行計算來進(jìn)一步提高矩陣乘法器的計算性能.

［1］AMIRA A，BENSAALI F.An FPGA based parameterizable system for matrix product implementation［C］//IEEE Workshop on Signal Processing Systems(SPIS’02).San Diego，2002:75-79.

［2］JANG J，CHOI S，PRASANNA V K K.Area and time efficient implementations of matrix multiplication on FPGAs［C］//2002 IEEE International Conference on Field Programmable Technology.Seoul，Korea，2002:93-100.

［3］CAMPBELL S J，KHATRI S P.Resource and delay efficient matrix multiplication using newer FPGA devices［C］//Proceedings of the 16th ACM Great Lakes Symposium on VLSI.Philadelphia，USA，2006:308-311.

［4］田翔，周凡.基于FPGA的實時雙精度浮點矩陣乘法器設(shè)計［J］.浙江大學(xué)學(xué)報:工學(xué)版，2008，42(9):1611-1615.

TIAN Xiang，ZHOU Fan.Design of field programmable gate array based real-time double precision floating-point matrix multiplier［J］.Journal of Zhejiang University:Engineering Science，2008，42(9):1611-1615.

［5］LEISERSON C，ROSE F，SAXE J.Optimizing synchronous circuitry by retiming［C］//Proceedings of the 3rd Caltech Conference On VLSI.Rockville，Maryland，1983:87-116.

［6］SU Ming，ZHOU Lili.Maximizing the throughput-area efficiency of fully-parallel low-density parity-check decoding with C-slow retiming and asynchronous deep pipelining［C］//The 25th International Conference on Computer Design.Washington，DC，USA，2007:93-100.

［7］肖宇，王建業(yè)，張偉.基于IP核的數(shù)選式浮點矩陣相乘設(shè)計［J］.集成電路應(yīng)用，2011，37(6):52-55.

XIAO Yu，WANG Jianye，ZHANG Wei.Floating-point matrix multiplication design based on IP core［J］.Application of Integrated Circuits，2011，37(6):52-55.

［8］許芳，席毅，陳虹.基于FPGA Nios-Ⅱ的矩陣運算硬件加速器設(shè)計［J］.電子測量與儀器學(xué)報，2011，25(4):376-383.

XU Fang，XI Yi，CHEN Hong.Design and implementation of matrix hardware acceleration based on FPGA/Nios-II［J］，Journal of Electronic Measurement and Instrument，2011，25(4):376-383.

［9］黎鐵軍，李秋亮，徐煒遐.一種128位高性能全流水浮點乘加部件［J］.國防科技大學(xué)學(xué)報，2010，32(2):56-60.

LI Tiejun，LI Qiuliang，XU Weixia.A high performance pipeline architecture of 128 bit floating-point fused multiply add unit［J］.Journal of National University of Defense Technology，2010，32(2):56-60.