FPGA 架構上面向稀疏矩陣求解的靜態調度算法

王晞陽，陳繼林，李猛，劉首文

（1.國家超級計算無錫中心，江蘇 無錫 214072；2.中國電力科學研究院有限公司，北京 100192；3.國網湖北省電力有限公司，武漢 430070）

0 概述

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

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

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

1 FPGA 稀疏矩陣求解器硬件設計

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

圖1 有向無環圖Fig.1 Directed acyclic graph

由于軟件預處理及劃分完成后的稀疏下三角各個部分（對角元及非對角元之間）存在相互依賴，因此需要為多個處理單元間提供硬件連接通路。在實際的FPGA 實現過程中，使用二維單向環網來實現多個處理單元間的通信連接，從而完成各單元間的數據傳輸［22］。FPGA 稀疏矩陣求解器的硬件結構如圖2 所示。

圖2 求解器的硬件結構Fig.2 Hardware structure of solver

從圖2 可以看出，系統中存在多個處理單元，每個處理單元都有一個指令控制部件M，用來存儲PE內部控制部件的指令碼。另外，每個處理單元內還有3 類操作部件，分別是Mul 部件、Add 部件和T/R部件，其中：Mul 部件用來處理稀疏矩陣對角元和非對角元的雙精度浮點復數乘法（包括除法）操作；Add部件用來處理非對角元對同行對角元的更新；T/R部件用來實現相鄰PE 間的數據傳輸。

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

1）使用Mul 部件，從緩沖A中讀取對角元系數，從緩沖B中讀取右端項，求解獲得未知數x，再存放到緩沖B中。這個x不被其他PE 上的非對角元需求，無需傳輸。

2）使用Mul 部件，從緩沖A中讀取對角元系數，從緩沖B中讀取右端項，求解獲得未知數x，存放到緩沖B中。如果這個x被其他PE 的非對角元需求，則同時將其發送到R部件。

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

4）使用Mul 部件，從緩沖A'中讀取非對角元系數，從緩沖B"中讀取求得的x，計算獲得非對角元乘積，將其存放在緩沖B'中。本次讀取到的x對應第2步操作中傳輸到R部件中的x。

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

6）使用T部件，將左方或者下方傳輸來的x通過網絡傳輸到下一個PE，并保存到R部件和緩沖B"中。

根據上述稀疏矩陣求解過程，雙精度浮點復數乘法單元和加法單元的數據調度過程分別如圖3和圖4所示。

圖3 乘法單元的數據調度Fig.3 Data scheduling of the multiplication unit

圖4 加法單元的數據調度Fig.4 Data scheduling of the addition unit

2 軟件靜態調度算法

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

2.1 稀疏矩陣預處理

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

2.2 矩陣劃分

針對上述硬件架構，軟件調度算法需要為其每個PE 進行任務分派和調度，將求解下三角稀疏矩陣L需要用到的數據分布到多個計算單元上。

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

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

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

2.3 傳輸算法

在實際的硬件實現中，當PE 數較多時，PE 間直接相連會帶來巨大的網絡資源消耗，因此，硬件設計通常使用二維mesh 網絡來替代PE 間的直接連接。但是，二維mesh 網絡帶來的問題是每個PE 只和周圍2 個維度上的4 個PE 直接相連，當要和遠端的PE通信時，需要經過中間的PE 進行連接，此時需要多個時鐘周期來完成傳輸。

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

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

實際上，并非所有的PE 都需要獲得當前正在傳輸的x，廣播算法會占用不需要x的PE 的傳輸部件T/R。為了解決這一問題，可以將廣播算法改為點對點傳輸算法，點對點算法的路由與廣播算法一致，但只占用發送PE 和接收PE 之間的通路，空出的通路可以同時傳輸其他x，這使得系統中能夠同時傳輸更多的x，進一步提升了傳輸效率。

2.4 指令排布和生成

在數據分布完成后，就可以對指令進行排布，需要完整地規劃每個PE 上的部件操作，從而充分利用每個PE 的操作能力。在指令排布時需要確定以下信息：

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

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

2）需要明確系統中各個部件間的依賴關系。

如上文所述，系統中主要存在2 種依賴關系，即PE 內部對角元之間的依賴和跨PE 的對角元之間的依賴，分別被保存在pe_inside_diag_dep和pe_outside_diag_dep 數組中，這里直接使用這2 個數組即可。

3）需要明確系統中各個部件間的沖突關系。

通過分析可以發現，當前硬件系統中沒有訪存沖突，所有的內存模塊至多只有2 個寫入端口和2 個讀出端口，可以通過分頻實現復用。然而，從上述求解過程可以看出，系統中存在不同操作對部件的爭用沖突，主要包括以下3 種：

（1）對Mul 部件的爭用。當Mul 部件用于對角元求解時，就無法用于非對角元的乘法，反之亦然。

（2）對R部件的沖突。如果T部件在某一時刻被占用，那么相應的R部件也被占用，無法再進行傳輸。

（3）對Add 部件的爭用。當Add 部件正在處理某一行的右端項更新時，同行的非對角元不能同時更新右端項，否則會造成讀寫沖突。

為了解決沖突問題，需要為上述可能發生沖突的部件設置部件占用標記。某一時刻，當部件占用標記被置位時，意味著某個部件已經被占用，無法再用其排布指令。本文將上述3 個部件占用標記分別設置為Mul_occupy、R_occupy 和Add_occupy，其中，Add_occupy 是一個長度為Add_lat 的列表，用來標記當前加法部件正在處理的行。

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

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

算法1指令排布算法

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

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

2.5 硬件空間限制

在硬件實現的過程中，通常都存在一定的空間限制。由于緩沖區A、A'、B、B'、B"都使用片上內存實現，導致其空間相對受限，也意味著求解的矩陣大小是受限的。

緩沖區A和B分別存儲分布在某個PE 上的對角元系數和右端項（包括求解得到的x）。假設對圖的分割相對均勻，則A和B的大小約為（N/PE_num），其中，N表示矩陣階數，PE_num 表示PE 的個數。

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

B"存儲所有PE 間求解得到并且需要傳輸的x。B''的大小受矩陣稀疏度和PE_num 影響，矩陣越稠密，PE 的個數越多（劃分越多），B''則越大。B"最大不會超過N，實際大小依賴于具體矩陣，一般來說會比N小得多。

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

3 性能測試分析

算法求解的過程可以分為3 個階段：

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

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

3）第三階段是矩陣求解后期的稠密部分，此時只有極少部分PE（小于10%）能夠運行乘法部件，多個未知數x求解時存在強相關性，求解過程串行化明顯增加。

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

圖7 FPGA 加速器邏輯結構Fig.7 Logic structure of FPGA accelerator

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

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

本文對2 個典型算例進行測試，2 個算例均來源于實際的電網模型。算例1 的矩陣大小為10 188×10 188，非零元為25 720 個；算例2 的矩陣大小為21 464×21 464，非零元為121 890 個。算例測試結果如表3 所示。

表3 2 個典型算例測試結果Table 3 Test results of two typical examples

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

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

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

4 結束語

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