文敏華，劉永志，鮑 華，胡 躍，沈泳星，韋建文，林新華+

1.上海交通大學(xué) 高性能計(jì)算中心，上海 200240

2.上海交通大學(xué) 密西根學(xué)院，上海 200240

1 引言

DGX-2是NVIDIA家族最強(qiáng)大的人工智能系統(tǒng)，基于NVSwitch技術(shù)，整合16塊完全互聯(lián)的TeslaV100GPU構(gòu)建的可擴(kuò)展框架。整個(gè)系統(tǒng)顯存達(dá)到512 GB，共計(jì)擁有40 960個(gè)雙精度CUDA（compute unified device architecture）核心。除了多塊GPU卡帶來(lái)的強(qiáng)大計(jì)算能力，DGX-2中采用的NVSwitch技術(shù)提供了GPU之間的高速互聯(lián)，NVSwitch 擁有18 個(gè)51.5 GB/s 的NVLink 端口。基于其強(qiáng)大的計(jì)算能力以及GPU 間的高速通信網(wǎng)絡(luò)，DGX-2 在通用計(jì)算領(lǐng)域也擁有著巨大潛力。

聲子玻爾茲曼輸運(yùn)方程（Boltzmann transport equa-tion，BTE）能夠在亞連續(xù)尺度下模擬平衡或非平衡態(tài)的熱傳導(dǎo)現(xiàn)象，并且由于在大范圍長(zhǎng)尺度上的有效性持續(xù)受到歡迎[1]。在最近四十年中，半導(dǎo)體技術(shù)快速發(fā)展不僅使得芯片計(jì)算能力大大加強(qiáng)，也使得功耗越來(lái)越高；當(dāng)前熱問(wèn)題正成為決定芯片性能、可靠性和成本的主要因素。因此通過(guò)聲子BTE方程模擬計(jì)算半導(dǎo)體器件的熱傳導(dǎo)對(duì)了解其基本物理機(jī)制和發(fā)展散熱策略具有重要意義。但是當(dāng)前最為成熟的有限體積法求解聲子BTE方程來(lái)模擬工程實(shí)際問(wèn)題時(shí)，仍然存在計(jì)算量大、計(jì)算時(shí)間長(zhǎng)的困難。

基于將DGX-2強(qiáng)大的計(jì)算能力以及高速的通信速率應(yīng)用于聲子BTE 方程求解的初衷，本文做了如下工作：（1）使用CUDA框架對(duì)采用有限體積法的聲子BTE求解器的迭代部分進(jìn)行了針對(duì)GPU移植與并行優(yōu)化。在GPU 上實(shí)現(xiàn)了迭代計(jì)算的主要步驟，如聲子散射項(xiàng)的計(jì)算、線性方程組的求解、晶格溫度平衡態(tài)分布的計(jì)算等。（2）采用3 種方式，MPI+CUDA，NCCL（NVIDIA collective communications library）函數(shù)以及CUDA-AwareMPI，實(shí)現(xiàn)了對(duì)聲子BTE方程的多GPU的并行求解?？偠灾?，本文有以下兩點(diǎn)貢獻(xiàn)：

使用CUDA框架首次在GPU上實(shí)現(xiàn)了聲子BTE求解的迭代過(guò)程，并在單塊V100 上，較Intel Xeon Gold 6248單核性能提高了5.3倍至31.5倍。

在DGX-2 平臺(tái)上實(shí)現(xiàn)了對(duì)聲子BTE 的多GPU并行求解，測(cè)試了3 種并行方式對(duì)程序性能的影響，其中性能最優(yōu)的NCCL 庫(kù)函數(shù)版本在8 臺(tái)DGX-2 共128塊V100上實(shí)現(xiàn)了83%的強(qiáng)擴(kuò)展并行效率。

2 相關(guān)工作

有限體積法等確定性方法求解聲子玻爾茲曼輸運(yùn)方程較為困難。因?yàn)橹挥袑?duì)角度做到足夠數(shù)量的離散才能夠滿足網(wǎng)格無(wú)關(guān)性的驗(yàn)證，而這就使得該方法會(huì)產(chǎn)生數(shù)量龐大的離散方程，對(duì)計(jì)算量的需求較為龐大。針對(duì)這種情況，Ali 等[2]提出了求解聲子BTE的幾種大規(guī)模并行計(jì)算的策略和算法：（1）基于聲子模式；（2）基于角方向；（3）混合聲子模式和網(wǎng)格單元的并行方式。對(duì)三維器件類硅結(jié)構(gòu)進(jìn)行非平衡態(tài)的瞬態(tài)模擬，將計(jì)算域離散為604 054 個(gè)網(wǎng)格單元，使用400 個(gè)角度對(duì)角方向進(jìn)行離散，40 個(gè)聲子模式對(duì)聲子散射曲線進(jìn)行離散，使用400個(gè)進(jìn)程計(jì)算一個(gè)時(shí)間步長(zhǎng)都需要1 h。Ni等[3]也對(duì)有限體積法求解聲子BTE 方程進(jìn)行了并行計(jì)算研究，提出了基于空間域和聲子模式的并行策略。其中空間域較為適用全反射性的聲子BTE模型，能夠較好地?cái)U(kuò)展到128核并行。以上都是對(duì)聲子BTE 進(jìn)行CPU 并行求解，可以看出并行的核數(shù)為數(shù)百核左右，并行規(guī)模并不大。

