一、背景

    在工業(yè)仿真領(lǐng)域，對各種現(xiàn)實世界的問題進行數(shù)值模擬時，如流體動力學(xué)分析、電磁場仿真、結(jié)構(gòu)力學(xué)應(yīng)力應(yīng)變分析等，其控制方程通常是偏微分方程組，在經(jīng)過不同方法的隱式離散之后最終都可轉(zhuǎn)化為大型稀疏線性方程組。隨著人們對計算精度要求的不斷提高，方程組的階數(shù)也從上千階、幾十萬階提高到百萬、千萬階甚至更高，所需的計算量以及存儲需求也隨之迅速膨脹。根據(jù)一般經(jīng)驗，方程組求解時間會占總計算時間的70%以上，往往是整個計算過程中的性能瓶頸。如果說求解器是工業(yè)CAE軟件的核心模塊，那么大型稀疏線性方程組的求解技術(shù)將毫無疑問是底層求解器的核心。   

   

NASA翼型網(wǎng)格經(jīng)過離散得到的稀疏矩陣（素材來源于網(wǎng)絡(luò)）    

二、方法

    眾所周知，稀疏線性方程組的求解方法可以分為直接法和迭代法 ，兩類方法各有優(yōu)劣，特點比較如下：

    迭代法[1]：   

1. 對于不同類型稀疏矩陣表現(xiàn)差異較大，存在收斂性與收斂速度問題，催生了許多預(yù)處理技術(shù)（Preconditioners）；
2. 對原矩陣的編輯很少，SpMV(Sparse matrix-vector multiplication)是其核心運算；
3. 內(nèi)存需求小，求解速度較快，算法復(fù)雜度低；
4. 較易實現(xiàn)并行化。

    直接法[2]：   

1. 通用、穩(wěn)定；通過前后處理，能夠保證計算的收斂性與精度；
2. 對原矩陣的編輯多（分解、排序、縮放等）；
3. 內(nèi)存需求大，求解速度慢，算法復(fù)雜度更高；
4. 并行度有限。

   其中迭代法的種類很多，可以分為定常（Stationary）迭代法與非定常迭代法[3]。經(jīng)典的定常迭代法有Jacobi迭代法、Gauss-Seidel迭代法、SOR迭代法等，這些方法均可基于矩陣分裂推導(dǎo)得到；而在數(shù)值模擬中，非定常迭代法則顯得更加重要，常見的有共軛梯度法（Conjugate Gradient CG）、廣義最小殘量方法（Generalized Minimal Residual GMRES）、穩(wěn)定雙共軛梯度法（Biconjugate Gradient Stabilized Bi-CGSTAB）等，這些算法都屬于Krylov子空間方法。 如今高性能共軛梯度（HPCG）測試包已成為國際上評測超級計算機性能的主要工具[4]，Krylov子空間方法也被評為20世紀(jì)最偉大的十大算法之一[5]。 還有一些迭代算法則比較特殊，它們基于特殊結(jié)構(gòu)和性質(zhì)，例如我們耳熟能詳?shù)亩嘀鼐W(wǎng)格（Multigrid MG）方法，它更多地被稱為數(shù)值計算領(lǐng)域中一種加速迭代收斂的技術(shù)，而不僅僅是一種單純的算法。   

  

  

2014-2021年HPCG性能評測結(jié)果對比（素材來源于網(wǎng)絡(luò)）   

  

    直接法的基礎(chǔ)是矩陣的分解，常見的分解形式有LU分解、Cholesky分解、QR分解等。稀疏線性方程組的兩類常見直接求解算法分別為超節(jié)點（Supernodal）方法和多波前（Multifrontal）法，其主要思想是將完整的稀疏矩陣的分解任務(wù)轉(zhuǎn)化成許多個相對稠密的子矩陣的分解任務(wù)，任務(wù)間的依賴關(guān)系由消去樹（Elimination tree）或其他類似的數(shù)據(jù)結(jié)構(gòu)來確定。直接法的求解步驟通常分為矩陣重排、符號分解、數(shù)值分解與回代求解四個部分。   

