王晞陽(yáng)，陳繼林，李猛，劉首文

（1.國(guó)家超級(jí)計(jì)算無(wú)錫中心，江蘇 無(wú)錫 214072；2.中國(guó)電力科學(xué)研究院有限公司，北京 100192；3.國(guó)網(wǎng)湖北省電力有限公司，武漢 430070）

0 概述

在電力系統(tǒng)仿真中，電磁暫態(tài)仿真是主要的系統(tǒng)應(yīng)用之一，也是電力系統(tǒng)安全分析和運(yùn)行的關(guān)鍵組件［1］，該應(yīng)用的核心算法是對(duì)大規(guī)模線性方程組Ax=b進(jìn)行求解。通過(guò)對(duì)實(shí)際數(shù)據(jù)的分析可知，線性方程組中的系數(shù)矩陣A通常為稀疏矩陣，具體到電力系統(tǒng)，其稠密度通常小于1%［2-3］，如果不能有效利用矩陣的稀疏性，在使用計(jì)算機(jī)處理大型稀疏矩陣時(shí)，大量的存儲(chǔ)和計(jì)算資源將會(huì)浪費(fèi)在無(wú)效的零元上，導(dǎo)致存儲(chǔ)空間不足和處理效率低下［4］。因此，在實(shí)現(xiàn)大型稀疏矩陣存儲(chǔ)和計(jì)算時(shí)，需要利用專用的算法和數(shù)據(jù)結(jié)構(gòu)，基于特殊設(shè)計(jì)的硬件架構(gòu)，最大化地利用計(jì)算系統(tǒng)的算力和效能［5-6］。綜上，電力系統(tǒng)的電磁暫態(tài)加速求解問(wèn)題最終聚焦于稀疏線性方程組的加速求解問(wèn)題。

稀疏線性方程組的求解方法包括直接法和迭代法兩大類。由于迭代法存在迭代次數(shù)不可控、結(jié)果精度較低等問(wèn)題，因此一般使用直接法進(jìn)行求解。直接法是指在不考慮計(jì)算舍入誤差的情況下，通過(guò)矩陣分解和三角方程［7-8］進(jìn)行求解。下三角稀疏矩陣求解是求解稀疏線性方程組的核心算法的組成部分［9］，而層次集合法［10-11］是最重要的并行性分析方法之一。近年來(lái)，越來(lái)越多的算法開始使用眾核加速硬件來(lái)實(shí)現(xiàn)并行加速。NAUMOV［12］提出基于層次集合的方法，其將多個(gè)小的層次集合進(jìn)行合并，以減少同步運(yùn)算。PARK等［13］通過(guò)對(duì)同步進(jìn)行剪枝來(lái)提升效率。LIU等［14-15］提出利用圖形處理器（Graphics Processing Unit，GPU）的原子操作來(lái)實(shí)現(xiàn)無(wú)同步的算法。SU等［16］介紹大規(guī)模的線程級(jí)無(wú)同步算法。LU等［17］介紹基于循環(huán)塊優(yōu)化的算法。上述算法主要在GPU 上進(jìn)行優(yōu)化，WANG等［18］則提出了在國(guó)產(chǎn)神威眾核上實(shí)現(xiàn)的加速算法。

相比眾核架構(gòu)（GPU 或神威眾核），現(xiàn)場(chǎng)可編程門陣列（Field Programmable Gate Array，F(xiàn)PGA）具有片上緩沖大、傳輸帶寬可定制、調(diào)度靈活等優(yōu)勢(shì)，可以最大化地利用計(jì)算系統(tǒng)的算力［19］。吳志勇等［20］設(shè)計(jì)了一種面向FPGA 求解稀疏矩陣的硬件架構(gòu)，本文基于這一硬件結(jié)構(gòu)，提出軟件映射和調(diào)度求解算法。該算法給出數(shù)據(jù)分布和指令排布的過(guò)程，通過(guò)將下三角稀疏矩陣的求解過(guò)程靜態(tài)映射到多個(gè)FPGA 片上的處理單元（Process Elements，PEs），以軟硬件協(xié)同的方法實(shí)現(xiàn)下三角稀疏矩陣在定制化FPGA 架構(gòu)上的高速求解。

1 FPGA 稀疏矩陣求解器硬件設(shè)計(jì)

