張 琨，賈金芳，嚴(yán)文昕，黃建強(qiáng)，2，王曉英

（1.青海大學(xué)計(jì)算機(jī)技術(shù)與應(yīng)用系，西寧 810016；2.清華大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)系，北京 100084）

0 概述

全球區(qū)域同化預(yù)報(bào)系統(tǒng)（Global and Regional Assimilation Prediction System，GRAPES）是我國(guó)自主研發(fā)的新一代數(shù)值天氣預(yù)報(bào)系統(tǒng)，能夠改善天氣預(yù)報(bào)模式的性能并有效提高天氣預(yù)報(bào)精度。GRAPES 動(dòng)力框架的計(jì)算核心是赫姆霍茲方程的求解，該方程是原始大氣方程組經(jīng)過(guò)一系列離散化處理之后形成的大規(guī)模稀疏線性方程組。目前對(duì)線性方程組的求解方法主要有直接法和迭代法［1-2］兩種。由于計(jì)算機(jī)硬件資源的限制，使用直接法并不能完成有效求解，而迭代法具有硬件存儲(chǔ)空間要求較小、原始矩陣在計(jì)算過(guò)程中保持不變等優(yōu)點(diǎn)。因此，目前求解大規(guī)模稀疏線性方程組主要選擇迭代法。

GRAPES 數(shù)值天氣預(yù)報(bào)系統(tǒng)采用的迭代法為廣義共軛余差（Generalized Conjugate Residual，GCR）法。文獻(xiàn)［3］利用基準(zhǔn)測(cè)試程序?qū)RAPES 并行應(yīng)用進(jìn)行分析，發(fā)現(xiàn)在迭代法計(jì)算過(guò)程中全局通信對(duì)整個(gè)系統(tǒng)的性能影響較大，指出在進(jìn)行性能優(yōu)化時(shí)需要考慮系統(tǒng)通信所帶來(lái)的性能下降問(wèn)題。文獻(xiàn)［4］針對(duì)赫姆霍茲方程求解過(guò)程中通信占比高的現(xiàn)象，引入短向量計(jì)算替代長(zhǎng)向量計(jì)算，以減少通信開(kāi)銷。文獻(xiàn)［5］在使用重啟動(dòng)和截?cái)鄡煞N方式進(jìn)行性能優(yōu)化的同時(shí)，還使用一種新的稀疏近似逆預(yù)條件子來(lái)加快算法收斂速度。

研究發(fā)現(xiàn)，增加重啟動(dòng)和截?cái)嗵幚淼膹V義共軛余差法能夠減少數(shù)據(jù)存儲(chǔ)量，同時(shí)改善計(jì)算的局部性。文獻(xiàn)［6］為改善GRAPES 模式動(dòng)力框架中存在的極點(diǎn)問(wèn)題，利用陰陽(yáng)網(wǎng)格設(shè)計(jì)新的動(dòng)力框架，與原GRAPES 模式框架相比，新的動(dòng)力框架表現(xiàn)出較好的數(shù)值穩(wěn)定性及計(jì)算性能。文獻(xiàn)［7］利用MPI 和OpenMP 混合并行的方式提升GRAPES 模式的并行度，以適應(yīng)主流硬件架構(gòu)既有分布式內(nèi)存又有共享內(nèi)存的情況。文獻(xiàn)［8］結(jié)合國(guó)產(chǎn)高性能平臺(tái)實(shí)現(xiàn)了基于MPI+共享變量編程模型的眾核線程級(jí)并行的多級(jí)并行方案，提升了方程求解的運(yùn)行效率。文獻(xiàn)［9］基于PETSc 科學(xué)計(jì)算庫(kù)和Hypre 預(yù)條件庫(kù)實(shí)現(xiàn)了廣義最小殘差（Generalized Minimal Residual，GMRES）法，與GRAPES 模式現(xiàn)有的GCR 算法相比，GMRES 算法在更高分辨率的條件下具有較好的性能。文獻(xiàn)［10］則利用CPU 設(shè)備結(jié)合線程級(jí)并行、矢量化和數(shù)據(jù)重用技術(shù)提升了算法計(jì)算性能。當(dāng)前GRAPES 模式動(dòng)力框架中廣義共軛余差法的優(yōu)化工作主要集中在通用CPU 及國(guó)產(chǎn)高性能機(jī)器領(lǐng)域，部分學(xué)者討論了算法本身的改進(jìn)方案，但大多是與CPU 設(shè)備進(jìn)行兼容，而結(jié)合GPU 的性能優(yōu)化相對(duì)較少［11］。