一個稀疏矩陣與其對應(yīng)的消去樹（來自文獻6）    

三、挑戰(zhàn)

    當(dāng)前，國產(chǎn)超級計算機的峰值性能已達每秒十億億次量級，不久便將進入百億億次（E級）計算時代，我國的神威E級計算機和天河E級計算機已經(jīng)蓄勢待發(fā)。這些國際領(lǐng)先的超級計算機為我國科學(xué)與工程計算應(yīng)用邁進超大規(guī)模計算時代、實現(xiàn)更高精細(xì)度的數(shù)值模擬提供了強力支撐。然而，超大規(guī)模計算也給高實用性與高性能的大型稀疏線性方程組求解的算法設(shè)計與優(yōu)化帶來了巨大挑戰(zhàn)。

1.高實用

   當(dāng)前，國產(chǎn)超級計算機的峰值性能已達每秒十億億次量級，不久便將進入百億億次（E級）計算時代，我國的神威E級計算機和天河E級計算機已經(jīng)蓄勢待發(fā)。這些國際領(lǐng)先的超級計算機為我國科學(xué)與工程計算應(yīng)用邁進超大規(guī)模計算時代、實現(xiàn)更高精細(xì)度的數(shù)值模擬提供了強力支撐。然而，超大規(guī)模計算也給高實用性與高性能的大型稀疏線性方程組求解的算法設(shè)計與優(yōu)化帶來了巨大挑戰(zhàn)。    

以SiP封裝芯片的電磁-熱-力耦合數(shù)值模擬為例，其稀疏矩陣具有明顯的病態(tài)特征（來自文獻7）

 

2.高性能

       隨著計算機硬件性能的提升，超級計算機呈現(xiàn) “多級嵌套并行、異構(gòu)眾核加速” 的復(fù)雜體系結(jié)構(gòu)特征，會導(dǎo)致大型稀疏線性方程組求解器的實現(xiàn)效率急劇下降。從下圖也可以發(fā)現(xiàn)，隨著高性能計算機系統(tǒng)變得更加復(fù)雜，特別是眾核架構(gòu)采用后，在每一代世界性能最為強大的超級計算機上，應(yīng)用的實際求解能力變得更加低效，即解決問題時間（time-to-solution）變得越來越長。如何設(shè)計能 匹配機器體系結(jié)構(gòu)特征 的算法與性能優(yōu)化技術(shù)，是大型稀疏線性方程組求解技術(shù)以及其他科學(xué)計算核心算法中當(dāng)前亟待解決的關(guān)鍵問題。  

      

解決問題時間與超級計算機性能趨勢對比

    對于大規(guī)模稀疏線性方程組，原有串行和小規(guī)模并行模式下的數(shù)據(jù)結(jié)構(gòu)和算法容易導(dǎo)致并行求解性能低下或失敗。在分布式并行層面，需要解決以下幾個問題：一是在數(shù)據(jù)和任務(wù)分解方面，如何設(shè)計良好的負(fù)載均衡策略、稀疏矩陣的高效存儲格式以及計算通信重疊等優(yōu)化策略；二是在負(fù)載均衡的前提下，如何設(shè)計以盡力避免節(jié)點間的通信；三是在內(nèi)在串行特性導(dǎo)致并行化困難的算法方面，如何改進數(shù)據(jù)的分布方式以增加并行性。

    在共享內(nèi)存環(huán)境中，稀疏線性方程組求解算法的可擴展性問題也需要特別關(guān)注。因為現(xiàn)代多核/眾核處理器上的核數(shù)在可預(yù)見的未來也將越多，在單個CPU上封裝數(shù)十甚至上百個有較強處理能力的核心，或是在GPU上封裝成千上萬個輕量級處理單元將變得非常普遍。如何在這種共享內(nèi)存節(jié)點上實現(xiàn)細(xì)粒度的并行仍然是很有挑戰(zhàn)的研究內(nèi)容。  

    

  左圖：AMD霄龍CPU，64核128線程， 右圖：英偉達Hopper架構(gòu)GPU，1.8萬核心  

