賈 迅 鄔貴明 謝向輝 吳 東

(數(shù)學(xué)工程與先進(jìn)計算國家重點實驗室 江蘇無錫 214125)

矩陣乘運算核心廣泛應(yīng)用于數(shù)值計算、數(shù)字信號分析、圖像處理、人工智能等領(lǐng)域，其計算形式為

[i,j]+=A[i,k]×B[k,j],

其中，A∈M×K,B∈K×N,C∈M×N.循環(huán)k的計算存在寫后讀數(shù)據(jù)相關(guān)，無法并行；循環(huán)i和循環(huán)j的計算不存在真相關(guān)，可以并行.另外，計算過程中的存儲訪問均比較規(guī)整.由于計算復(fù)雜度為O(n3)，雙精度浮點矩陣乘往往是大規(guī)模數(shù)值計算應(yīng)用的性能瓶頸[1]，其性能直接反映了系統(tǒng)的計算能力.以性能測試程序Linpack[2]為例，程序?qū)崿F(xiàn)了基于雙精度浮點矩陣乘的線性方程求解，并以單位時間內(nèi)系統(tǒng)執(zhí)行浮點操作的次數(shù)(FLOPS)來衡量系統(tǒng)的計算能力.因此，為雙精度浮點矩陣乘運算提供高性能且高效的計算平臺對大規(guī)模數(shù)值計算應(yīng)用的性能發(fā)揮至關(guān)重要.

目前，矩陣乘的計算平臺主要包括CPU，GPGPU，F(xiàn)PGA.其中，CPU和GPGPU平臺上矩陣乘的實現(xiàn)和應(yīng)用較為成熟.處理器廠商會提供針對各自產(chǎn)品架構(gòu)進(jìn)行性能優(yōu)化的計算函數(shù)接口庫，如Intel MKL[3]，IBM ESSL[4]，NVIDIA cuBLAS[5]等.由于處理器架構(gòu)面向通用計算，CPU和GPGPU上矩陣乘的計算效率無法達(dá)到最優(yōu)，如Intel Knights Corner眾核處理器矩陣乘的計算效率僅為89%[6]，NVIDIA P100的計算效率約為93%[7].矩陣乘效率對系統(tǒng)計算效率的影響幾乎是線性的，矩陣乘效率的損失必然導(dǎo)致系統(tǒng)計算效率的降低[8].

FPGA計算平臺具有靈活的可編程性、細(xì)粒度的并行能力和豐富的邏輯資源.同時，矩陣乘運算數(shù)據(jù)并行性好和存儲訪問規(guī)整的算法特點適合并行計算結(jié)構(gòu)的設(shè)計[9]，因此基于FPGA平臺實現(xiàn)針對矩陣乘的并行計算結(jié)構(gòu)成為學(xué)術(shù)研究的熱點，并且已經(jīng)取得一些研究成果[10].以浮點乘加運算部件作為計算單元(processing elements， PE)的核心，多個計算單元組成廣播[11]或陣列結(jié)構(gòu)[12-14]進(jìn)行加速計算是典型的實現(xiàn)方式.其中文獻(xiàn)[12]提出的線性陣列計算結(jié)構(gòu)可以處理任意規(guī)模的矩陣乘運算，且數(shù)據(jù)傳輸和存儲的效率較高.由于矩陣乘的性能受限于FPGA的可用的邏輯資源和時鐘頻率，即使以資源利用和時鐘頻率為優(yōu)化目標(biāo)的結(jié)構(gòu)設(shè)計，其計算性能也只能達(dá)到203.1 GFLOPS[15](252個PE，時鐘頻率為403 MHz).顯然，目前基于FPGA的矩陣乘無法滿足大規(guī)模數(shù)值計算應(yīng)用對計算能力的需求.

近年來，面向特定應(yīng)用的硬件定制加速成為高性能、高效計算結(jié)構(gòu)設(shè)計的一個重要趨勢[16-17].Google公開了TPU(tensor processing unit)芯片的設(shè)計[18]，其核心是整數(shù)矩陣乘脈動陣列計算結(jié)構(gòu).NVIDIA最新發(fā)布的V100 GPGPU[19]芯片中引入了稱為“Tensor Core”的16 b，32 b混合精度矩陣乘陣列計算結(jié)構(gòu)，單芯片集成了672個TensorCore，可提供120 TFLOPS的計算能力.由于TPU和V100芯片中矩陣乘定制加速結(jié)構(gòu)的設(shè)計面向深度學(xué)習(xí)類應(yīng)用，結(jié)構(gòu)支持的計算精度較低[20-21]，因而無法直接應(yīng)用于對計算精度要求較高的數(shù)值計算領(lǐng)域.