本文針對(duì)GRAPES 數(shù)值天氣預(yù)報(bào)系統(tǒng)中的赫姆霍茲方程求解問(wèn)題，分別通過(guò)MPI、MPI+OpenMP、CUDA 三種并行方式實(shí)現(xiàn)求解大規(guī)模稀疏線性方程組的廣義共軛余差法，并對(duì)測(cè)試結(jié)果進(jìn)行對(duì)比分析，評(píng)估優(yōu)化性能。

1 GRAPES模式中大規(guī)模稀疏線性方程組求解

GRAPES 數(shù)值天氣預(yù)報(bào)系統(tǒng)主要包括動(dòng)力框架和物理過(guò)程參數(shù)化方案兩部分［12-13］，其中動(dòng)力框架為核心部分，將動(dòng)力框架方程組進(jìn)行簡(jiǎn)化和離散后可得到GRAPES 模式的基本預(yù)報(bào)方程［14-16］。對(duì)該方程進(jìn)行處理后，動(dòng)力框架的主要計(jì)算內(nèi)容就變成對(duì)一個(gè)包含大型稀疏矩陣的赫姆霍茲方程的求解［17-19］。如圖1 所示，赫姆霍茲方程的求解過(guò)程計(jì)算量龐大，占據(jù)整個(gè)過(guò)程求解計(jì)算量的30%以上。隨著模式分辨率的提高，方程的計(jì)算量還會(huì)增加。因此，對(duì)赫姆霍茲方程進(jìn)行高效求解成為動(dòng)力框架性能提升的關(guān)鍵［20］。

圖1 GRAPES 模式核心計(jì)算Fig.1 GRAPES core computing

如圖2 所示，GRAPES 模式動(dòng)力框架部分的赫姆霍茲方程中計(jì)算每一個(gè)格點(diǎn)的方程需要選取空間中與其相關(guān)的19 個(gè)格點(diǎn)作為系數(shù)［21］。

圖2 赫姆霍茲方程系數(shù)空間分布Fig.2 Spatial distribution of the coefficients of Helmholtz equation

空間格點(diǎn)的計(jì)算公式如下：

將每一個(gè)空間格點(diǎn)的計(jì)算進(jìn)行組合后可以得到一個(gè)大規(guī)模的稀疏線性方程組，該方程組可簡(jiǎn)化為：

其中：A是一個(gè)大規(guī)模稀疏矩陣；x是規(guī)模為空間總格點(diǎn)數(shù)的解向量；b是與x相對(duì)應(yīng)的方程右端向量。因此，GRAPES 模式動(dòng)力框架中赫姆霍茲方程的求解就是對(duì)一個(gè)大規(guī)模稀疏線性方程組的求解。在當(dāng)前系統(tǒng)中，方程求解所采用的算法為廣義共軛余差法，在滿足精度要求的前提下，廣義共軛余差法算法收斂較快，整體算法易于實(shí)現(xiàn)，綜合表現(xiàn)較好［22］。廣義共軛余差法的具體內(nèi)容如算法1 所示。

算法1廣義共軛余差法

2 廣義共軛余差法的并行實(shí)現(xiàn)

2.1 問(wèn)題模型

本文主要考慮問(wèn)題模型Ax=b，等式中的A為稀疏矩陣，x為需要求解的向量，b為右端向量。系數(shù)矩陣具有明顯的稀疏性，在矩陣存儲(chǔ)格式的選擇上，統(tǒng)一使用CSR（Compressed Sparse Row）格式［23-24］。方程求解算法選擇廣義共軛余差法，算法主要計(jì)算內(nèi)容包括稀疏矩陣向量乘、向量?jī)?nèi)積、向量數(shù)乘。

2.2 系數(shù)矩陣預(yù)處理優(yōu)化