FPGA 稀疏矩陣求解器的硬件設(shè)計(jì)在文獻(xiàn)［20］中已經(jīng)詳細(xì)說(shuō)明，本文只簡(jiǎn)述該求解器的基本原理和結(jié)構(gòu)，以便進(jìn)一步介紹本文所提軟件調(diào)度算法。直接法需要遵循下三角稀疏矩陣內(nèi)部的對(duì)角元依賴關(guān)系，其求解路徑構(gòu)成一個(gè)有向無(wú)環(huán)圖（Directed Acyclic Graph，DAG）［21］，如圖1 所示。由于DAG 的條件路徑依賴關(guān)系無(wú)法直接映射為并行化算法，因此無(wú)法直接在算法層進(jìn)行并行優(yōu)化。下三角稀疏矩陣具有稀疏性，構(gòu)成的DAG 存在隱式的數(shù)據(jù)并行，意味著多個(gè)對(duì)角元和非對(duì)角元之間存在并行處理的可能。整個(gè)系統(tǒng)的基本思想是設(shè)計(jì)多個(gè)硬件處理單元，每個(gè)處理單元都能夠有效地處理對(duì)角元和非對(duì)角元。通過(guò)軟件對(duì)該稀疏矩陣進(jìn)行預(yù)處理，找到稀疏下三角的并行性，并且將劃分完成的下三角稀疏矩陣映射到這些硬件處理單元上，從而實(shí)現(xiàn)軟硬件協(xié)同的并行求解。

圖1 有向無(wú)環(huán)圖Fig.1 Directed acyclic graph

由于軟件預(yù)處理及劃分完成后的稀疏下三角各個(gè)部分（對(duì)角元及非對(duì)角元之間）存在相互依賴，因此需要為多個(gè)處理單元間提供硬件連接通路。在實(shí)際的FPGA 實(shí)現(xiàn)過(guò)程中，使用二維單向環(huán)網(wǎng)來(lái)實(shí)現(xiàn)多個(gè)處理單元間的通信連接，從而完成各單元間的數(shù)據(jù)傳輸［22］。FPGA 稀疏矩陣求解器的硬件結(jié)構(gòu)如圖2 所示。

圖2 求解器的硬件結(jié)構(gòu)Fig.2 Hardware structure of solver

從圖2 可以看出，系統(tǒng)中存在多個(gè)處理單元，每個(gè)處理單元都有一個(gè)指令控制部件M，用來(lái)存儲(chǔ)PE內(nèi)部控制部件的指令碼。另外，每個(gè)處理單元內(nèi)還有3 類操作部件，分別是Mul 部件、Add 部件和T/R部件，其中：Mul 部件用來(lái)處理稀疏矩陣對(duì)角元和非對(duì)角元的雙精度浮點(diǎn)復(fù)數(shù)乘法（包括除法）操作；Add部件用來(lái)處理非對(duì)角元對(duì)同行對(duì)角元的更新；T/R部件用來(lái)實(shí)現(xiàn)相鄰PE 間的數(shù)據(jù)傳輸。

每個(gè)PE 內(nèi)部包括5 個(gè)專用緩沖，分別是A、A'、B、B'和B"，其中：緩沖A用來(lái)存儲(chǔ)稀疏矩陣對(duì)角元系數(shù)；緩沖A'用來(lái)存儲(chǔ)稀疏矩陣非對(duì)角元系數(shù)；緩沖B用來(lái)存儲(chǔ)稀疏矩陣右端項(xiàng)和求得的未知數(shù)；緩沖B'用來(lái)存儲(chǔ)非對(duì)角元和同列未知數(shù)的乘積；緩沖B"用來(lái)存儲(chǔ)其他PE 傳輸過(guò)來(lái)的已求解的未知數(shù)x。由于每個(gè)PE對(duì)B"緩沖的寫操作一定都指向不同的地址，因此B"緩沖可以在2 個(gè)相鄰單元間共享，以減少對(duì)FPGA 資源的占用，有利于FPGA 綜合實(shí)現(xiàn)。處理單元間則通過(guò)二維單向環(huán)網(wǎng)相連來(lái)進(jìn)行通信?；谠撚布軜?gòu)，原本串行的稀疏矩陣求解過(guò)程被分解為以下操作步驟：

1）使用Mul 部件，從緩沖A中讀取對(duì)角元系數(shù)，從緩沖B中讀取右端項(xiàng)，求解獲得未知數(shù)x，再存放到緩沖B中。這個(gè)x不被其他PE 上的非對(duì)角元需求，無(wú)需傳輸。