本文以線性陣列計算結(jié)構(gòu)為基礎(chǔ)，對其進(jìn)行結(jié)構(gòu)優(yōu)化，并采用硬件定制的方法設(shè)計實現(xiàn)了支持雙精度浮點運算的矩陣乘協(xié)處理器.同時，本文建立了性能模型并深入分析了各結(jié)構(gòu)設(shè)計參數(shù)對協(xié)處理器實際計算性能和效率的影響.此外，本文還驗證了矩陣乘協(xié)處理器的功能正確性并評估了其硬件實現(xiàn)的開銷.本文探索了硬件定制結(jié)構(gòu)設(shè)計在雙精度浮點矩陣乘加速計算中的應(yīng)用，研究成果對提升現(xiàn)有計算系統(tǒng)的性能和效率有一定的借鑒意義.

1 線性陣列計算結(jié)構(gòu)及大規(guī)模矩陣乘算法

線性陣列計算結(jié)構(gòu)[12]如圖1所示，多個計算單元線性互連，每個計算單元包含局部存儲器c以存儲矩陣C的分塊數(shù)據(jù)，寄存器a和b分別存儲矩陣A和矩陣B的數(shù)據(jù)元素，一個雙精度浮點乘加運算部件和一個計數(shù)器CNT.

Fig. 1 Linear array computing architecture圖1 線性陣列計算結(jié)構(gòu)

數(shù)據(jù)b在PE之間進(jìn)行傳輸，其有效時驅(qū)動PE執(zhí)行一次乘加計算：c[CNT]=c[CNT]+a×b.

算法1給出了基于上述線性陣列計算結(jié)構(gòu)的大規(guī)模分塊矩陣乘算法[12].參數(shù)Sp，St分別是矩陣C行數(shù)M和列數(shù)N的分塊粒度，同時也表示線性陣列計算結(jié)構(gòu)中PE的數(shù)目和各PE集成的局部存儲器c的深度.算法中索引為p和t的最內(nèi)層循環(huán)對應(yīng)線性陣列計算結(jié)構(gòu)執(zhí)行的浮點矩乘迭代計算，即Sp個PE實現(xiàn)矩陣A中規(guī)模為Sp×1的矩陣子塊與矩陣B中規(guī)模為1×St的矩陣子塊乘計算，并將計算結(jié)果與矩陣C中規(guī)模為Sp×St的子塊相加.陣列計算結(jié)構(gòu)在單次迭代計算中共執(zhí)行2SpSt次雙精度浮點運算.

Fig. 2 Data block mapping for the linear array圖2 矩陣分塊數(shù)據(jù)到陣列計算結(jié)構(gòu)的映射

算法1. 基于線性陣列結(jié)構(gòu)的大規(guī)模矩陣乘算法.

輸入：雙精度浮點矩陣A、矩陣B、矩陣C；

輸出：雙精度浮點矩陣C.

① ForTp=0 toM-1,Sp

② ForTt=0 toN-1,St

③ LoadC[Tp:Tp+Sp-1,Tt:Tt+St-1]；

④ Fork=0 toK-1

⑤ LoadA[Tp:Tp+Sp-1,k]；

⑥ LoadB[k,Tt:Tt+St-1]；

⑦ Forp=TptoTp+Sp-1

⑧ Fort=TttoTt+St-1

⑨C[p,t]+=A[p,k]×B[k,t]；

⑩ End For

算法1中索引為Tp，Tt，k的3層循環(huán)在線性陣列計算結(jié)構(gòu)進(jìn)行矩陣乘迭代計算的同時，處理矩陣分塊數(shù)據(jù)的加載與寫回，其功能在硬件上由數(shù)據(jù)傳輸控制部件DTC(data transfer controller)實現(xiàn).各矩陣分塊數(shù)據(jù)到線性陣列計算結(jié)構(gòu)的映射如圖2所示：

線性陣列結(jié)構(gòu)經(jīng)過K輪迭代計算才能得到矩陣C中一個分塊數(shù)據(jù)的計算結(jié)果.因此，在K輪迭代計算開始前，數(shù)據(jù)傳輸控制部件需要先將矩陣C中規(guī)模為Sp×St的分塊數(shù)據(jù)C[Tp:Tp+Sp-1,Tt:Tt+St-1]分別加載到Sp個計算單元的深度為St的局部存儲器c中，并在K輪迭代計算完成后將對應(yīng)的分塊計算結(jié)果寫回主存.另外，每次迭代計算中各PE需要訪問矩陣A中相同的數(shù)據(jù)元素.因此，每次迭代計算開始前，數(shù)據(jù)傳輸控制部件還需要將矩陣A中規(guī)模為Sp×1的分塊數(shù)據(jù)A[Tp:Tp+Sp-1,k]加載到Sp個計算單元的寄存器a中.