關(guān)于使用GPU 進(jìn)行科學(xué)計(jì)算方面，在電動(dòng)力學(xué)領(lǐng)域，Priimak[4]在GPU 上實(shí)現(xiàn)了描述半導(dǎo)體超晶格中電子傳輸?shù)幕谟邢薏罘址ǖ牟柶澛斶\(yùn)方程，在GTX 680上相較于CPU串行版本獲得了118倍的加速。在流體力學(xué)領(lǐng)域，Calore等[5]在GPU上實(shí)現(xiàn)并優(yōu)化了可求解大規(guī)模湍流的格子玻爾茲曼方法，在K80 上，相較Xeon-Phi 獲得了4 倍加速比，并且最多擴(kuò)展到32個(gè)GPU并行計(jì)算，并行效率約90%，總的性能接近20 Tflops。在科學(xué)計(jì)算領(lǐng)域，Bell 等[6]基于不同稀疏矩陣存儲(chǔ)格式實(shí)現(xiàn)其相應(yīng)的矩陣向量乘，并在結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格上進(jìn)行測(cè)試，相較于4 核CPU 并行，GPU 計(jì)算達(dá)到了10 倍以上的加速比。Anzt 等[7]研究了預(yù)處理的BiCGSTAB（biconjugate gradient stabilized method）、CGS（conjugate gradient squared method）和QMR（quasi-minimal residual）等Krylov求解器在GPU上的效果，證明了在性能方面，相比于ILU預(yù)處理器，Jacobi預(yù)處理通常會(huì)有更高的效率。

從以上工作看出，GPU 在通用計(jì)算領(lǐng)域應(yīng)用較廣，但是目前還沒(méi)有在GPU 上對(duì)有限體積法求解聲子BTE方程進(jìn)行相關(guān)研究。本文不僅嘗試在GPU上對(duì)該方程進(jìn)行求解，且使用角方向的并行策略擴(kuò)展到128塊GPU上進(jìn)行計(jì)算，并取得83%的并行效率。

3 背景介紹

3.1 聲子玻爾茲曼輸運(yùn)方程

聲子玻爾茲曼輸運(yùn)方程是一個(gè)7 維非線性積分微分方程。目前求解該方程的方法基本有3種：（1）隨機(jī)性方法；（2）格子玻爾茲曼方法；（3）確定性方法。隨機(jī)性方法如蒙特卡羅法能有效求解這類高維偏微分方程[8]。但它最大的問(wèn)題在于求解實(shí)際工程問(wèn)題時(shí)需要避免結(jié)果的統(tǒng)計(jì)波動(dòng)，但這樣的代價(jià)過(guò)高[9]。格子玻爾茲曼方法僅用于在簡(jiǎn)單的二維結(jié)構(gòu)中求解聲子BTE方程[10-11]。同時(shí)，格子玻爾茲曼方法角度與空間離散不能相互獨(dú)立，離散方式不合適會(huì)極大影響結(jié)果的準(zhǔn)確性[12-13]?；陔x散的確定性方法，如離散坐標(biāo)法（discrete ordinate method，DOM）[14]、有限體積法（finite volume method，F(xiàn)VM）[15]和有限單元法（finite element method，F(xiàn)EM）[16]等廣泛用于聲子BTE的求解。采用離散方法產(chǎn)生的離散方程數(shù)目龐大，對(duì)求解的計(jì)算量需求也非常龐大。

聲子是晶格振動(dòng)的能量量子，其能量與頻率有關(guān)。頻率和波矢之間的關(guān)系為色散曲線，在計(jì)算中，連續(xù)色散曲線難以直接處理，因此采用兩種簡(jiǎn)化的計(jì)算模型：灰體模型和非灰模型。灰體模型將所有聲子歸于一個(gè)群體；而非灰模型，將色散曲線離散成若干個(gè)聲子模式。以下以非灰模型為例介紹聲子玻爾茲曼輸運(yùn)方程。通過(guò)結(jié)合聲子散射項(xiàng)和基于能量的聲子玻爾茲曼輸運(yùn)方程，可得到如下所示的在弛豫時(shí)間近似下的聲子能量玻爾茲曼輸運(yùn)方程：

3.2 基于有限體積法的聲子BTE求解

采用有限體積法對(duì)聲子BTE 方程進(jìn)行離散，將其變成一系列可求解的低維方程。對(duì)角度、計(jì)算域及色散曲線進(jìn)行離散，可得到如下形式的離散方程：