四、我們的探索——UNAP

    面對來自“應(yīng)用與機器”的雙重挑戰(zhàn)，神工坊團隊與國產(chǎn)異構(gòu)眾核平臺體系架構(gòu)緊密耦合，發(fā)展了一套大型稀疏線性方程組求解庫UNAP。該求解庫是早期應(yīng)基于非結(jié)構(gòu)網(wǎng)格的仿真需求而開發(fā)的，其全稱為“UNstructured Algebra Package”。

    為了解決高實用性的挑戰(zhàn)，UNAP結(jié)合異構(gòu)眾核處理器多級并行的特點和稀疏矩陣迭代解法的需求，初步探索了各種預(yù)處理方法在眾核異構(gòu)平臺上的并行實現(xiàn)技術(shù)。UNAP已實現(xiàn)PCG、PBiCGStab、GMRES等Krylov子空間方法以及AMG代數(shù)多重網(wǎng)格求解器，預(yù)處理器包含Jacobi、DIC、DILU等，未來還將開發(fā)直接求解模塊，以滿足來自各領(lǐng)域的復(fù)雜應(yīng)用需求。

    為了達到高性能的目標(biāo)，UNAP根據(jù)國產(chǎn)超算的異構(gòu)特點，結(jié)合其處理器的多級緩沖區(qū)，實現(xiàn)了計算/通信混合的并行迭代算法；由于迭代算法的并發(fā)度天然較高，UNAP主要通過算法調(diào)整、浮點運算替換內(nèi)存訪問、通信同步等手段[9]，充分利用了申威主從眾核浮點計算性能高而帶寬受限的特點，且通過增加計算比例降低了全局集合通信代價。在共享內(nèi)存層面，UNAP還可以調(diào)用太湖之光超級計算機上的加速工具套件UNAT和向量計算加速庫swArrays，充分發(fā)揮出從核陣列的計算能力，達到進一步的性能提升。  

    

代數(shù)求解庫UNAP的組織架構(gòu)  

    UNAP的整體結(jié)構(gòu)如上圖所示， 底層并行環(huán)境主要包含對“神威·太湖之光”超級計算機申威眾核異構(gòu)芯片SW26010的從核并行計算支持和普通的進程間MPI層級并行支持。采用以容器模板為核心的架構(gòu)體系， 基本容器層主要包含Vector分布式向量容器以及Matrix分布式矩陣容器。在容器層向下可調(diào)用向量計算加速庫和非結(jié)構(gòu)計算加速庫，實現(xiàn)在國產(chǎn)申威眾核芯片上的高效并行，同時，向上可擴展預(yù)條件子、迭代解法以及直接解法。這種架構(gòu)可以將元素的內(nèi)部處理通過模板泛化，從而降低求解算法實現(xiàn)的復(fù)雜度，也具有較好的 可移植性 ，易于兼容x86-CPU平臺。核心的代數(shù)求解器模塊包含了Krylov子空間迭代求解器、預(yù)條件子、代數(shù)多重網(wǎng)格求解器以及開發(fā)中的直接求解模塊。  

幾個案例說明UNAP的應(yīng)用情況：

•  發(fā)動機燃燒室大渦模擬： 在某航空發(fā)動機全環(huán)燃燒室的大渦模擬中，網(wǎng)格量達10億，采用UNAP作為核心求解模塊后，最終的并行規(guī)模達到1.6萬進程，稀疏線性方程組求解部分加速達到20倍，收斂速度顯著提升。
•  船舶水動力學(xué)應(yīng)用： 在某水動力學(xué)軟件中，通過UNAP的接入（主要使用了求解壓力和壓力修正方程的   代數(shù)多重網(wǎng)格算法 和求解速度等方程 的預(yù)條件穩(wěn)定雙共軛梯度法 ），獲得了在代數(shù)求解過程中的自動多級并行能力，在神威·太湖之光超級計算機上完成了千萬級網(wǎng)格對標(biāo)算例計算， 并行規(guī)模達到萬核級別，相對百進程 并行效率不低于50% ，經(jīng)測試與商業(yè)軟件FLUENT相當(dāng)。
•  電機設(shè)備電磁場分析： 在某公司核心求解器向神威平臺的移植部署中，采用UNAP代替原有直接求解庫，進行了網(wǎng)格量為千萬級的電機模型有限元算例并行計算測試，并對比了替換前后算例的節(jié)點磁通密度計算結(jié)果，UNAP表現(xiàn)良好。