矩陣A和矩陣C的分塊數(shù)據(jù)全部加載完成后，數(shù)據(jù)傳輸控制部件再加載矩陣B的數(shù)據(jù)元素以驅(qū)動陣列計算結(jié)構(gòu)進(jìn)行迭代計算.矩陣乘運算執(zhí)行完成共需K×(MSp)×(NSt)次迭代計算.對計算過程中數(shù)據(jù)傳輸控制部件傳輸矩陣A，B，C分塊數(shù)據(jù)的數(shù)據(jù)量、時間間隔以及次數(shù)的量化分析如表1所示：

Table 1 Analysis on Data Transfer for Matrix Multiplication表1 針對矩陣乘運算過程中數(shù)據(jù)傳輸?shù)姆治?/p>

文獻(xiàn)[12]基于Xilinx公司Virtex-5系列的FPGA實現(xiàn)了上述線性陣列計算結(jié)構(gòu).PE個數(shù)Sp設(shè)置為96、局部存儲器深度St設(shè)置為1 024的情況下，F(xiàn)PGA的最大運行頻率為158.08 MHz.輸入矩陣的規(guī)模M，N，K均為8 192時，陣列計算結(jié)構(gòu)完成雙精度浮點矩陣乘運算所需的時間為40.48 s，即結(jié)構(gòu)的計算性能為27.16 GFLOPS，效率為89.49%.

線性陣列計算結(jié)構(gòu)雖然存儲效率較高，但每次迭代計算中各PE需要等待計算所需的矩陣分塊數(shù)據(jù)全部加載完成后才能開始計算，計算與存儲訪問在時間上是串行的，計算資源的閑置必然導(dǎo)致計算效率的損失，高效計算結(jié)構(gòu)的設(shè)計必須對此進(jìn)行改進(jìn).由于結(jié)構(gòu)設(shè)計面向FPGA計算平臺，文獻(xiàn)[12]僅以FPGA可用的邏輯資源作為設(shè)計約束，并未考慮局部存儲器的深度、輸入矩陣的規(guī)模等參數(shù)對線性陣列計算結(jié)構(gòu)實際計算性能和效率的影響，而這些參數(shù)在定制協(xié)處理器的設(shè)計中尤為重要.另外，本文協(xié)處理器的設(shè)計除了考慮核心的計算結(jié)構(gòu)外，還需要考慮其與主機(jī)端的接口，以實現(xiàn)控制信息和計算數(shù)據(jù)的傳輸，從而構(gòu)成完整的加速計算系統(tǒng).

2 協(xié)處理器的設(shè)計

2.1 線性陣列計算結(jié)構(gòu)的優(yōu)化

根據(jù)算法1給出的大規(guī)模分塊矩陣乘運算算法，線性陣列計算結(jié)構(gòu)在進(jìn)行矩陣乘迭代計算前，需要先從主存讀入矩陣A和C的分塊數(shù)據(jù)；K輪迭代計算完成后，需要先將矩陣C的分塊計算結(jié)果寫回主存，才能開始下一輪的迭代計算.矩陣分塊數(shù)據(jù)在主存和線性陣列計算結(jié)構(gòu)之間進(jìn)行傳輸?shù)倪^程中，結(jié)構(gòu)中各計算單元無法進(jìn)行計算，計算資源被閑置，從而導(dǎo)致線性陣列計算結(jié)構(gòu)整體的效率較低.針對這個問題，本文采用雙緩沖技術(shù)對陣列計算結(jié)構(gòu)中各計算單元的設(shè)計進(jìn)行了優(yōu)化.優(yōu)化后，計算單元的結(jié)構(gòu)如圖3所示：

Fig. 3 Structure of PE optimized with double buffering圖3 采用雙緩沖優(yōu)化技術(shù)的PE結(jié)構(gòu)

圖3所示的PE結(jié)構(gòu)中除了已有的寄存器a，b和局部存儲器c之外，還包括雙緩沖優(yōu)化引入的存儲矩陣A，B數(shù)據(jù)元素的寄存器a_buf，b_buf以及存儲矩陣C分塊數(shù)據(jù)的局部存儲器c_buf.

應(yīng)用雙緩沖優(yōu)化技術(shù)后，陣列計算結(jié)構(gòu)可以在使用當(dāng)前緩沖中的數(shù)據(jù)進(jìn)行迭代計算的同時，提前將下一輪迭代計算所需的矩陣分塊數(shù)據(jù)裝入到另一緩沖，實現(xiàn)計算與訪存的重疊，隱藏存儲訪問的開銷.此時，高效的存儲訪問調(diào)度成為應(yīng)用雙緩沖優(yōu)化提升計算結(jié)構(gòu)資源利用率的關(guān)鍵.

2.2 存儲訪問調(diào)度

