賈 迅 鄔貴明 謝向輝 吳 東

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

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

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

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

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

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

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

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

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

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

Fig. 1 Linear array computing architecture圖1 線性陣列計算結構

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

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

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

算法1. 基于線性陣列結構的大規(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)在線性陣列計算結構進行矩陣乘迭代計算的同時，處理矩陣分塊數(shù)據(jù)的加載與寫回，其功能在硬件上由數(shù)據(jù)傳輸控制部件DTC(data transfer controller)實現(xiàn).各矩陣分塊數(shù)據(jù)到線性陣列計算結構的映射如圖2所示：

線性陣列結構經(jīng)過K輪迭代計算才能得到矩陣C中一個分塊數(shù)據(jù)的計算結果.因此，在K輪迭代計算開始前，數(shù)據(jù)傳輸控制部件需要先將矩陣C中規(guī)模為Sp×St的分塊數(shù)據(jù)C[Tp:Tp+Sp-1,Tt:Tt+St-1]分別加載到Sp個計算單元的深度為St的局部存儲器c中，并在K輪迭代計算完成后將對應的分塊計算結果寫回主存.另外，每次迭代計算中各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ù)元素以驅動陣列計算結構進行迭代計算.矩陣乘運算執(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>

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

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

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

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

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

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

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

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

2.2 存儲訪問調度

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

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

算法2. 基于雙緩沖的存儲訪問調度算法.

輸入：雙精度浮點矩陣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中，陣列計算結構進行第一次迭代計算所需的矩陣A，C的分塊數(shù)據(jù)先被讀入.讀取矩陣B的數(shù)據(jù)元素進行迭代計算時，優(yōu)先利用剩余的存儲訪問帶寬提前讀取下一次迭代計算所需的矩陣A的分塊數(shù)據(jù)；待矩陣A的分塊數(shù)據(jù)讀入完成，再利用剩余帶寬寫回前K輪迭代計算得到的矩陣C的分塊計算結果；待計算結果寫回完成，最后利用剩余帶寬裝入后K輪迭代計算所需的矩陣C的分塊數(shù)據(jù).所有迭代計算完成后，矩陣C中最后一個分塊計算結果被寫回.

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

2.3 系統(tǒng)結構

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

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

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

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

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

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

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

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

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

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

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

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

3 性能模型

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

當計算過程中的存儲訪問延遲完全被計算隱藏，協(xié)處理器的計算性能可以充分發(fā)揮，此時的計算效率最高.就本文提出的存儲訪問調度算法，每次迭代計算提前讀入矩陣A中規(guī)模為Sp的分塊數(shù)據(jù)的延遲與讀入矩陣B中規(guī)模為St的分塊數(shù)據(jù)的延遲之和應當小于等于對應的計算延遲，即St個時鐘周期.同理，K輪迭代計算的總存儲訪問延遲應當小于等于對應的總計算延遲，即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)的約束條件可進一步改寫為

(2)

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

(3)

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

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

(4)

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

(5)

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

4 實驗與結果

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

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

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

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

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

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

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

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

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

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

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

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

Table 2 Effect of Matrix Size on Coprocessor Performance

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

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

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

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

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

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

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

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

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

Table 4 Coprocessor versus K40 GPGPU on Performance

5 結束語

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

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

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