下標(biāo)i和向量S分別代表空間單元和角度方向，因此對(duì)于每一個(gè)聲子模式的每一個(gè)角方向及其每一個(gè)空間上的單元來(lái)說(shuō)，采用一階迎風(fēng)格式，以方向s1為例，式（2）將會(huì)有一個(gè)如下形式的離散方程：

在假定當(dāng)前的平衡態(tài)分布函數(shù)后，該離散方程會(huì)變成一個(gè)線性方程，而對(duì)每個(gè)角方向以及整體的區(qū)域來(lái)說(shuō)，將會(huì)有如下線性方程組：

系數(shù)矩陣中Kij表示第j個(gè)單元對(duì)第i個(gè)單元的影響系數(shù)。對(duì)該線性方程的求解可得到該方向下的能量密度分布。關(guān)聯(lián)角度之間及聲子模式的能量密度，可以求解溫度和平衡能量密度分布函數(shù)。

4 聲子BTE方程在DGX-2集群的移植優(yōu)化

在使用有限體積法求解聲子BTE 方程中，根據(jù)對(duì)色散曲線的簡(jiǎn)化不同分為灰體與非灰模型，本文采用非灰模型來(lái)求解聲子BTE方程。為求解非灰模型穩(wěn)態(tài)的能量密度分布和溫度分布，需要將角度間的能量密度進(jìn)行關(guān)聯(lián)，不僅需要考慮角度間的影響，還要考慮每個(gè)聲子模式間的影響。

4.1 聲子BTE方程求解過(guò)程在GPU上的實(shí)現(xiàn)

聲子BTE方程在GPU求解流程如圖1所示。算法的關(guān)鍵步驟如下：

（1）讀入網(wǎng)格文件、邊界條件以及參數(shù)文件，獲取所有聲子模式的遷移速度、比熱容和弛豫時(shí)間。

（2）假定每個(gè)聲子模式上所有方向的能量密度分布的值，初始化為0；計(jì)算離散后的角度控制方向。

（3）根據(jù)一階迎風(fēng)格式得到式（4）中每個(gè)聲子模式及其所有角方向上的系數(shù)矩陣。

Fig.1 Process of solving non-gray phonon BTE model圖1 聲子BTE非灰模型求解過(guò)程

（4）進(jìn)行CPU內(nèi)存與GPU顯存之間數(shù)據(jù)傳輸，主要數(shù)據(jù)為每個(gè)單元之間的影響系數(shù)。

（5）計(jì)算聲子散射項(xiàng)，使用BiCGSTAB 算法求解線性方程組，求解每個(gè)網(wǎng)格單元的能量密度分布。

（6）求解新的溫度分布和平衡態(tài)分布函數(shù)。

（7）計(jì)算能量密度分布新舊值之差，驗(yàn)證是否滿足收斂條件。若滿足收斂條件則停止計(jì)算，輸出結(jié)果，否則繼續(xù)（5）、（6）步驟。

由于整個(gè)迭代過(guò)程均在GPU 上實(shí)現(xiàn)，因此迭代過(guò)程所需數(shù)據(jù)一次性傳入顯存中，之后在GPU 上完成迭代過(guò)程中數(shù)據(jù)的更新，這樣減少了GPU 與CPU間的數(shù)據(jù)傳輸，提高了計(jì)算效率。通過(guò)上述步驟得到穩(wěn)態(tài)的能量密度的分布和溫度分布，之后可以求得熱流量，從而可以得到如熱導(dǎo)率這樣的宏觀特性，完成對(duì)介觀尺度下傳熱問(wèn)題的模擬。

4.1.1 線程分配策略

在GPU 上實(shí)現(xiàn)的迭代計(jì)算部分的函數(shù)，如聲子散射項(xiàng)函數(shù)、聲子模式溫度函數(shù)、晶格溫度函數(shù)以及平衡態(tài)分布函數(shù)，不同網(wǎng)格單元間沒(méi)有數(shù)據(jù)依賴。因此可使每個(gè)CUDA 線程對(duì)應(yīng)一個(gè)相應(yīng)的網(wǎng)格單元，GPU并行計(jì)算過(guò)程中，每個(gè)線程利用上一次迭代計(jì)算的結(jié)果或者是上一階段計(jì)算函數(shù)輸出的結(jié)果，完成相應(yīng)網(wǎng)格單元中的數(shù)據(jù)的計(jì)算，進(jìn)而完成整個(gè)迭代過(guò)程的計(jì)算。

4.1.2 數(shù)據(jù)存儲(chǔ)結(jié)構(gòu)