由表1列出的矩陣乘運算過程中矩陣分塊數(shù)據(jù)的傳輸可知，線性陣列計算結(jié)構(gòu)在進(jìn)行矩陣乘運算時，每個時鐘周期僅需訪問矩陣B的一個數(shù)據(jù)元素，其對存儲帶寬的需求較低；而矩陣A，C的訪問時間間隔分別為St和K×St個時鐘周期，因此可以利用讀取矩陣B數(shù)據(jù)元素剩余的存儲訪問帶寬提前加載矩陣A，C的分塊數(shù)據(jù)至PE的緩沖中，或?qū)⒕仃嘋的分塊計算結(jié)果從緩沖寫回主存.

算法2給出了基于雙緩沖的存儲訪問調(diào)度算法.數(shù)據(jù)傳輸控制部件在硬件上實現(xiàn)了該算法，并支持雙緩沖讀寫訪問的控制.

算法2. 基于雙緩沖的存儲訪問調(diào)度算法.

輸入：雙精度浮點矩陣A、矩陣B、矩陣C；

輸出：雙精度浮點矩陣C.

① LoadC[0:Sp-1,0:St-1]；

② LoadA[0:Sp-1,0]；

③ ForTp=0 toM-1,Sp

④ ForTt=0 toN-1,St

⑤ Fork=0 toK-1

⑥ LoadB[k,Tt:Tt+St-1]；

LoadA[Tp:Tp+Sp-1,k+1] or

StoreC[Tp:Tp+Sp-1,Tt-St:Tt-1] or

LoadC[Tp:Tp+Sp-1,Tt+St-1:Tt+2St-1]；

⑦ Forp=TptoTp+Sp-1

⑧ Fort=TttoTt+St-1

⑨C[p,t]+=A[p,k]×B[k,t]；

⑩ End For

算法2中，陣列計算結(jié)構(gòu)進(jìn)行第一次迭代計算所需的矩陣A，C的分塊數(shù)據(jù)先被讀入.讀取矩陣B的數(shù)據(jù)元素進(jìn)行迭代計算時，優(yōu)先利用剩余的存儲訪問帶寬提前讀取下一次迭代計算所需的矩陣A的分塊數(shù)據(jù)；待矩陣A的分塊數(shù)據(jù)讀入完成，再利用剩余帶寬寫回前K輪迭代計算得到的矩陣C的分塊計算結(jié)果；待計算結(jié)果寫回完成，最后利用剩余帶寬裝入后K輪迭代計算所需的矩陣C的分塊數(shù)據(jù).所有迭代計算完成后，矩陣C中最后一個分塊計算結(jié)果被寫回.

理想情況下，線性陣列計算結(jié)構(gòu)在開始新一次的迭代計算時，各PE所需的矩陣數(shù)據(jù)均已在緩沖中，即存儲訪問的延遲可以有效地被計算隱藏，結(jié)構(gòu)中的計算資源可以得到充分利用.而實際的存儲訪問調(diào)度受限于存儲帶寬、計算單元數(shù)、局部存儲器的深度以及矩陣規(guī)模等參數(shù)，其對矩陣乘協(xié)處理器實際性能的影響將在2.3節(jié)進(jìn)行分析.

2.3 系統(tǒng)結(jié)構(gòu)

基于矩陣乘協(xié)處理器的加速計算系統(tǒng)的整體結(jié)構(gòu)如圖4所示.系統(tǒng)由主機(jī)端處理器芯片和雙精度浮點矩陣乘協(xié)處理器芯片2部分組成.其中，矩陣乘協(xié)處理器通過IO總線與主機(jī)端處理器相連.系統(tǒng)采用PCIe-3.0×16作為芯片間的互連接口，雙向通信帶寬最高為31.5 GBps.

Fig. 4 System architecture based on matrix multiplication coprocessor圖4 基于矩陣乘協(xié)處理器計算系統(tǒng)的架構(gòu)

主機(jī)端和協(xié)處理器端存儲空間相互獨立，待計算的矩陣數(shù)據(jù)存儲在協(xié)處理器端的內(nèi)存.本文協(xié)處理器的設(shè)計面向大規(guī)模矩陣乘，支持的矩陣規(guī)模至少為16 384×16 384.因此，主機(jī)端與協(xié)處理器端的內(nèi)存容量配置為8 GB，并采用DDR3-2133[22]規(guī)格的存儲芯片，接口的數(shù)據(jù)寬度為128 b，帶寬為34.13 GBps.

系統(tǒng)對主機(jī)端處理器的類型沒有限制，只要支持PCIe傳輸協(xié)議即可，其指令集系統(tǒng)架構(gòu)可以是x86，Power，ARM，RISC-V[23]等.矩陣乘協(xié)處理器由應(yīng)用雙緩沖優(yōu)化的線性陣列計算結(jié)構(gòu)、存儲訪問控制器、協(xié)處理器與主機(jī)端接口(coprocessor-host interface， CHI)3部分構(gòu)成.乘加運算部件是計算單元的核心，本文采用全流水且兼容IEEE-754浮點數(shù)標(biāo)準(zhǔn)[24]的乘加運算部件設(shè)計，以保證協(xié)處理器的計算性能.CHI接口實現(xiàn)主機(jī)端與協(xié)處理器之間控制和數(shù)據(jù)的傳輸.主機(jī)端向矩陣乘協(xié)處理器發(fā)送的主要命令如下：