2）使用Mul 部件，從緩沖A中讀取對(duì)角元系數(shù)，從緩沖B中讀取右端項(xiàng)，求解獲得未知數(shù)x，存放到緩沖B中。如果這個(gè)x被其他PE 的非對(duì)角元需求，則同時(shí)將其發(fā)送到R部件。

3）使用Mul 部件，從緩沖A'中讀取非對(duì)角元系數(shù)，從緩沖B中讀取求得的x，計(jì)算獲得非對(duì)角元乘積，將其存放在緩沖B'中。

4）使用Mul 部件，從緩沖A'中讀取非對(duì)角元系數(shù)，從緩沖B"中讀取求得的x，計(jì)算獲得非對(duì)角元乘積，將其存放在緩沖B'中。本次讀取到的x對(duì)應(yīng)第2步操作中傳輸?shù)絉部件中的x。

5）使用Add 部件，從緩沖B'中讀取非對(duì)角元乘積，從緩沖B中讀取對(duì)角元，計(jì)算更新對(duì)角元，并將結(jié)果存放到緩沖B中。

6）使用T部件，將左方或者下方傳輸來(lái)的x通過(guò)網(wǎng)絡(luò)傳輸?shù)较乱粋€(gè)PE，并保存到R部件和緩沖B"中。

根據(jù)上述稀疏矩陣求解過(guò)程，雙精度浮點(diǎn)復(fù)數(shù)乘法單元和加法單元的數(shù)據(jù)調(diào)度過(guò)程分別如圖3和圖4所示。

圖3 乘法單元的數(shù)據(jù)調(diào)度Fig.3 Data scheduling of the multiplication unit

圖4 加法單元的數(shù)據(jù)調(diào)度Fig.4 Data scheduling of the addition unit

2 軟件靜態(tài)調(diào)度算法

基于上文FPGA 稀疏矩陣求解器的硬件架構(gòu)，本文實(shí)現(xiàn)了稀疏矩陣求解算法和求解過(guò)程的映射，充分挖掘了稀疏矩陣求解過(guò)程中的并行性，這一映射過(guò)程通過(guò)軟件靜態(tài)調(diào)度算法來(lái)實(shí)現(xiàn)。中央處理器（Central Processing Unit，CPU）一般通過(guò)順序求解所有解向量來(lái)實(shí)現(xiàn)，而GPU 則會(huì)對(duì)解向量進(jìn)行分塊，在分塊后的解向量之間通過(guò)引用計(jì)數(shù)構(gòu)建依賴關(guān)系，然后逐塊觸發(fā)從而完成求解。無(wú)論是CPU 還是GPU，都可以通過(guò)動(dòng)態(tài)尋址來(lái)定位解向量，但是，F(xiàn)PGA 不具備動(dòng)態(tài)尋址的能力，因此，必須構(gòu)建靜態(tài)指令流，指示每個(gè)操作部件在特定步時(shí)的訪存地址和處理操作。

2.1 稀疏矩陣預(yù)處理

在正式求解之前，需要對(duì)矩陣進(jìn)行預(yù)處理。稀疏矩陣存儲(chǔ)通常使用壓縮稀疏列矩陣（Compressed Sparse Column matrix，CSC）或壓縮稀疏行矩陣（Compressed Sparse Row matrix，CSR）［23］。首先，將矩陣進(jìn)行重排序，以減少LU 分解增加的額外非零元，本文使用metis［24］工具進(jìn)行重排序，其能提供一組可以獨(dú)立運(yùn)行的命令行程序，同時(shí)也提供應(yīng)用程序接口（Application Programming Interface，API），方便集成到C/C++或Fortran 程序中，該程序可以得到上述算法目標(biāo)的一個(gè)近似解；其次，對(duì)完成重排序的稀疏矩陣進(jìn)行LU 分解，獲得右下局部稠密的下三角稀疏矩陣，如果沒(méi)有特殊說(shuō)明，下文以L代指這里的下三角稀疏矩陣。

2.2 矩陣劃分

針對(duì)上述硬件架構(gòu)，軟件調(diào)度算法需要為其每個(gè)PE 進(jìn)行任務(wù)分派和調(diào)度，將求解下三角稀疏矩陣L需要用到的數(shù)據(jù)分布到多個(gè)計(jì)算單元上。