迭代法具有依賴系數(shù)矩陣條件數(shù)的特點(diǎn)，可能存在收斂速度慢甚至不收斂的情況，因此，通常使用預(yù)處理技術(shù)對(duì)算法進(jìn)行優(yōu)化。目前的預(yù)處理技術(shù)主要通過(guò)構(gòu)建預(yù)條件子將原方程Ax=b轉(zhuǎn)化為M(-1)Ax=M(-1)b等形式求解。常用的預(yù)條件子構(gòu)造方法有稀疏近似逆預(yù)條件、不完全分解預(yù)條件等。本實(shí)驗(yàn)所采用的預(yù)條件技術(shù)為不完全分解LU 預(yù)條件子（ILU）。預(yù)處理過(guò)程將原始矩陣A分為上三角矩陣U、下三角矩陣L及殘差矩陣R，分解之后的矩陣可以滿足不同條件P，例如與原系數(shù)矩陣的稀疏結(jié)構(gòu)保持一致。構(gòu)造ILU 預(yù)條件子的過(guò)程如算法2 所示。

算法2ILU 分解

2.3 MPI 并行及MPI+OpenMP 混合并行

通信開(kāi)銷是MPI 程序中不可忽視的部分，尤其是隨著進(jìn)程規(guī)模的增大，在計(jì)算時(shí)間減少的同時(shí)也會(huì)增加程序的通信開(kāi)銷［25-26］。因此，在設(shè)計(jì)MPI 并行算法時(shí)，需要平衡計(jì)算任務(wù)劃分與通信函數(shù)使用之間的關(guān)系。在廣義共軛余差法中，每一個(gè)迭代步驟之間都存在數(shù)據(jù)依賴，并行任務(wù)主要選擇在每一個(gè)迭代步驟內(nèi)部進(jìn)行。MPI 并行方式將計(jì)算任務(wù)根據(jù)分配的進(jìn)程數(shù)量進(jìn)行劃分，每個(gè)進(jìn)程處理不同的子任務(wù)，最后利用進(jìn)程間通信進(jìn)行數(shù)據(jù)同步。

本文實(shí)驗(yàn)所測(cè)試的MPI+OpenMP 混合并行方式，在MPI 并行方式的基礎(chǔ)上，利用多線程實(shí)現(xiàn)了更細(xì)粒度的任務(wù)劃分。OpenMP 是共享內(nèi)存并行編程模式，與MPI 并行模式相比可減少數(shù)據(jù)通信的時(shí)間。在矩陣向量乘部分的并行處理中，MPI+OpenMP 混合并行任務(wù)劃分如圖3 所示。矩陣數(shù)據(jù)按行子域劃分為不同的row_part，每個(gè)MPI 進(jìn)程只計(jì)算所屬row_part 與向量的乘積，各進(jìn)程并行執(zhí)行計(jì)算任務(wù)。同時(shí)在MPI 進(jìn)程內(nèi)部利用OpenMP 模式將row_part劃分為更細(xì)粒度的subrow_part，結(jié)合線程級(jí)并行完成子任務(wù)計(jì)算。

圖3 MPI+OpenMP 混合并行結(jié)構(gòu)Fig.3 MPI+OpenMP hybrid parallel structure

2.4 GPU 并行優(yōu)化

2.4.1 數(shù)據(jù)傳輸優(yōu)化

利用CUDA并行求解方程時(shí)主要開(kāi)銷在于GPU與CPU 之間的數(shù)據(jù)傳輸，為減少數(shù)據(jù)傳輸時(shí)間開(kāi)銷，避免PCI-E 總線上的數(shù)據(jù)傳輸成為程序性能提升的瓶頸，將矩陣A以及各個(gè)中間向量等都保存在GPU 上。在計(jì)算過(guò)程中利用顯存進(jìn)行通信，避免和主機(jī)端進(jìn)行頻繁的數(shù)據(jù)傳輸。無(wú)法避免的數(shù)據(jù)傳輸過(guò)程只出現(xiàn)在需要進(jìn)行收斂判斷及主機(jī)端計(jì)算的情況下，同時(shí)在數(shù)據(jù)傳輸時(shí)采用頁(yè)鎖定內(nèi)存（pinned memory）進(jìn)一步提高數(shù)據(jù)傳輸速度。如圖4 所示，與可分頁(yè)內(nèi)存（pageable memory）相比，頁(yè)鎖定內(nèi)存能夠減少一次數(shù)據(jù)傳輸操作，提高數(shù)據(jù)傳輸效率。