聲子BTE 方程維度較高，迭代計(jì)算過(guò)程所需的數(shù)據(jù)，如能量密度值，除空間維度，還額外擁有聲子模式以及角方向這兩個(gè)維度。其中空間維度使用網(wǎng)格編號(hào)表示，因此數(shù)據(jù)一般存儲(chǔ)在三維數(shù)組當(dāng)中。由于GPU 函數(shù)是以網(wǎng)格單元進(jìn)行的CUDA 線程分配，因此計(jì)算過(guò)程中相鄰線程間訪問(wèn)的是數(shù)組中同一聲子模式和角方向下，相鄰網(wǎng)格單元間的數(shù)據(jù)。為了能夠利用GPU合并訪存特性，簡(jiǎn)化CPU、GPU端數(shù)據(jù)傳輸，數(shù)據(jù)采用如圖2所示的一維數(shù)組進(jìn)行存儲(chǔ)。

Fig.2 Data storage format圖2 數(shù)據(jù)存儲(chǔ)結(jié)構(gòu)

4.2 BiCGSTAB算法實(shí)現(xiàn)流程

上一節(jié)描述的算法中最關(guān)鍵和耗時(shí)的步驟是通過(guò)求解離散后的聲子BTE線性方程來(lái)確定聲子能量密度分布（4.1節(jié)步驟5）。因?yàn)楸仨氠槍?duì)每個(gè)角方向和聲子模式求解聲子BTE 線性方程組，每一次迭代需進(jìn)行大量線性方程組求解。在本研究中，使用預(yù)處理穩(wěn)定雙共軛梯度法（BiCGSTAB）[17]求解線性方程。該算法如下所示：

Jacobi預(yù)處理方法易于在GPU并行計(jì)算，相較于ILU 預(yù)處理方法有更好的性能[7,17]，因此本文采用Jacobi 預(yù)處理方式，該方式為最簡(jiǎn)單的預(yù)處理形式，其預(yù)處理矩陣為原矩陣的對(duì)角矩陣。

整個(gè)求解過(guò)程流程如圖3所示。算法的主要部分是稀疏矩陣向量乘（sparse matrix-vector multiplication，SpMV）、點(diǎn)積（Dot）以及標(biāo)量向量乘和向量加（alpha X plus Y，AXPY）。其中CPU 負(fù)責(zé)求解過(guò)程的流程控制，核心計(jì)算部分在GPU 中進(jìn)行。上述算法在計(jì)算過(guò)程中，若出現(xiàn)關(guān)于r0的點(diǎn)積值為零，需要對(duì)x0、r0以及其他參數(shù)如ρ、α、ω等重新進(jìn)行初始化并繼續(xù)進(jìn)行迭代計(jì)算。因此在計(jì)算過(guò)程中，需要將一些參數(shù)傳輸?shù)絻?nèi)存，使用CPU 判斷是否需要重新初始化以及計(jì)算是否收斂。

4.2.1 GPU上稀疏矩陣向量乘的實(shí)現(xiàn)

SpMV 在線性方程組求解算法中占有重要的地位，許多學(xué)者對(duì)其進(jìn)行了研究。Bell 等[6,11]在GPU 實(shí)現(xiàn)了不同稀疏矩陣格式的SpMV計(jì)算，發(fā)現(xiàn)SpMV的性能主要受稀疏矩陣存儲(chǔ)格式以及核函數(shù)線程分配方式的影響。

由于本文研究的應(yīng)用中網(wǎng)格單元的系數(shù)，如3.2節(jié)式（4）所示，僅受相鄰3個(gè)網(wǎng)格單元的影響，即每行影響系數(shù)最多為3 個(gè)，因此整個(gè)系數(shù)矩陣為稀疏矩陣。由于ELL（Ellpack）格式適用于分布較為均衡的稀疏矩陣，且SpMV 實(shí)現(xiàn)方式簡(jiǎn)單，因此本文采用ELL格式稀疏矩陣。

Fig.3 Implementation of BiCGSTAB on GPU圖3 BiCGSTAB算法在GPU上的實(shí)現(xiàn)流程

Fig.4 ELL sparse matrix storage format圖4 ELL稀疏矩陣存儲(chǔ)格式

ELL格式存儲(chǔ)方式如圖4所示，采用兩個(gè)二維數(shù)組來(lái)存儲(chǔ)一個(gè)n×k的矩陣（k為包含非零元素最多行的非零元素?cái)?shù)目）。在實(shí)際使用中，使用兩個(gè)一維向量按列方向?qū)仃囘M(jìn)行存儲(chǔ)。使用ELL 格式的SpMV 算法在GPU 內(nèi)易于并行實(shí)現(xiàn)，計(jì)算過(guò)程中每個(gè)CUDA 線程計(jì)算矩陣的一行，由ELL 存儲(chǔ)方式可知，CUDA線程對(duì)矩陣值以及列索引的訪問(wèn)均是連續(xù)的，能夠充分利用GPU合并訪存特性。

4.2.2 GPU上點(diǎn)積計(jì)算的實(shí)現(xiàn)