數(shù)據(jù)分布的設(shè)計(jì)目標(biāo)是使得盡可能多的計(jì)算單元盡可能飽和地運(yùn)行。第一個(gè)目標(biāo)是使得分布到多個(gè)PE 上的數(shù)據(jù)盡可能均勻，以保證每個(gè)PE 需要求解的數(shù)據(jù)量接近。由于稀疏矩陣的非零元計(jì)算過(guò)程實(shí)際上是一個(gè)有向無(wú)環(huán)圖，因此將一個(gè)稀疏矩陣映射到多個(gè)PE 上的問(wèn)題可以認(rèn)為是一個(gè)圖劃分（Graph Partition）問(wèn)題。第二個(gè)目標(biāo)是使得分割后的各個(gè)PE 間通信盡可能少，從而減小通信開銷，使計(jì)算單元盡可能飽和地運(yùn)行。

圖的均勻劃分是一個(gè)NP 困難（NP-hard）問(wèn)題，同樣使用metis 提供的接口，對(duì)L矩陣的對(duì)角元進(jìn)行劃分。劃分完成的對(duì)角元分布在不同的分塊中，這些分塊一一映射到多個(gè)不同的PE 上。原有的各個(gè)對(duì)角元之間的依賴關(guān)系被分成了PE 內(nèi)部的對(duì)角元依賴關(guān)系和跨PE 的對(duì)角元依賴關(guān)系，分別保存在pe_inside_diag_dep 和pe_outside_diag_dep 數(shù)組中。

非對(duì)角元被分布到其所在行的對(duì)角元所在的PE 上，這樣操作的好處是使得非對(duì)角元就緒后，可以直接更新其所在行的對(duì)角元，傳輸部件只需要傳輸對(duì)角元，能夠大幅簡(jiǎn)化與傳輸部件相關(guān)的算法設(shè)計(jì)。

2.3 傳輸算法

在實(shí)際的硬件實(shí)現(xiàn)中，當(dāng)PE 數(shù)較多時(shí)，PE 間直接相連會(huì)帶來(lái)巨大的網(wǎng)絡(luò)資源消耗，因此，硬件設(shè)計(jì)通常使用二維mesh 網(wǎng)絡(luò)來(lái)替代PE 間的直接連接。但是，二維mesh 網(wǎng)絡(luò)帶來(lái)的問(wèn)題是每個(gè)PE 只和周圍2 個(gè)維度上的4 個(gè)PE 直接相連，當(dāng)要和遠(yuǎn)端的PE通信時(shí)，需要經(jīng)過(guò)中間的PE 進(jìn)行連接，此時(shí)需要多個(gè)時(shí)鐘周期來(lái)完成傳輸。

當(dāng)前某個(gè)PE 計(jì)算獲得的x有可能會(huì)被最遠(yuǎn)端的PE 訪問(wèn)，因此，在軟件調(diào)度算法上，一開始選擇廣播算法進(jìn)行數(shù)據(jù)共享。廣播算法可以保證當(dāng)前PE 計(jì)算求得的未知數(shù)x最終一定能夠被所有的PE 訪問(wèn)到，并且對(duì)于二維單向環(huán)網(wǎng)而言，廣播算法非常容易設(shè)計(jì)和實(shí)現(xiàn)，如圖5 所示。二維單向環(huán)網(wǎng)中所有的PE 對(duì)于網(wǎng)絡(luò)是完全對(duì)稱的，從非0 號(hào)PE 發(fā)起的廣播傳輸過(guò)程也是類似的。在使用廣播算法時(shí)，同一時(shí)刻多個(gè)PE 發(fā)起的消息只要滿足PE 的行號(hào)與列號(hào)之和不同于PE 數(shù)量（mod PE），就可以在系統(tǒng)中并行傳輸。

圖5 從0 號(hào)PE 開始的16PE 廣播算法Fig.5 16PE broadcast algorithm starting from PE 0

實(shí)際上，并非所有的PE 都需要獲得當(dāng)前正在傳輸?shù)膞，廣播算法會(huì)占用不需要x的PE 的傳輸部件T/R。為了解決這一問(wèn)題，可以將廣播算法改為點(diǎn)對(duì)點(diǎn)傳輸算法，點(diǎn)對(duì)點(diǎn)算法的路由與廣播算法一致，但只占用發(fā)送PE 和接收PE 之間的通路，空出的通路可以同時(shí)傳輸其他x，這使得系統(tǒng)中能夠同時(shí)傳輸更多的x，進(jìn)一步提升了傳輸效率。