圖4 內(nèi)存模式對(duì)比Fig.4 Comparison of memory modes

2.4.2 訪存優(yōu)化

GPU 對(duì)全局存儲(chǔ)器的訪問(wèn)速度也會(huì)對(duì)程序性能造成一定影響，利用CUDA 模型可以實(shí)現(xiàn)線程塊間及線程塊內(nèi)兩級(jí)并行，從而完成對(duì)全局存儲(chǔ)器的合并訪問(wèn)。初始的數(shù)據(jù)存儲(chǔ)格式由于將需要連續(xù)訪問(wèn)的數(shù)據(jù)離散存儲(chǔ)，不滿足線程合并訪問(wèn)的條件，因此需要對(duì)數(shù)據(jù)存儲(chǔ)格式進(jìn)行轉(zhuǎn)換。將三維格式的矩陣轉(zhuǎn)化為一維格式存儲(chǔ)之后，每一個(gè)網(wǎng)格點(diǎn)在計(jì)算過(guò)程中同一線程束內(nèi)的線程所訪問(wèn)的內(nèi)存空間是連續(xù)的，滿足合并訪存的要求。如圖5 所示，當(dāng)同一線程束內(nèi)的線程訪存連續(xù)時(shí)，能夠?qū)⒍啻卧L存操作減少為一次訪存操作，從而減少計(jì)算過(guò)程中的訪存開(kāi)銷。

圖5 合并訪存Fig.5 Coalesced access

2.4.3 存儲(chǔ)器優(yōu)化

在GPU 計(jì)算過(guò)程中，存在一部分?jǐn)?shù)據(jù)只需要線程塊內(nèi)部的線程訪問(wèn)。對(duì)于這一部分?jǐn)?shù)據(jù)，利用共享存儲(chǔ)器進(jìn)行保存。共享存儲(chǔ)器的訪問(wèn)延遲比全局存儲(chǔ)器低，在計(jì)算過(guò)程中線程塊內(nèi)部的線程只需要訪問(wèn)共享存儲(chǔ)器中，而不需要去訪問(wèn)全局存儲(chǔ)器，從而進(jìn)一步提高訪存的速度。

在向量?jī)?nèi)積的計(jì)算過(guò)程中，并行方式如圖6 所示。內(nèi)積計(jì)算過(guò)程可以表示成形式，向量a和向量b分別對(duì)應(yīng)圖中Vector a 和Vector b。當(dāng)核函數(shù)啟動(dòng)后，將向量子任務(wù)subVector 分配到不同的線程塊中。由于向量?jī)?nèi)積涉及到對(duì)同一內(nèi)存地址的修改操作，因此先利用線程塊訪問(wèn)各自私有的共享存儲(chǔ)器進(jìn)行局部規(guī)約操作，這一部分計(jì)算內(nèi)容可以并行執(zhí)行。當(dāng)所有subVector 都完成局部歸約后，再進(jìn)行最后的規(guī)約操作，從而提高整體的計(jì)算速度。

圖6 CUDA 內(nèi)積并行結(jié)構(gòu)Fig.6 CUDA dot parallel structure

3 實(shí)驗(yàn)與分析

3.1 實(shí)驗(yàn)環(huán)境及測(cè)試數(shù)據(jù)

實(shí)驗(yàn)中共使用5 個(gè)CPU 節(jié)點(diǎn)，節(jié)點(diǎn)處理器為Intel?Xeon?CPU E5-2692 v2 @ 2.20 GHz。每個(gè)節(jié)點(diǎn)兩塊CPU（共24 核心），節(jié)點(diǎn)內(nèi)存為64 GB；實(shí)驗(yàn)使用GPU 為T(mén)esla T4，顯存容量為16 GB；CUDA 版本為v10.1。實(shí)驗(yàn)所使用的測(cè)試數(shù)據(jù)包含系數(shù)矩陣A、初始解x0、右端向量b，其中系數(shù)矩陣A規(guī)模為360×180×38，共2 462 400 個(gè)網(wǎng)格點(diǎn)。