點(diǎn)積計(jì)算需要進(jìn)行GPU規(guī)約計(jì)算。為了規(guī)約的方便，劃分塊時(shí)，塊內(nèi)線程數(shù)為2的指數(shù)倍。該方式可使每個(gè)線程訪問(wèn)共享內(nèi)存時(shí)不會(huì)發(fā)生bankconflict?？紤]到GPU 中warp 內(nèi)線程是同步的，不必進(jìn)行顯式線程同步。因此當(dāng)活動(dòng)線程均在同一warp 時(shí)，可直接進(jìn)行warp 內(nèi)規(guī)約計(jì)算，進(jìn)而減少了線程同步的次數(shù)，提高了計(jì)算效率。

由于采用多個(gè)塊（Block）進(jìn)行點(diǎn)積計(jì)算，因此整個(gè)規(guī)約的過(guò)程分為兩個(gè)階段：第一階段采用多個(gè)Block，首先計(jì)算出向量中每個(gè)值的乘積，之后進(jìn)行塊內(nèi)規(guī)約計(jì)算，并將當(dāng)前塊內(nèi)規(guī)約結(jié)果順序存入中間數(shù)組中；第二階段僅使用1個(gè)Block，對(duì)中間數(shù)組進(jìn)行規(guī)約求和，規(guī)約方式同上，最終可求得向量點(diǎn)積的值。

4.2.3 GPU上AXPY實(shí)現(xiàn)

在GPU 上并行計(jì)算AXPY 較為簡(jiǎn)單，為向量中每個(gè)數(shù)據(jù)分配一個(gè)線程即可。同時(shí)可將前后步驟間有數(shù)據(jù)依賴且并行方式相同的函數(shù)合并到一個(gè)CUDA 內(nèi)核中，這樣就可避免雙倍的數(shù)據(jù)加載，減少了內(nèi)核數(shù)量，節(jié)約內(nèi)核調(diào)用開(kāi)銷。例如在求解線性方程組中，會(huì)有如下所示的相鄰計(jì)算。

其中，K為Jacobi 預(yù)處理矩陣，實(shí)際計(jì)算過(guò)程為向量乘，并行過(guò)程與AXPY 一致，因此可以將兩步操作融合成一個(gè)CUDA 函數(shù)。該方式不僅簡(jiǎn)化了代碼，而且更為重要的是利用了數(shù)據(jù)的局部性，每一個(gè)線程計(jì)算的結(jié)果可直接應(yīng)用到下一階段計(jì)算中。

4.3 基于角方向的大規(guī)模并行策略

聲子玻爾茲曼輸運(yùn)方程通常有3 種并行策略：（1）計(jì)算域分解，該方式能做到較大規(guī)模擴(kuò)展，但是本文應(yīng)用針對(duì)的是二維計(jì)算域，網(wǎng)格規(guī)模較小，若將其分解成更小的網(wǎng)格，則難以發(fā)揮GPU 多線程的優(yōu)勢(shì)；（2）基于聲子模式的并行，該方式雖然較為容易實(shí)現(xiàn)，但是聲子模式的數(shù)量一般為數(shù)十個(gè)，不能做到較大規(guī)模的并行；（3）基于角方向的并行，角方向的數(shù)量比聲子模式大一個(gè)數(shù)量級(jí)，采用角方向并行能夠進(jìn)行數(shù)百個(gè)進(jìn)程并行計(jì)算。因此為了做到較大規(guī)模的并行計(jì)算，本研究采用基于角方向的并行策略。本文使用3種不同方法實(shí)現(xiàn)了該并行策略，其中第一種實(shí)現(xiàn)方式如圖5所示。

Fig.5 Parallel strategy based on angle directions圖5 基于角方向的并行策略

4.3.1 MPI+CUDA

如圖5 所示使用的是MPI+CUDA 的并行策略，即一個(gè)MPI 進(jìn)程控制一個(gè)GPU 設(shè)備。在初始化階段，獲取MPI進(jìn)程號(hào)rank，以及節(jié)點(diǎn)上GPU設(shè)備數(shù)量gpu_count，并將MPI進(jìn)程映射到GPU設(shè)備上gpu_id=rank%gpu_count。在DGX-2 節(jié)點(diǎn)上，一個(gè)節(jié)點(diǎn)有16個(gè)V100，因此一個(gè)節(jié)點(diǎn)分配16 個(gè)MPI 進(jìn)程，每個(gè)進(jìn)程對(duì)應(yīng)一個(gè)GPU設(shè)備。

求解聲子BTE 方程的迭代過(guò)程中，不同角方向的聲子散射所需數(shù)據(jù)為上一次迭代計(jì)算結(jié)果，與本次迭代計(jì)算無(wú)關(guān)。因此不同角方向可以并行執(zhí)行。不同MPI 進(jìn)程計(jì)算不同角方向，并在GPU 上完成該角方向線性方程組求解。由于在計(jì)算聲子模式的溫度分布中，每個(gè)網(wǎng)格單元需要對(duì)其所有角方向的能量密度進(jìn)行求和。因此每個(gè)進(jìn)程求得該角方向的能量密度后，需將結(jié)果數(shù)據(jù)從設(shè)備端傳輸?shù)街鳈C(jī)端，然后通過(guò)MPI_Allgather函數(shù)使得每個(gè)進(jìn)程均可獲得所有角方向的能量密度；最后將所有角方向能量密度傳輸?shù)紾PU中完成對(duì)聲子模式溫度的計(jì)算。然后繼續(xù)計(jì)算其他聲子模式下角方向的相關(guān)數(shù)據(jù)，最終求得當(dāng)前迭代過(guò)程中的晶格溫度和平衡態(tài)分布。