一個典型算例展示UNAP的功能與性能：     

•  算例 ： 方腔驅(qū)動流，不可壓，上壁面滑移速度為1，其他為固壁邊界，方腔的邊長為1，雷諾數(shù)100。
•  并行計算結(jié)果： 以4進程為例，下圖給出了方腔在各個進程的分割情況，流線圖中可以清楚看到主渦和二次渦。

   

方腔驅(qū)動流的流線圖    

   

簡單的并行計算效率測試：            

•  強擴展性測試：

    固定算例網(wǎng)格整體求解規(guī)模，通過比較不同并行核數(shù)下的計算時間可獲得 強擴展性 并行計算效率。本測試采用的網(wǎng)格規(guī)模數(shù)量為2千萬，網(wǎng)格類型為六面體，進程數(shù)量分別為64、256、1024，并行核數(shù)規(guī)模分別為4160、16640和66560，統(tǒng)計時間為前20步計算時間。      

•  弱擴展性測試：

    固定每個進程上的網(wǎng)格規(guī)模，通過比較不同并行核數(shù)下的計算時間來獲得 弱擴展性 并行計算效率。本測試采用的網(wǎng)格規(guī)模數(shù)量分別為2千萬、4千萬、8千萬，網(wǎng)格類型為六面體，進程數(shù)量分別為64進程、128進程和256進程，并行核數(shù)規(guī)模分別為4160、16640和66560，單個進程網(wǎng)格數(shù)量約為31.25萬，統(tǒng)計時間為前20步計算時間。（本文作者：趙程鵬）

    

參考文獻：

[1] Saad Y. Iterative methods for sparse linear systems[M]. Society for Industrial and Applied Mathematics, 2003.

[2] Davis T A, Rajamanickam S, Sid-Lakhdar W M. A survey of direct methods for sparse linear systems[J]. Acta Numerica, 2016, 25: 383-566.

[3] Barrett R, Berry M, Chan T F, et al. Templates for the solution of linear systems: building blocks for iterative methods[M]. Society for Industrial and Applied Mathematics, 1994.

[4] Marjanovi? V, Gracia J, Glass C W. Performance modeling of the HPCG benchmark[C]//International Workshop on Performance Modeling, Benchmarking and Simulation of High Performance Computer Systems. Springer, Cham, 2014: 172-192.

[5] Cipra B A. The best of the 20th century: Editors name top 10 algorithms[J]. SIAM news, 2000, 33(4): 1-2.

[6] Gupta A, Karypis G, Kumar V. Highly scalable parallel algorithms for sparse matrix factorization[J]. IEEE Transactions on Parallel and Distributed systems, 1997, 8(5): 502-520.

[7] Wang W, Liu Y, Zhao Z, et al. Parallel Multiphysics Simulation of Package Systems Using an Efficient Domain Decomposition Method[J]. Electronics, 2021, 10(2): 158.

[8] 劉偉峰. 高可擴展, 高性能和高實用的稀疏矩陣計算研究進展與挑戰(zhàn)[J]. 數(shù)值計算與計算機應(yīng)用, 2020, 41(4): 259.

[9] GU H, REN H U, LIU C, et al. An optimized Chebyshev smoother In GAMG solver of openfoam on sunway Taihulight supercomputer[C]//The 13th OpenFOAM Workshop. 2018.