3.2 結(jié)果分析

實(shí)驗(yàn)測(cè)試了數(shù)據(jù)在不同計(jì)算方式下的計(jì)算殘差。結(jié)果表明，使用3 種并行方式實(shí)現(xiàn)的廣義共軛余差法與串行方式能夠得到一致的計(jì)算結(jié)果，驗(yàn)證了不同并行方式求解稀疏線性方程組的正確性和有效性。

3.2.1 預(yù)處理優(yōu)化結(jié)果

對(duì)原始GCR 算法進(jìn)行預(yù)處理優(yōu)化后，方程的收斂速度獲得大幅提升。如圖7 所示，在誤差精度要求為1E-10 的條件下，預(yù)處理算法收斂所需的迭代次數(shù)從134 次下降至21 次。ILU 預(yù)處理的引入能夠改善方程的求解性能，同時(shí)也會(huì)增加新的計(jì)算內(nèi)容。但是在并行算法中，計(jì)算量已經(jīng)不是限制性能的最大因素，例如在MPI 程序中，隨著并行度的增加，進(jìn)程之間的通信成為限制程序性能的最大因素。在利用MPI 實(shí)現(xiàn)的并行GCR 算法中，每一次迭代過(guò)程都會(huì)進(jìn)行數(shù)據(jù)通信，所以，利用預(yù)處理技術(shù)加快算法收斂速度之后，整體的計(jì)算性能可以獲得提升。

圖7 預(yù)處理優(yōu)化對(duì)計(jì)算性能的影響Fig.7 Effect of preprocessing optimization on computational performance

3.2.2 MPI 并行結(jié)果

雖然使用MPI 并行方式能夠通過(guò)提高進(jìn)程數(shù)目減少子任務(wù)的計(jì)算量，從而有效提升計(jì)算部分速度，但是也會(huì)增加數(shù)據(jù)通信的時(shí)間開(kāi)銷。MPI 并行算法在實(shí)現(xiàn)過(guò)程中避免了冗余的通信開(kāi)銷，只對(duì)部分必要的向量更新操作進(jìn)行通信。MPI 并行算法在不同進(jìn)程規(guī)模下的運(yùn)行時(shí)間如圖8 所示?？梢钥闯?，隨著進(jìn)程規(guī)模的增加，運(yùn)行時(shí)間呈先下降后上升的趨勢(shì)。這是因?yàn)镸PI 算法的進(jìn)程數(shù)量會(huì)影響進(jìn)程之間的通信開(kāi)銷，當(dāng)進(jìn)程數(shù)為32 時(shí)，并行計(jì)算任務(wù)所帶來(lái)的性能提升不足以抵消增加的通信開(kāi)銷，程序整體性能就會(huì)下降。由于MPI 并行方式存在可擴(kuò)展性問(wèn)題，因此只通過(guò)增加進(jìn)程數(shù)量的方式來(lái)提高性能并不能高效地利用計(jì)算資源。

圖8 MPI 運(yùn)行時(shí)間Fig.8 Runtime of MPI

3.2.3 MPI+OpenMP 混合并行結(jié)果

與MPI 并行方式相比，MPI+OpenMP 混合并行方式充分結(jié)合了共享內(nèi)存的優(yōu)勢(shì)，在不改變?cè)糓PI 并行方式粗粒度通信的情況下，提高了進(jìn)程內(nèi)部的線程級(jí)并行度。圖9為MPI并行與MPI+OpenMP混合并行在不同核數(shù)規(guī)模下的計(jì)算性能對(duì)比?？梢钥闯?，當(dāng)計(jì)算所使用的核數(shù)過(guò)少時(shí)，利用OpenMP 進(jìn)行更細(xì)粒度的并行劃分之后性能并沒(méi)有明顯改善，但隨著核數(shù)的增加，OpenMP 的共享內(nèi)存優(yōu)勢(shì)得以發(fā)揮，計(jì)算性能相比MPI 并行方式有一定提升。當(dāng)核數(shù)為64 時(shí)，MPI 并行方式已經(jīng)出現(xiàn)性能下降，但是MPI+OpenMP 混合并行方式的性能依然較穩(wěn)定。因此，在MPI 并行方式由于可擴(kuò)展性原因造成計(jì)算效率下降的情況下，使用MPI+OpenMP 混合并行方式可以獲得更好的性能。