1)Initialization對協(xié)處理器的控制與狀態(tài)寄存器進(jìn)行初始化；

2)DataIn將待計算矩陣數(shù)據(jù)從主機(jī)端內(nèi)存加載到協(xié)處理器端內(nèi)存；

3)MatrixMultiply啟動協(xié)處理器上的雙精度浮點矩陣乘運算；

4)DataOut將矩陣乘運算的結(jié)果從協(xié)處理器端內(nèi)存讀出到主機(jī)端內(nèi)存.

接收到來自主機(jī)端發(fā)送的命令，CHI接口對命令進(jìn)行分析并生成線性陣列計算結(jié)構(gòu)相關(guān)的控制信號.協(xié)處理器在進(jìn)行矩陣乘運算的過程中，主機(jī)端與協(xié)處理器之間不需要進(jìn)行數(shù)據(jù)傳輸，從而簡化了協(xié)處理器的設(shè)計與調(diào)試.

應(yīng)用程序在計算系統(tǒng)上加速矩陣乘運算的過程如圖5所示.當(dāng)矩陣乘協(xié)處理器連接至主機(jī)端時，主機(jī)端發(fā)送Initialization命令對協(xié)處理器進(jìn)行初始化，如設(shè)置浮點計算的舍入模式、溢出控制等.應(yīng)用程序運行在主機(jī)端，當(dāng)其調(diào)用雙精度浮點矩陣乘運算核心，如基礎(chǔ)線性代數(shù)子程序庫(basic linear algebra subroutines， BLAS)中的函數(shù)blas_dgemm(M,N,K,A,B,C,…)時，主機(jī)端通過CHI接口向協(xié)處理器發(fā)送DataIn命令和函數(shù)調(diào)用傳入的矩陣規(guī)模、矩陣數(shù)據(jù)在內(nèi)存中的起始地址等參數(shù).協(xié)處理器的數(shù)據(jù)傳輸控制部件將待計算的矩陣數(shù)據(jù)從主機(jī)端內(nèi)存以直接內(nèi)存訪問(direct memory access， DMA)方式傳輸至協(xié)處理器端內(nèi)存.數(shù)據(jù)傳輸完成后，主機(jī)端發(fā)送MatrixMultiply命令以啟動協(xié)處理器上的矩陣乘運算.運算過程中，數(shù)據(jù)傳輸控制部件按算法2設(shè)計的存儲訪問調(diào)度實現(xiàn)矩陣分塊數(shù)據(jù)在協(xié)處理器端內(nèi)存和線性陣列計算結(jié)構(gòu)之間的傳輸.

Fig. 5 Computation based on matrix multiplication coprocessor圖5 基于矩陣乘協(xié)處理器的計算過程

協(xié)處理器完成矩陣乘運算后，CHI接口向主機(jī)端發(fā)出中斷信號.主機(jī)端收到中斷后，向CHI接口發(fā)送DataOut命令，協(xié)處理器的數(shù)據(jù)傳輸控制部件以DMA的方式將矩陣乘運算的結(jié)果從協(xié)處理器端內(nèi)存?zhèn)鬏敾刂鳈C(jī)端內(nèi)存.運算結(jié)果傳輸完成后，矩陣乘運算調(diào)用結(jié)束，應(yīng)用程序繼續(xù)執(zhí)行.

3 性能模型

為了量化分析雙精度浮點矩陣乘協(xié)處理器的各結(jié)構(gòu)設(shè)計參數(shù)對其實際計算性能和效率的影響，本文基于算法2設(shè)計的存儲訪問調(diào)度算法構(gòu)建了協(xié)處理器的性能模型.模型中的基本參數(shù)包括線性陣列計算結(jié)構(gòu)中PE的個數(shù)Sp、PE集成的局部存儲器的深度St以及輸入矩陣的規(guī)模M，N，K.此外，本文將模型中協(xié)處理器的時鐘頻率記為Freq(單位為GHz)，協(xié)處理器端內(nèi)存的訪問帶寬記為Band(單位為GBps)，矩陣乘運算核心在協(xié)處理器上執(zhí)行完成所需的時間記為T(單位為s)，協(xié)處理器的實際計算性能和效率記為Perf(單位為GFLOPS)和Eff.