4.3.2 CUDA-AwareMPI

如4.3.1 小節(jié)所述的常規(guī)的MPI 實(shí)現(xiàn)，需首先將數(shù)據(jù)傳輸至CPU內(nèi)存中，之后使用MPI函數(shù)；而使用CUA-AwareMPI 實(shí)現(xiàn)，可直接傳輸GPU 中的數(shù)據(jù)。從圖6、圖7 兩種實(shí)現(xiàn)方式可以看出，使用CUDA-AwareMPI的實(shí)現(xiàn)，其代碼更簡(jiǎn)潔，編程更容易。

Fig.6 Traditional MPI send and receive scheme圖6 傳統(tǒng)MPI send receive模式

Fig.7 CUDA-AwareMPI implementation圖7 CUDA-AwareMPI實(shí)現(xiàn)方式

除此之外，使用CUDA-AwareMPI 能夠使數(shù)據(jù)傳輸?shù)牟僮鞫紩?huì)被流水化且能夠透明化地使用GPUDirect-RDMA的加速技術(shù)。該技術(shù)使得GPU中的數(shù)據(jù)可以直接被送到網(wǎng)卡進(jìn)行傳輸，從而消除了GPU到CPU設(shè)備的時(shí)間消耗。這樣也就顯著增大了GPU和其他節(jié)點(diǎn)的通信效率。

使用CUDA-AwareMPI 與4.3.1 小節(jié)所述的并行方式類似，僅需進(jìn)行兩點(diǎn)改動(dòng)即可，一是省略CPU和GPU 間數(shù)據(jù)通信的過(guò)程，二是將MPI 函數(shù)傳輸?shù)臄?shù)據(jù)改為GPU緩沖區(qū)數(shù)據(jù)。

4.3.3 NCCL函數(shù)

NCCL[18]是NVIDIA提供的GPU間的集合通信函數(shù)庫(kù)，實(shí)現(xiàn)了all-reduce、all-gather 和reduce-scatter 等集合通信函數(shù)。可以在使用PCIe、NVLink 和NVSwitch的平臺(tái)上實(shí)現(xiàn)高通信速率。其中2.0 版本可以進(jìn)行跨節(jié)點(diǎn)間的GPU 通信，并可進(jìn)行網(wǎng)絡(luò)拓?fù)錂z測(cè)以及自動(dòng)使用GPU DirectRDMA。為了在原并行程序的基礎(chǔ)上使用NCCL函數(shù)，采用一個(gè)進(jìn)程分配一個(gè)GPU的方式。

NCCL庫(kù)的使用方式為，首先完成對(duì)MPI的初始化，之后在0 號(hào)進(jìn)程上創(chuàng)建unique ID 并廣播到其他進(jìn)程，最后創(chuàng)建NCCL 的通信子。上述步驟完成后，即可使用NCCL中的集合通信函數(shù)。

使用NCCL 函數(shù)對(duì)聲子BTE 進(jìn)行并行求解，僅需在4.3.1小節(jié)所述方式上進(jìn)行如下3點(diǎn)修改：（1）增加NCCL相關(guān)初始化函數(shù)；（2）略去CPU與GPU數(shù)據(jù)傳輸過(guò)程；（3）使用ncclAllGather函數(shù)進(jìn)行集合通信。

5 實(shí)驗(yàn)結(jié)果及分析

5.1 實(shí)驗(yàn)環(huán)境

本研究實(shí)驗(yàn)所采用硬件配置如下，CPU 為Intel Xeon Gold 6248,架構(gòu)為Cascade Lake，主頻為2.5 GHz，核心數(shù)為20，雙精度浮點(diǎn)性能為1.6 Tflops。節(jié)點(diǎn)內(nèi)存大小為192 GB；GPU版本測(cè)試使用的是NVIDIA的DGX-2平臺(tái)，其中該平臺(tái)采用的GPU為V100，5 120個(gè)單精度CUDA核心，2 560個(gè)雙精度CUDA核心，雙精度浮點(diǎn)性能達(dá)到7.5 Tflops，顯存帶寬可達(dá)900 GB/s；DGX-2 內(nèi)部通過(guò)NVSwtich 橋接器支持16 個(gè)全互聯(lián)的V100 GPU卡。