2.4 指令排布和生成

在數(shù)據(jù)分布完成后，就可以對(duì)指令進(jìn)行排布，需要完整地規(guī)劃每個(gè)PE 上的部件操作，從而充分利用每個(gè)PE 的操作能力。在指令排布時(shí)需要確定以下信息：

1）需要明確硬件時(shí)鐘的頻率和延遲。

上述硬件設(shè)計(jì)中的操作部件并不是立刻得出結(jié)果的，當(dāng)頻率不同時(shí)，各個(gè)操作部件可能會(huì)產(chǎn)生一定的延遲。在一般情況下，頻率越高，延遲越大。當(dāng)頻率控制在300～400 MHz 時(shí)，Mul 部件的延遲是5 拍，而Add 部件的延遲是3 拍。在指令排布時(shí)，設(shè)定Mul_lat 和Add_lat 分別表示Mul 部件和Add 部件的延遲。

2）需要明確系統(tǒng)中各個(gè)部件間的依賴關(guān)系。

如上文所述，系統(tǒng)中主要存在2 種依賴關(guān)系，即PE 內(nèi)部對(duì)角元之間的依賴和跨PE 的對(duì)角元之間的依賴，分別被保存在pe_inside_diag_dep和pe_outside_diag_dep 數(shù)組中，這里直接使用這2 個(gè)數(shù)組即可。

3）需要明確系統(tǒng)中各個(gè)部件間的沖突關(guān)系。

通過(guò)分析可以發(fā)現(xiàn)，當(dāng)前硬件系統(tǒng)中沒(méi)有訪存沖突，所有的內(nèi)存模塊至多只有2 個(gè)寫入端口和2 個(gè)讀出端口，可以通過(guò)分頻實(shí)現(xiàn)復(fù)用。然而，從上述求解過(guò)程可以看出，系統(tǒng)中存在不同操作對(duì)部件的爭(zhēng)用沖突，主要包括以下3 種：

（1）對(duì)Mul 部件的爭(zhēng)用。當(dāng)Mul 部件用于對(duì)角元求解時(shí)，就無(wú)法用于非對(duì)角元的乘法，反之亦然。

（2）對(duì)R部件的沖突。如果T部件在某一時(shí)刻被占用，那么相應(yīng)的R部件也被占用，無(wú)法再進(jìn)行傳輸。

（3）對(duì)Add 部件的爭(zhēng)用。當(dāng)Add 部件正在處理某一行的右端項(xiàng)更新時(shí)，同行的非對(duì)角元不能同時(shí)更新右端項(xiàng)，否則會(huì)造成讀寫沖突。

為了解決沖突問(wèn)題，需要為上述可能發(fā)生沖突的部件設(shè)置部件占用標(biāo)記。某一時(shí)刻，當(dāng)部件占用標(biāo)記被置位時(shí)，意味著某個(gè)部件已經(jīng)被占用，無(wú)法再用其排布指令。本文將上述3 個(gè)部件占用標(biāo)記分別設(shè)置為Mul_occupy、R_occupy 和Add_occupy，其中，Add_occupy 是一個(gè)長(zhǎng)度為Add_lat 的列表，用來(lái)標(biāo)記當(dāng)前加法部件正在處理的行。

圖6 所示為指令排布的算法流程，算法的偽代碼如算法1 所示。在排布時(shí)，乘法單元會(huì)優(yōu)先處理非對(duì)角元，只有當(dāng)非對(duì)角元全部處理完成后才會(huì)處理對(duì)角元，這樣做是為了盡可能多地在PE 上生成可供處理的對(duì)角元，以便保持PE 的滿載，提高整個(gè)系統(tǒng)的并行度。

圖6 指令排布算法流程Fig.6 The procedure of instruction layout algorithm

算法1指令排布算法

指令排布完成后就可以生成指令。系統(tǒng)中主要有6 種不同的操作，其中，Mul 部件占據(jù)了4 種操作，Add 部件和T部件各占據(jù)1 種操作，每個(gè)時(shí)刻都有可能出現(xiàn)所有操作部件同時(shí)工作的情況。由于各操作的取數(shù)空間基址和偏移都不同，需要將多種不同的操作整合起來(lái)生成指令。生成的指令放在緩沖M中，每個(gè)存儲(chǔ)單元的位寬為128 位，指令格式中的每個(gè)字段的位、名稱和含義如表1 所示。