當(dāng)計算過程中的存儲訪問延遲完全被計算隱藏，協(xié)處理器的計算性能可以充分發(fā)揮，此時的計算效率最高.就本文提出的存儲訪問調(diào)度算法，每次迭代計算提前讀入矩陣A中規(guī)模為Sp的分塊數(shù)據(jù)的延遲與讀入矩陣B中規(guī)模為St的分塊數(shù)據(jù)的延遲之和應(yīng)當(dāng)小于等于對應(yīng)的計算延遲，即St個時鐘周期.同理，K輪迭代計算的總存儲訪問延遲應(yīng)當(dāng)小于等于對應(yīng)的總計算延遲，即K×St個時鐘周期.存儲訪問的數(shù)據(jù)粒度為8 B的雙精度浮點數(shù)，上述2個約束條件可表示為

(1)

協(xié)處理器端DDR3-2133存儲芯片的數(shù)據(jù)傳輸率為2.133 GTs，單次傳輸?shù)臄?shù)據(jù)寬度為128 b.若協(xié)處理器的工作時鐘頻率Freq=1.5 GHz，則其每個時鐘周期可訪問2個雙精度浮點數(shù).此時，式(1)的約束條件可進(jìn)一步改寫為

(2)

若參數(shù)Sp，St，K可以滿足式(2)中的約束條件，則計算過程中的存儲訪問延遲可以完全被計算隱藏.此時，矩陣乘運算核心在協(xié)處理器上的執(zhí)行完成所需的時間T最小.其由3部分時間構(gòu)成，包括讀入矩陣A和C第1個子塊數(shù)據(jù)的時間、迭代計算的時間以及將矩陣C中最后1個子塊的計算結(jié)果寫回內(nèi)存的時間.計算時間T可表示為

(3)

若參數(shù)Sp，St，K不滿足式(2)的約束條件，存在St1-SpSt這2種可能.下文將分別分析這2種情況下矩陣乘運算核心在協(xié)處理器上的執(zhí)行時間.

對于第1種情況，提前讀入矩陣A子塊數(shù)據(jù)的存儲訪問延遲無法完全被計算時間隱藏，每次迭代計算的時間增加(Sp-St)2個時鐘周期.根據(jù)算法2，當(dāng)前迭代計算執(zhí)行過程中，寫回前K輪迭代計算中矩陣C子塊計算結(jié)果和提前讀入后K輪迭代計算中矩陣C子塊計算數(shù)據(jù)的存儲訪問延遲均無法被計算隱藏，每K輪迭代計算的時間增加Sp×St個時鐘周期.此時，矩陣乘運算核心在協(xié)處理器上執(zhí)行完成所需的時間T可表示為

(4)

對于第2種情況，當(dāng)前迭代計算的執(zhí)行時間可以完全隱藏提前讀入矩陣A子塊數(shù)據(jù)的存儲訪問延遲，但只能部分隱藏寫回前K輪迭代計算中矩陣C子塊計算結(jié)果和提前讀入后K輪迭代計算中矩陣C子塊計算數(shù)據(jù)的存儲訪問延遲，其導(dǎo)致每K輪迭代計算的時間增加Sp×St-K×(St-Sp)2個時鐘周期.此時，矩陣乘運算核心在協(xié)處理器上執(zhí)行完成所需的時間T可表示為

(5)

協(xié)處理器完成矩陣乘核心的加速計算共需執(zhí)行2MNK次雙精度浮點操作.基于對3種情況下矩陣乘運算核心在協(xié)處理器上執(zhí)行時間T的分析結(jié)果，可以構(gòu)建式(6)所示的協(xié)處理器的性能模型.此時，協(xié)處理器的計算效率可表示為實際計算性能與峰值計算性能的比值，即Eff=Perf(2×Sp×Freq).

4 實驗與結(jié)果

本節(jié)以第3節(jié)提出的性能模型為基礎(chǔ)，分析矩陣乘協(xié)處理器結(jié)構(gòu)設(shè)計中雙緩沖優(yōu)化技術(shù)、局部存儲器的深度和輸入矩陣規(guī)模等參數(shù)對協(xié)處理器實際計算性能和效率的影響.實驗假設(shè)線性陣列計算結(jié)構(gòu)的PE數(shù)為1 024，協(xié)處理器的時鐘頻率為1.5 GHz，待計算的矩陣均為方陣，即M=N=K=n.此外，本節(jié)還驗證了矩陣乘協(xié)處理器功能實現(xiàn)的正確性；在主流工藝下評估了其硬件實現(xiàn)的開銷，并與同等工藝的GPGPU芯片進(jìn)行了對比，以探索其硬件可實現(xiàn)性和在實際應(yīng)用中的優(yōu)勢.

4.1 雙緩沖優(yōu)化技術(shù)對矩陣乘協(xié)處理器性能的影響

實驗在不同的矩陣規(guī)模下，分別統(tǒng)計了應(yīng)用雙緩沖優(yōu)化技術(shù)前后矩陣乘協(xié)處理器的性能，以分析應(yīng)用雙緩沖優(yōu)化技術(shù)對協(xié)處理器性能的影響.實驗中，各PE局部存儲器的深度均設(shè)置為2 048，實驗的結(jié)果如圖6所示.