性能測(cè)試對(duì)比采用CPU-GPU 并行計(jì)算和僅采用CPU 串行計(jì)算的加速效果。其中GPU 計(jì)算采用的CUDA 框架版本為9.2，串行CPU版本中進(jìn)行線性方程求解采用的是PETSc（portable，extensible toolkit for scientific computation）[19]中函數(shù)，其中PETSc版本為3.10。對(duì)4個(gè)算例進(jìn)行了GPU加速性能的測(cè)試，其中每個(gè)測(cè)試算例均為二維計(jì)算域，且計(jì)算域大小、邊界條件、遷移速度和弛豫時(shí)間等參數(shù)均相同，網(wǎng)格及網(wǎng)格數(shù)量不同。算例網(wǎng)格數(shù)量分別為6 282、11 974、18 060及24 912。其中網(wǎng)格數(shù)量為24 912的算例，在實(shí)際二維問(wèn)題的研究中屬于規(guī)模較大的算例。串行版本中不同網(wǎng)格算例內(nèi)存使用量分別為0.90 GB、1.85 GB、3.68GB、6.70GB，GPU顯存使用量約為0.30GB、0.81GB、1.00 GB、1.24 GB。內(nèi)存數(shù)據(jù)量較大主要是因?yàn)镃PU計(jì)算過(guò)程中需要使用二維矩陣，并將其壓縮為ELL格式稀疏矩陣，而GPU顯存中僅存儲(chǔ)壓縮后的ELL稀疏矩陣。聲子BTE方程求解過(guò)程中，離散聲子色散曲線的聲子模式個(gè)數(shù)為1，離散的角方向的個(gè)數(shù)為256。

5.2 GPU加速性能測(cè)試

圖8 所示為使用單卡V100 相較于Intel Xeon Gold 6248串行版本的迭代過(guò)程的加速效果?？v坐標(biāo)為單次迭代計(jì)算的時(shí)間，橫坐標(biāo)為不同網(wǎng)格大小的算例。由于本文旨在在GPU上實(shí)現(xiàn)BTE迭代求解過(guò)程，因此僅比較了GPU 版本與CPU 串行版本間的加速性能，未與CPU并行進(jìn)行對(duì)比。

Fig.8 Speedup of GPU for different mesh counts圖8 不同網(wǎng)格數(shù)目的GPU加速比

圖8中顯示的加速比針對(duì)的是迭代部分的時(shí)間，這是因?yàn)檎麄€(gè)迭代的步數(shù)在103至104量級(jí)之間，迭代過(guò)程之外部分的耗時(shí)相對(duì)較少，對(duì)加速比的影響較小，因此僅比較迭代過(guò)程的加速比。迭代過(guò)程的CPU-GPU 混合計(jì)算需考慮數(shù)據(jù)傳輸?shù)拈_(kāi)銷，本文的實(shí)現(xiàn)方式是在迭代計(jì)算前將GPU迭代所需數(shù)據(jù)全部傳入顯存中，迭代過(guò)程中的數(shù)據(jù)在GPU 顯存中獲取并計(jì)算更新。因此迭代過(guò)程中，GPU 僅需將與計(jì)算流程有關(guān)的系數(shù)傳入內(nèi)存。迭代過(guò)程中主機(jī)端與設(shè)備端數(shù)據(jù)傳輸為迭代時(shí)間的1.1%～1.2%。這是由于迭代過(guò)程中傳輸?shù)臄?shù)據(jù)多為單個(gè)參數(shù)，即單個(gè)浮點(diǎn)數(shù)，數(shù)據(jù)量小。

從圖8 中可以看出，隨著測(cè)試算例網(wǎng)格數(shù)量增加，GPU加速效果越來(lái)越明顯。算例1中網(wǎng)格數(shù)量為6 282，算例4為24 912，算例4網(wǎng)格數(shù)量約為算例1的4 倍，其串行CPU 版本迭代時(shí)間為算例1 的15.6 倍，而其GPU 版本迭代時(shí)間僅為算例1 的2.6 倍。這是由于GPU 并行過(guò)程中，多數(shù)函數(shù)中CUDA 線程數(shù)與網(wǎng)格數(shù)有關(guān)，而當(dāng)網(wǎng)格數(shù)量較小時(shí)，其使用的CUDA線程（thread）數(shù)也較少，未能充分利用GPU 的性能。同時(shí)GPU可以通過(guò)線程塊的切換掩蓋數(shù)據(jù)的訪存延遲，因此當(dāng)網(wǎng)格數(shù)量較大，線程塊數(shù)量較多時(shí)，可以充分利用該特性，提高并行效率，因而在本例中GPU對(duì)越大規(guī)模數(shù)據(jù)進(jìn)行并行加速，其效果越明顯。

5.3 多GPU并行性能測(cè)試