表1 指令及其含義Table 1 Instructions and their meanings

2.5 硬件空間限制

在硬件實(shí)現(xiàn)的過(guò)程中，通常都存在一定的空間限制。由于緩沖區(qū)A、A'、B、B'、B"都使用片上內(nèi)存實(shí)現(xiàn)，導(dǎo)致其空間相對(duì)受限，也意味著求解的矩陣大小是受限的。

緩沖區(qū)A和B分別存儲(chǔ)分布在某個(gè)PE 上的對(duì)角元系數(shù)和右端項(xiàng)（包括求解得到的x）。假設(shè)對(duì)圖的分割相對(duì)均勻，則A和B的大小約為（N/PE_num），其中，N表示矩陣階數(shù)，PE_num 表示PE 的個(gè)數(shù)。

A'存儲(chǔ)非對(duì)角元系數(shù)，B'暫存乘以同列未知數(shù)x后的非對(duì)角元。非對(duì)角元不作劃分，而是直接分配到所在行對(duì)角元的PE 上。假設(shè)分布均勻，則A'和B'的大小約為（nnz/PE_num)，其中，nnz 表示非零元的個(gè)數(shù)。

B"存儲(chǔ)所有PE 間求解得到并且需要傳輸?shù)膞。B''的大小受矩陣稀疏度和PE_num 影響，矩陣越稠密，PE 的個(gè)數(shù)越多（劃分越多），B''則越大。B"最大不會(huì)超過(guò)N，實(shí)際大小依賴于具體矩陣，一般來(lái)說(shuō)會(huì)比N小得多。