從圖6可以看出，不同矩陣規(guī)模下，雙緩沖優(yōu)化技術(shù)的應(yīng)用均可以提升矩陣乘協(xié)處理器的計算性能，協(xié)處理器的性能平均提升30.06%.矩陣規(guī)模為4 096時，性能提升最為明顯，達(dá)到了45.45%.可見，本文設(shè)計的存儲訪問調(diào)度算法可以在計算過程中利用雙緩沖優(yōu)化有效隱藏存儲訪問的延遲.

同時，雙緩沖優(yōu)化對矩陣乘協(xié)處理器的性能提升受矩陣規(guī)模的影響.對此分析如下：當(dāng)矩陣規(guī)模較小時，計算無法有效隱藏存儲訪問的延遲，應(yīng)用雙緩沖優(yōu)化對性能的提升作用有限；隨著矩陣規(guī)模逐漸增加，計算可以隱藏更多的訪存延遲，此時應(yīng)用雙緩沖優(yōu)化可以顯著提升協(xié)處理器的計算性能；當(dāng)矩陣規(guī)模較大時，計算可以很好隱藏存儲訪問的延遲，而實驗中局部存儲器的深度保持不變，雙緩沖優(yōu)化對性能的提升作用逐漸降低.

Fig. 6 Effect of double-buffering optimization technique on coprocessor performance圖6 雙緩沖優(yōu)化技術(shù)對矩陣乘協(xié)處理器性能的提升

4.2 局部存儲器深度對矩陣乘協(xié)處理器性能的影響

Fig. 7 Effect of local memory depth on coprocessor performance圖7 局部存儲器深度對矩陣乘協(xié)處理器性能的影響

局部存儲器的深度直接影響協(xié)處理器硬件實現(xiàn)的開銷，因而是矩陣乘協(xié)處理器結(jié)構(gòu)設(shè)計中一個重要的參數(shù).實驗在128，1 024，8 192，16 384這4種矩陣規(guī)模下，分別將計算單元局部存儲器的深度h從128增加至8 192，統(tǒng)計矩陣乘協(xié)處理器實際計算性能的變化情況，實驗結(jié)果如圖7所示:

從圖7可以看出，4種矩陣規(guī)模下，矩陣乘協(xié)處理器的實際計算性能均隨局部存儲器深度的增加而提升，局部存儲器深度增加大于2 048時，矩陣乘協(xié)處理器的性能基本保持不變.協(xié)處理器實際計算性能隨局部存儲器深度的變化可以根據(jù)式(2)中計算完全隱藏存儲訪問延遲的約束條件進(jìn)行分析.當(dāng)St

4.3 矩陣規(guī)模對矩陣乘協(xié)處理器計算效率的影響

盡管本文矩陣乘協(xié)處理器的設(shè)計面向數(shù)值計算領(lǐng)域的大規(guī)模應(yīng)用，矩陣規(guī)模對協(xié)處理器實際計算性能和效率的影響仍然值得關(guān)注.實驗在應(yīng)用雙緩沖優(yōu)化、局部存儲器深度設(shè)置為2 048時，統(tǒng)計了不同矩陣規(guī)模下，矩陣乘協(xié)處理器實際的計算性能和效率，以及存儲訪問延遲占矩陣乘運算核心總執(zhí)行時間的比例，實驗結(jié)果如表2所示:

Table 2 Effect of Matrix Size on Coprocessor Performance

從表2中的數(shù)據(jù)可以看出，隨著矩陣規(guī)模的增加，矩陣乘協(xié)處理器的實際計算效率不斷提高：矩陣規(guī)模為128時，協(xié)處理器的計算效率僅為6.25%；矩陣規(guī)模增加至4 096時，協(xié)處理器的計算效率達(dá)到96.97%；矩陣規(guī)模大于等于8 192時，協(xié)處理器的計算效率可以達(dá)到99%.

對協(xié)處理器計算效率與矩陣規(guī)模相關(guān)性的分析如下：當(dāng)PE數(shù)Sp=1 024、局部存儲器深度St=2 048時，式(3)表示的計算完全隱藏訪存延遲的約束條件是否滿足僅取決于計算中的矩陣規(guī)模K.隨著矩陣規(guī)模的增大，計算可以更好地隱藏存儲訪問的延遲，協(xié)處理器的計算效率也進(jìn)一步提升.表2中存儲訪問延遲占矩陣乘總執(zhí)行時間的比例驗證了上述分析，即較大的矩陣規(guī)模對應(yīng)較低的存儲訪問延遲占比.

4.4 功能正確性驗證與硬件實現(xiàn)開銷的評估

