劉沛華，魯華祥，龔國(guó)良，劉文鵬

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

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

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

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

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

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

1.1 處理單元設(shè)計(jì)

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

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

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

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

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

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

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

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

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

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

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

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

根據(jù)矩陣乘法的簡(jiǎn)單并行算法，在FPGA芯片上實(shí)現(xiàn)P個(gè)PE單元，這些PE單元按照?qǐng)D3所示的1×P陣列形式排列.PE單元之間不存在信息交互，它們獨(dú)立地完成各自的計(jì)算任務(wù).由式(2)和(3)可知，每個(gè)PE單元進(jìn)行計(jì)算時(shí)要用到A的n行和B的某1列數(shù)據(jù)，整個(gè)PE陣列一次計(jì)算需要用到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的功能分別對(duì)Ai和Bj進(jìn)行處理，使得它們分別滿足式(2)和(3)中d_inA［t］、d_inB［t］的要求并將其送入各個(gè) PE單元的a、b端口進(jìn)行計(jì)算.

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

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

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

2 矩陣乘法器性能分析

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

2.1 峰值計(jì)算性能

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

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

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

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

2.2 平均計(jì)算性能

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

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

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

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

表3 平均計(jì)算性能對(duì)比Table 3 Comparison of computation performance

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

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

3 結(jié)束語(yǔ)

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

［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í)時(shí)雙精度浮點(diǎn)矩陣乘法器設(shè)計(jì)［J］.浙江大學(xué)學(xué)報(bào):工學(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ù)選式浮點(diǎn)矩陣相乘設(shè)計(jì)［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-Ⅱ的矩陣運(yùn)算硬件加速器設(shè)計(jì)［J］.電子測(cè)量與儀器學(xué)報(bào)，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位高性能全流水浮點(diǎn)乘加部件［J］.國(guó)防科技大學(xué)學(xué)報(bào)，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.