硬件總的空間占用約為(N+nnz+N×PE_num）個(gè)實(shí)數(shù)或復(fù)數(shù)，在使用算法前，需要評(píng)估FPGA 上提供的空間是否能滿足實(shí)際的空間需求。

3 性能測(cè)試分析

算法求解的過(guò)程可以分為3 個(gè)階段：

1）第一階段是矩陣求解初期的稀疏部分，此時(shí)絕大部分PE（大于80%）能夠找到就緒的對(duì)角元或非對(duì)角元，用于運(yùn)行乘法部件，即使有依賴元暫未求解或硬件資源存在沖突，通常也只需等待1～2 拍就能夠繼續(xù)運(yùn)行。

2）第二階段是矩陣求解中期的相對(duì)稀疏部分，此時(shí)只有部分PE（大于10%而小于80%）能夠運(yùn)行乘法部件，剩余的PE 由于對(duì)角元未就緒（需要等待同行非對(duì)角元處理完成）或找不到非對(duì)角元（需要等待求解出的未知數(shù)x）進(jìn)行處理，因此只能等待。

3）第三階段是矩陣求解后期的稠密部分，此時(shí)只有極少部分PE（小于10%）能夠運(yùn)行乘法部件，多個(gè)未知數(shù)x求解時(shí)存在強(qiáng)相關(guān)性，求解過(guò)程串行化明顯增加。

圖7 所示為FPGA 加速器的邏輯結(jié)構(gòu)。本文對(duì)應(yīng)的項(xiàng)目在Xilinx XCVU37P（8 GB HBM、4 200 萬(wàn)等效邏輯門）FPGA 上實(shí)現(xiàn)硬件加速器，其與主機(jī)間采用PCIE4.0 X16 接口（帶寬為64 GB/s）連接。主機(jī)采用20 核Intel Xeon Gold 5320H（4.3 GHz）處理器，128 GB 3 200 Mb/s DDR4配置，運(yùn) 行Linux 操作系統(tǒng)，移植和適配了電力系統(tǒng)仿真所需的庫(kù)環(huán)境。FPGA 加速器實(shí)物如圖8 所示。本文硬件結(jié)構(gòu)設(shè)計(jì)開銷情況如表2 所示。

圖7 FPGA 加速器邏輯結(jié)構(gòu)Fig.7 Logic structure of FPGA accelerator

圖8 FPGA 加速器實(shí)物示意圖Fig.8 Physical diagram of FPGA accelerator

表2 硬件資源利用情況Table 2 Utilization of hardware resource

本文對(duì)2 個(gè)典型算例進(jìn)行測(cè)試，2 個(gè)算例均來(lái)源于實(shí)際的電網(wǎng)模型。算例1 的矩陣大小為10 188×10 188，非零元為25 720 個(gè)；算例2 的矩陣大小為21 464×21 464，非零元為121 890 個(gè)。算例測(cè)試結(jié)果如表3 所示。

表3 2 個(gè)典型算例測(cè)試結(jié)果Table 3 Test results of two typical examples

從表3 可以看出，在64 個(gè)PE 配置的FPGA 加速器中，算例1 約耗時(shí)1 251 拍，則每拍能夠處理約20 個(gè)非零元，算例2 約耗時(shí)4 672 拍，每拍能夠處理約25 個(gè)非零元。

將當(dāng)前國(guó)網(wǎng)電力系統(tǒng)電磁暫態(tài)仿真程序中計(jì)算最為密集的稀疏矩陣消元求解核心段作為實(shí)際測(cè)試對(duì)象，對(duì)基于FPGA 和基于傳統(tǒng)CPU/GPU 的2 種環(huán)境進(jìn)行對(duì)比測(cè)試和分析。當(dāng)FPGA 加速器系統(tǒng)在300 MHz 頻率、256 個(gè)PE 配置下時(shí)，針對(duì)上述核心段的處理效率能達(dá)到30.84 GFLOPs。與此同時(shí)，將該電磁暫態(tài)仿真程序移植到20 核Intel Xeon Gold 5320H（4.3 GHz）處理器平臺(tái)上，進(jìn)行多核上的多線程并行優(yōu)化，同時(shí)降低大數(shù)據(jù)量離散訪問(wèn)導(dǎo)致的CACHE 緩存顛簸，減少訪存延遲。處理器中4 核用于操作系統(tǒng)和驅(qū)動(dòng)運(yùn)行，其余16 核用于優(yōu)化后的稀疏矩陣并行運(yùn)算，其處理效率僅為5.22 GFLOPs?；贔PGA 算法的實(shí)測(cè)性能是傳統(tǒng)CPU/GPU 求解算法的5.9 倍，加速效果顯著。

雖然本文所提算法相較傳統(tǒng)CPU/GPU 算法有明顯的加速效果，但仍然存在優(yōu)化的空間。從上述3 個(gè)求解階段來(lái)看，第一階段已經(jīng)充分利用了所有的硬件計(jì)算資源，幾乎沒(méi)有可以優(yōu)化的余地；第二階段的節(jié)拍占比最大，隨著矩陣規(guī)模的增大，其增長(zhǎng)速度也最快，存在優(yōu)化算法調(diào)度的可能，主要優(yōu)化思路是更好地進(jìn)行矩陣的劃分和映射，進(jìn)一步減小通信代價(jià)，提升算法并行度；第三階段主要是矩陣的稠密部分，可以使用逆矩陣乘法進(jìn)行優(yōu)化。在LU 分解時(shí)，通過(guò)行列調(diào)整可以使得矩陣L右下角局部稠密，對(duì)于稠密的部分，可以求其逆矩陣，通過(guò)逆矩陣乘法進(jìn)行求解。由于矩陣乘法可以在多個(gè)PE 上完全并行，因此當(dāng)矩陣有一定的稠密度時(shí)，用逆矩陣乘法進(jìn)行求解相比直接法更具性能優(yōu)勢(shì)。

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

本文提出一種基于FPGA 硬件的靜態(tài)調(diào)度優(yōu)化算法，利用該算法實(shí)現(xiàn)了下三角稀疏矩陣的求解。通過(guò)對(duì)稀疏矩陣直接法求解步驟的分解和對(duì)稀疏矩陣的解析排布，設(shè)計(jì)一種節(jié)拍級(jí)的靜態(tài)調(diào)度流程，以充分利用FPGA 的硬件資源獲取較高的求解效率和硬件利用率。性能測(cè)試結(jié)果驗(yàn)證了該算法的高效性，對(duì)于在FPGA 上實(shí)現(xiàn)類似的基于圖劃分的隱式數(shù)據(jù)并行算法具有一定的參考意義。下一步擬對(duì)LU 分解中的右下角局部稠密矩陣進(jìn)行優(yōu)化，無(wú)需修改現(xiàn)有的硬件拓?fù)?，僅增加稠密懲罰操作指令，在軟件方面，對(duì)稠密部分進(jìn)行求逆并均勻劃分求得的逆矩陣，然后通過(guò)乘法來(lái)求解該矩陣。