本文基于Verilog硬件描述語言實現(xiàn)了雙精度浮點矩陣乘協(xié)處理器.其中，線性陣列計算結(jié)構(gòu)中各計算單元采用雙緩沖優(yōu)化技術(shù)，PE數(shù)設(shè)置為1 024、局部存儲器的深度設(shè)置為2 048.

對于矩陣乘協(xié)處理器功能正確性的驗證，本文以Linpack性能測試程序HPL[25]隨機(jī)生成的雙精度浮點矩陣數(shù)據(jù)作為測試數(shù)據(jù)集，通過調(diào)用OpenBLAS[26-27]線性代數(shù)計算函數(shù)庫實現(xiàn)了雙精度浮點矩陣乘運算核心.本文將OpenBLAS的計算結(jié)果與矩陣乘協(xié)處理器的仿真計算結(jié)果進(jìn)行了比對，從而驗證了協(xié)處理器功能實現(xiàn)的正確性.

對于矩陣乘協(xié)處理器硬件開銷的評估，本文基于Synopsys Design Compiler設(shè)計工具在1.5 GHz、主流工藝下對協(xié)處理器的硬件設(shè)計進(jìn)行了邏輯綜合，綜合過程中的參數(shù)配置如表3所示：

Table 3 Configurations for Synopsys Design Compiler表3 邏輯綜合的參數(shù)設(shè)置

根據(jù)Design Compiler輸出的綜合結(jié)果，矩陣乘協(xié)處理器的硬件設(shè)計完全滿足1.5 GHz頻率下的時序要求.其功耗為38.57 W，總面積為116.26 mm2.其中，組合邏輯的面積為16.78 mm2，非組合邏輯的面積為99.48 mm2.

本文將矩陣乘協(xié)處理器與采用同等工藝的NVIDIA K40 GPGPU[19]芯片進(jìn)行了對比，比較的內(nèi)容包括時鐘頻率、峰值計算性能，進(jìn)行矩陣乘運算的實際性能和效率以及硬件開銷等，結(jié)果如表4所示.針對雙精度浮點矩陣乘運算，本文協(xié)處理器的實際性能是K40 GPGPU芯片的1.95倍，而面積開銷僅為后者的21.05%(不考慮GPGPU芯片中PCIe、圖形顯示和存儲控制器等接口的硬件開銷).

Table 4 Coprocessor versus K40 GPGPU on Performance

5 結(jié)束語

矩陣乘運算廣泛應(yīng)用于科學(xué)與工程領(lǐng)域.因計算復(fù)雜度較高，且在CPU，GPGPU，F(xiàn)PGA等現(xiàn)有計算平臺上的性能和效率受限，雙精度浮點矩陣乘運算往往成為大規(guī)模數(shù)值計算應(yīng)用的性能瓶頸.近年來，面向應(yīng)用的硬件定制成為高性能、高效計算結(jié)構(gòu)設(shè)計的重要趨勢.本文將這一結(jié)構(gòu)設(shè)計方法應(yīng)用到矩陣乘的加速計算中，設(shè)計了雙精度浮點矩陣乘定制協(xié)處理器和基于協(xié)處理器的加速計算系統(tǒng).

本文雙精度浮點矩陣乘協(xié)處理器的設(shè)計以線性陣列計算結(jié)構(gòu)為基礎(chǔ)，應(yīng)用了雙緩沖技術(shù)并對存儲訪問調(diào)度進(jìn)行了優(yōu)化.同時，本文構(gòu)建了針對協(xié)處理器的性能模型，并基于性能模型深入分析了協(xié)處理器的各結(jié)構(gòu)設(shè)計參數(shù)對其計算性能和效率的影響.此外，本文還驗證了矩陣乘協(xié)處理器功能實現(xiàn)的正確性，并在主流工藝下評估了其硬件開銷.本文雙精度浮點矩陣乘協(xié)處理器的實際計算性能可達(dá)3 TFLOPS，計算效率最高可達(dá)99%，有效實現(xiàn)了大規(guī)模數(shù)值計算應(yīng)用的加速.與同等工藝的GPGPU計算平臺相比，本文設(shè)計的矩陣乘協(xié)處理器在計算性能和硬件開銷方面均具有明顯優(yōu)勢.

本文探索了硬件定制的結(jié)構(gòu)設(shè)計在高性能計算中的應(yīng)用.后續(xù)研究將基于BLAS庫和Linpack性能測試程序評估本文雙精度浮點矩陣乘協(xié)處理器的應(yīng)用對大規(guī)模計算系統(tǒng)性能和效率的提升效果；同時分析協(xié)處理器實際性能隨陣列計算結(jié)構(gòu)中計算單元數(shù)增加的可擴(kuò)展性，從而進(jìn)一步優(yōu)化協(xié)處理器的結(jié)構(gòu)設(shè)計.另外，協(xié)處理器的應(yīng)用適應(yīng)性也是未來研究工作的重點.