圖9 MPI 與MPI+OpenMP 運(yùn)行時(shí)間對(duì)比Fig.9 Runtime comparison between MPI and MPI+OpenMP

3.2.4 GPU 并行結(jié)果

基于GPU 的CUDA 并行方式在優(yōu)化措施上采用了共享存儲(chǔ)器優(yōu)化、全局存儲(chǔ)器訪存優(yōu)化以及主從設(shè)備數(shù)據(jù)傳輸優(yōu)化。CUDA 并行與MPI+OpenMP混合并行在算法求解過(guò)程中主要步驟的性能對(duì)比如圖10 所示。與MPI+OpenMP 混合并行相比，CUDA并行在矩陣向量乘部分（ar）有32%的性能提升，向量數(shù)乘部分（p&ap，x&r）有88%的性能提升，向量?jī)?nèi)積部分（beta，alpha，residual）有80%的性能提升。

圖10 CUDA 與MPI+OpenMP 主要步驟計(jì)算性能對(duì)比Fig.10 Performance comparison of main steps of CUDA and MPI+OpenMP

3.2.5 優(yōu)化性能對(duì)比

對(duì)不同優(yōu)化方式的優(yōu)化性能進(jìn)行對(duì)比，結(jié)果如圖11 所示?？梢钥闯觯ㄟ^(guò)預(yù)處理優(yōu)化減少算法迭代次數(shù)對(duì)整體計(jì)算性能提升明顯，加速比達(dá)到4.5。在預(yù)處理的基礎(chǔ)上，使用MPI 并行方式的性能獲得提升，加速比達(dá)到11.7。結(jié)合OpenMP 并行方式進(jìn)行細(xì)粒度任務(wù)劃分之后，混合并行計(jì)算性能相對(duì)于MPI 并行方式有一定提升，加速比達(dá)到15.8。同時(shí)數(shù)據(jù)也顯示CUDA 并行方式相對(duì)于MPI+OpenMP混合并行方式能夠獲得約50%的性能提升，加速比達(dá)到23.4?；贕PU 的并行方式在計(jì)算時(shí)間上的開(kāi)銷遠(yuǎn)小于基于CPU 的并行方式，這得益于GPU 數(shù)量眾多的計(jì)算核心。在數(shù)據(jù)通信方面，GPU 并行方式線程級(jí)的通信相比于CPU 并行方式也具有優(yōu)勢(shì)，不需要進(jìn)行各個(gè)進(jìn)程之間的數(shù)據(jù)通信，所以基于GPU的CUDA 并行方式能獲得較好的計(jì)算性能。

圖11 不同并行方式的優(yōu)化性能對(duì)比Fig.11 Optimization performance comparison of different parallel methods

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

大規(guī)模線性方程組的求解作為科學(xué)計(jì)算中的核心問(wèn)題，長(zhǎng)期以來(lái)都是研究者關(guān)注的重點(diǎn)，如何充分利用目前的高性能計(jì)算機(jī)資源對(duì)求解問(wèn)題的計(jì)算性能進(jìn)行優(yōu)化是重要的研究方向。本文從GRAPES 數(shù)值天氣預(yù)報(bào)系統(tǒng)動(dòng)力框架中的赫姆霍茲方程求解問(wèn)題出發(fā)，使用不同并行方式對(duì)廣義共軛余差法進(jìn)行并行及優(yōu)化，并對(duì)計(jì)算性能進(jìn)行對(duì)比分析。實(shí)驗(yàn)結(jié)果表明，在MPI 并行、MPI+OpenMP 混合并行及CUDA 并行3 種優(yōu)化方式中，基于GPU 的CUDA 并行方式能夠獲得更好的計(jì)算性能。后續(xù)將建立多機(jī)混合并行的優(yōu)化模型，測(cè)試MPI+CUDA 的多節(jié)點(diǎn)優(yōu)化性能，同時(shí)還將分析多重網(wǎng)格預(yù)處理、稀疏近似逆預(yù)處理等其他預(yù)處理方式對(duì)迭代法收斂速度的影響。