采用4.3節(jié)基于角方向的并行策略，對(duì)多GPU并行強(qiáng)擴(kuò)展性能進(jìn)行測(cè)試。由于本文旨在實(shí)現(xiàn)聲子BTE 程序的GPU 加速以及實(shí)現(xiàn)多GPU 版本的聲子BTE 程序，因此多GPU 并行的加速比是以單GPU 迭代時(shí)間為基準(zhǔn)，并未與CPU 并行進(jìn)行對(duì)比。測(cè)試使用算例為上述網(wǎng)格大小為24 912 的算例。每個(gè)DGX-2 節(jié)點(diǎn)有16 塊V100GPU，共有8 個(gè)節(jié)點(diǎn)。最多擴(kuò)展到128個(gè)進(jìn)程，其強(qiáng)擴(kuò)展性能如圖9所示。

Fig.9 Strong scalability for different parallel methods圖9 不同并行方式的強(qiáng)可擴(kuò)展性

圖9 中橫坐標(biāo)為GPU 個(gè)數(shù)，一個(gè)MPI 進(jìn)程對(duì)應(yīng)一個(gè)GPU設(shè)備，圖中標(biāo)注出了并行效率最好的NCCL庫(kù)函數(shù)方式，以及并行效率最差MPI+CUDA 方式的加速比。其中使用MPI+CUDA 的版本，擴(kuò)展到128卡上，相較于1 塊GPU 卡，其加速比為68，其并行效率為53%；使用CUDA-AwareMPI 技術(shù)的并行版本，128卡上加速比為97，并行效率為76%；使用NVIDIA的NCCL庫(kù)實(shí)現(xiàn)并行的版本，128卡上加速比為107，其并行效率為83%。采用NCCL 庫(kù)以及CUDA-AwareMPI的方式，在同一節(jié)點(diǎn)可利用NVSwitch實(shí)現(xiàn)GPU間直接高效通信，節(jié)點(diǎn)間通信則利用了GPUDirect-RDMA技術(shù)。而MPI+CUDA方式，首先需要通過(guò)PCIe（peripheral component interconnect express）將GPU 中的數(shù)據(jù)傳輸至CPU，之后進(jìn)行節(jié)點(diǎn)內(nèi)通信以及將數(shù)據(jù)傳入IB（infiniBand）網(wǎng)卡進(jìn)行節(jié)點(diǎn)間通信。由于通信方式的不同，使得使用NCCL函數(shù)相比MPI+CUDA實(shí)現(xiàn)方式性能提升57%。由此可以看出，基于NVLink和NVSwitch 技術(shù)的GPU 間通信相比傳統(tǒng)方式更快速和高效。

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

本文提出了采用非結(jié)構(gòu)化網(wǎng)格的有限體積法求解聲子玻爾茲曼輸運(yùn)方程（BTE）的GPU并行加速方法，并在8臺(tái)DGX-2集群上進(jìn)行了移植和并行優(yōu)化。在GPU上實(shí)現(xiàn)了整個(gè)迭代過(guò)程，減少了主機(jī)端與設(shè)備端的數(shù)據(jù)傳輸，同時(shí)在求解線性方程組時(shí)，通過(guò)規(guī)約過(guò)程的循環(huán)展開(kāi)以及內(nèi)核融合等優(yōu)化方法，在單塊V100GPU上，網(wǎng)格單元為2.4萬(wàn)的算例，相較于Intel Xeon Gold 6248上的串行版本獲得了31.5倍的性能提升。

同時(shí)，本文實(shí)現(xiàn)和評(píng)估了3 種多節(jié)點(diǎn)多GPU 并行的方法。采用MPI+CUDA的并行方式并行效率最低，而采用CUDA-AwareMPI以及使用NCCL函數(shù)均能充分利用NVSwitch和GPUDirectRDMA帶來(lái)的通信性能增益。其中使用NCCL 庫(kù)函數(shù)的并行方法最為高效，在8 臺(tái)DGX-2（128 塊NVIDIAV100）上獲得83%的并行效率，比純MPI版本提升達(dá)57%。

目前僅實(shí)現(xiàn)了基于角方向的并行策略，該方法需對(duì)網(wǎng)格能量密度進(jìn)行集合通信，由于當(dāng)前網(wǎng)格規(guī)模并不大，因此對(duì)性能影響較小。若對(duì)三維問(wèn)題并行求解，網(wǎng)格規(guī)模會(huì)提升一個(gè)數(shù)量級(jí)，則其對(duì)性能的影響會(huì)較為明顯。因此未來(lái)會(huì)考慮實(shí)現(xiàn)基于聲子模式的并行，該方式通信頻率較少，但是會(huì)有負(fù)載不均衡問(wèn)題，因此需采取相應(yīng)策略進(jìn)行優(yōu)化。此外本研究中僅采用Jacobi 預(yù)處理方式，雖然其較容易并行，但是收斂較差，因此接下來(lái)考慮測(cè)試不同的預(yù)處理矩陣的性能。同時(shí)將會(huì)測(cè)試更多的算例并且嘗試使用不同的矩陣格式與線性方程求解器來(lái)優(yōu)化單塊GPU上的性能。

致謝感謝上海交通大學(xué)高性能計(jì)算中心程盛淦老師在DGX-2設(shè)備的使用上提供的幫助。