基于FMM-PM方法的宇宙N體模擬在GPU上的實現(xiàn)和優(yōu)化

扶月月，王武，王喬

1.中國科學(xué)院計算機(jī)網(wǎng)絡(luò)信息中心，北京 100190

2.中國科學(xué)院大學(xué)，北京 100049

3.中國科學(xué)院國家天文臺，北京 100101

引言

N體問題計算相互作用的N個粒子的相互作用和運(yùn)動，是高性能計算中具有代表性和挑戰(zhàn)性的問題之一，在天體物理、分子動力學(xué)、磁流體動力學(xué)等領(lǐng)域具有廣泛的應(yīng)用。模擬N體問題的數(shù)值方法包括直接法（Particle-Particle，PP）、樹方法（Tree Method）和粒子網(wǎng)格方法(Particle Mesh，PM)等，其中樹方法有BH 方法(Barnes-Hut)[1]和快速多極子方法 (Fast Multipole Method，F(xiàn)MM)[2]。PP方法的復(fù)雜度為O(N2)，其中N為粒子規(guī)模，這限制了模擬問題的規(guī)模；BH 方法的復(fù)雜度為O(NlogN)，且計算精度可控；FMM的復(fù)雜度為O(N)，兼具精度和速度優(yōu)勢，被IEEE 計算機(jī)協(xié)會和美國計算物理學(xué)會列為20世紀(jì)十大算法之一[3]；基于快傅里葉變換（Fast Fourier Transform，F(xiàn)FT）的PM 方法[4]復(fù)雜度為O(NlogN)，雖然計算速度較快，但計算短程力時不精確，只適用于長程力計算。為了彌補(bǔ)PM方法計算精度不足的缺點，PM 混合算法應(yīng)運(yùn)而生，PM 方法結(jié)合直接法和樹方法的混合算法（P3M[5]、TreePM[6]）應(yīng)用廣泛。

隨著高性能計算技術(shù)的發(fā)展，以及天文學(xué)和分子動力學(xué)等領(lǐng)域?qū)Υ笠?guī)模計算需求的增長，僅從算法上降低計算復(fù)雜度或者增加CPU 核數(shù)量并不足以滿足領(lǐng)域科學(xué)家的計算需求。采用硬件加速來加快計算速度，從而擴(kuò)大計算規(guī)模也是較為可行的途徑。GPU的最新發(fā)展已經(jīng)能夠以低廉的成本提供高性能的通用計算，基于CPU+GPU的并行異構(gòu)架構(gòu)也已成為當(dāng)今大型超級計算機(jī)的主流體系結(jié)構(gòu)之一。例如，2019年11月全球TOP500 排行榜前10名的超級計算機(jī)中，有5臺使用了CPU+GPU的異構(gòu)體系結(jié)構(gòu)。

近些年N體問題成為了GPU 上的熱點應(yīng)用之一，國際上有不少針對各個算法使用GPU 來加速N體模擬的研究，很多顯著的研究成果都是對PP 方法的N體模擬軟件進(jìn)行GPU 優(yōu)化加速[7-9]，基于BH方法和FMM 方法在大規(guī)模GPU 集群上的性能提升的研究成果也不少[10-12]，但是利用CPU+GPU 異構(gòu)平臺對PM 混合算法進(jìn)行優(yōu)化加速的研究并不多，目前還沒有在GPU 集群上實現(xiàn)FMM-PM 耦合算法的N體模擬。本文采用FMM-PM 耦合的方法實現(xiàn)N體問題的模擬，既可兼顧計算速度和精度，又能重疊短程力的計算與長程力的通信，因此較適合大規(guī)模并行。我們的工作為今后對多體問題和異構(gòu)體系結(jié)構(gòu)的研究工作奠定了堅實的基礎(chǔ)。

宇宙模擬軟件PHoToNs是針對宇宙暗物質(zhì)大尺度結(jié)構(gòu)設(shè)計的N體模擬軟件[13]，其第二個版本改進(jìn)了引力的計算算法，采用了FMM-PM 耦合算法實現(xiàn)，其中FMM的核心函數(shù)在GPU 上實現(xiàn)，F(xiàn)MM 中計算量較小的算子、PM函數(shù)和通信函數(shù)在CPU上實現(xiàn)，本文給出了該軟件在CPU+GPU 架構(gòu)上的實現(xiàn)過程和性能優(yōu)化技術(shù)。第1 節(jié)介紹FMM-PM 耦合算法的基本原理和計算流程，第2 節(jié)介紹FMM 在GPU 上的實現(xiàn)策略，第3節(jié)介紹FMM 在GPU 上的三種性能優(yōu)化方法，第4 節(jié)給出測試結(jié)果。

1 FMM-PM 算法原理

1.1 快速多極子方法

天文N體問題主要計算引力場中的N個粒子（星體、星系或暗物質(zhì)都可以抽象為粒子）的相互作用和運(yùn)動，粒子所受作用力可表示為：

其中i= 1,2,...,N，G是引力常量，mi和mj是粒子i和j的質(zhì)量，rij=xi-xj為粒子i到j(luò)的位移矢量。采用該公式直接計算每對粒子之間的相互作用，其精度較高且易于并行實現(xiàn)，然而需要每次重復(fù)計算各個時間步，這個過程的時間復(fù)雜度為O(N2)，模擬問題的粒子規(guī)模大大受限。

FMM 將粒子的相互作用分為遠(yuǎn)場和近場作用，遠(yuǎn)場作用通過近似方法來計算，在近場作用的計算過程中，通過層次劃分和位勢函數(shù)的多極子展開計算各點位勢，然后將位勢轉(zhuǎn)換為各點所受的力。位勢和引力滿足：

其中i= 1,2,...N,mi,ri,F(ri)分別為粒子的質(zhì)量、位置、位勢和作用力。雖然近場部分使用直接法，但由于短程力作用范圍有限，近場作用的粒子對數(shù)量遠(yuǎn)小于遠(yuǎn)場作用，所以FMM 整體復(fù)雜度仍為O(N)。

FMM的核心思想是：首先將空間進(jìn)行多層組劃分，然后通過核函數(shù)多極子展開，將粒子的位勢聚集到所在盒子的中心，逐層向上遞歸計算各層盒子的展開系數(shù)，最后把粒子間的遠(yuǎn)場相互作用轉(zhuǎn)換為盒子中心之間的相互作用，逐層向下遞歸到所有葉盒子。FMM 樹結(jié)構(gòu)遞歸過程的示意圖見圖1，向上箭頭表示M2M 逐層聚集，向下箭頭表示逐層發(fā)散，水平虛線表示次相鄰轉(zhuǎn)移。

圖1 FMM 樹結(jié)構(gòu)算法示意圖Fig.1 FMM tree structure algorithm diagram

FMM 逐層遞推過程包括以下七個基本步驟：

(1)空間多層盒子劃分，自下而上構(gòu)造樹；

(2)計算每個葉盒子的多極子展開系數(shù)和粒子與葉盒子的相互作用（P2M)，然后逐層上聚到父盒子中心（M2M）；

(3)計算次相鄰盒子的相互作用，即級數(shù)展開系數(shù)的轉(zhuǎn)換（M2L）；

(4)逐層向下計算各盒子的局部展開系數(shù)以及父盒子對子盒子的作用（L2L）；

(5)采用直接法計算每個葉盒子和相鄰盒子里所有粒子的近相互作用（Fnear-field）；

(6)計算所有粒子的位勢和受到的長程力相互作用（Ffar-field）；

(7)計算所有粒子的加速度，基于蛙跳格式更新下一時刻的速度和位置。

1.2 粒子網(wǎng)格方法

粒子網(wǎng)格（Particle Mesh，PM）方法基于譜空間的位勢場泊松方程求解粒子系統(tǒng)的引力，它應(yīng)用于具有周期邊界的N體模擬。PM 方法根據(jù)粒子信息得到均勻網(wǎng)格上的密度分布函數(shù)，求解出網(wǎng)格點的位勢，然后根據(jù)粒子所處的網(wǎng)格單元位置，通過差分和插值得到每個粒子所在點的位勢和所受的力。PM 方法的基本步驟如下：

(1)構(gòu)造密度分布：基于粒子的分布情況計算構(gòu)造均勻網(wǎng)格，以及格點的密度分布函數(shù)；

(3)計算作用力：通過有限差分和插值方法從格點的位勢得到粒子所受的力；

(4)計算加速度：根據(jù)粒子質(zhì)量，將每個粒子的受力轉(zhuǎn)換為加速度；

(5)速度和位置更新：基于時間步進(jìn)積分（可采用蛙跳格式、龍格-庫塔格式等）計算各粒子下一時刻的速度和位置。

與PP 方法相比，PM 方法的網(wǎng)格點數(shù)量與粒子數(shù)相當(dāng)或者更少，而且基于規(guī)則網(wǎng)格的密度場分布通過譜空間的FFT 求解各點的引力勢，僅涉及到兩次FFT，其復(fù)雜度為O(NglogNg)，其中Ng為離散網(wǎng)格點數(shù)，所以PM 計算速度較快，但是如果粒子間的距離較小時，該方法算出的粒子受力計算精確不能滿足要求，因此PM 只適合計算粒子的長程力相互作用。

1.3 FMM-PM 耦合方法

耦合方法把引力勢拆分成長程和短程部分，其中長程部分通過譜空間位勢的傅立葉變換算出，短程部分通過位形空間泊松方程求解得到，計算公式如下：

俗話說：“愛子愛在心里?！比欢?，他愛子是愛在臉上的。有一次，我應(yīng)邀在其老家作客，吃飯時，大人還沒有動筷子，他的大兒子卻大口大口地先吃起來，還說：“這是我要吃的，你們不能吃！”一點規(guī)矩都沒有。然鄭君只是輕描淡寫地說一句：“老大！客人在，你收斂些好嗎！”飯后，我和鄭君正在品茶談天，突然，有一孩子哭著來告狀，說是他額頭上的烏青塊是鄭君的大兒子用磚頭砸的。鄭君連句安慰他的話都沒有，卻反問那個孩子：“他為什么砸你，而不砸別人呢?”這無疑是在默認(rèn)、縱容其子的錯誤行為。

Hockney 等人在1974年提出P3M[5]方法，采用PP 方法計算短程力，它早期用于等離子體模擬，后來應(yīng)用于宇宙模擬；Xu G 等在1995年提出TreePM[6]方法，采用樹方法計算短程力，TreePM 方法帶來計算性能的顯著提升，在天體物理和分子動力學(xué)中得到廣泛應(yīng)用。目前國際主流的宇宙N體模擬軟件就是采用TreePM 方法（如GADGET[14]）。

天文N體模擬軟件PHoToNs-2.0 采用FMM-PM耦合算法，長程力采用基于FFT的PM 方法求得，短程力通過FMM 計算。基于FFT的PM 方法雖然速度更快，但受限于精度，只適用于計算長程力，而且會帶來全局通信。FMM 不僅具有近似線性的復(fù)雜度和較高的精度，而且適合并行。其中近場相互作用是計算密集型的，可有效利用高性能異構(gòu)平臺的加速器件（如GPU、MIC 等）。采用FMM-PM 耦合算法來實現(xiàn)N體問題的大規(guī)模并行，既兼顧了計算速度和精度，又重疊了短程力的計算與長程力的通信。

圖2為FMM-PM 耦合算法示意圖，虛線內(nèi)為FMM 計算區(qū)域，其中粒子間的短程力相互作用（P2P）占FMM 計算時間的90%以上，適合在加速卡上實現(xiàn)；虛線外的區(qū)域基于三維FFT的PM 方法計算長程力，占整體運(yùn)行時間的比例約1%，因此不需要采用GPU 加速，可在CPU 端調(diào)用FFT 函數(shù)庫(如FFTW、P3DFFT 等)實現(xiàn)，PHoToNs-2.0 采用二維區(qū)域分解的三維并行FFT 函數(shù)庫2decomp-FFT[15]。

圖2 FMM-PM 耦合算法示意圖Fig.2 FMM-PM algorithm diagram

2 FMM 在GPU 上的實現(xiàn)

FMM的本地和非本地短程作用（P2P）的計算時間占FMM-PM 運(yùn)行時間的比例最高（對于大規(guī)模問題，通常在90%以上），考慮到P2P 算子是計算密集型和數(shù)據(jù)密集型，可采用GPU 進(jìn)行加速。FMM其它計算（構(gòu)造樹，盒子-粒子作用等）是訪存、傳輸密集型函數(shù)，計算量相對較小。核心函數(shù)P2P_kernel是計算短程力的主要函數(shù)，其作用是通過粒子的位置、質(zhì)量等信息算出粒子的受力和加速度。P2P_kernel 在GPU 上實現(xiàn)的過程如下：

(1)CPU 獲取任務(wù)隊列，做數(shù)據(jù)傳輸準(zhǔn)備，開辟顯存空間；

(2)CPU 把需要計算的粒子信息傳輸?shù)斤@存；

(3)GPU 上計算粒子的短程力相互作用；

圖3 計算任務(wù)（包）與GPU 線程映射關(guān)系示意圖Fig.3 Mapping relationship diagram between computing tasks (packages) and GPU threads

(5) CPU 端更新粒子信息。

在具體實現(xiàn)，把按順序存放的粒子信息以及其鄰居粒子分批存入CPU 中開辟的緩沖池中，之后一起存入GPU 顯存中保證完整性，而且使GPU有序訪問顯存，提高存取效率。如圖3所示，根據(jù)CUDA 編程模型，一個線程映射一個粒子信息，充分發(fā)揮GPU 計算線程較多的優(yōu)勢。

3 FMM 在GPU 上的性能優(yōu)化

3.1 核心函數(shù)的參數(shù)優(yōu)化

在FMM 中MaxPackage 參數(shù)至關(guān)重要，它表示樹結(jié)構(gòu)最細(xì)層粒子包中的粒子數(shù)上限，其大小與任務(wù)量和數(shù)據(jù)傳輸開銷有關(guān)。隨著MaxPackage 增大，雖然單次執(zhí)行核函數(shù)的任務(wù)量增大，但總?cè)蝿?wù)次數(shù)減小，總數(shù)據(jù)傳輸開銷也隨之減小。反之，MaxPackage 減小，總數(shù)據(jù)傳輸開銷增大。短程力數(shù)據(jù)傳輸?shù)膹?fù)雜度可表示為：

其中Npart表示每個進(jìn)程的粒子數(shù)，常數(shù)K表示每個盒子的近場相鄰盒子數(shù)，包括該盒子本身，取值與截斷半徑和判斷準(zhǔn)則有關(guān)。

MaxPackage 參數(shù)分別與數(shù)據(jù)傳輸量和數(shù)據(jù)傳輸時間關(guān)系如圖4和圖5所示。

圖4 MaxPackage與數(shù)據(jù)傳輸量的關(guān)系Fig.4 The relationship between MaxPackage and the amount of data transmitte

圖5 MaxPackage與數(shù)據(jù)傳輸時間的關(guān)系Fig.5 The relationship between MaxPackage and the Data transfer time

圖6為使用nvprof 工具分析內(nèi)存拷貝時間與核心函數(shù)p2p_kenel的執(zhí)行時間分別占P2P 計算總執(zhí)行時間的比重。其中，memcpyHtoD和memcpyDtoH分別表示從主機(jī)端到設(shè)備端的內(nèi)存拷貝和從設(shè)備端到主機(jī)端的內(nèi)存拷貝。可見，在多GPU 系統(tǒng)中，數(shù)據(jù)傳輸時間占總P2P 短程作用力的計算總時間的90%以上。雖然短程力計算量與MaxPackage 成正比，但測試表明在GPU 上傳輸比計算更耗時間，即P2P計算復(fù)雜度可表示為：

因此，在GPU 上通過增大MaxPackage 參數(shù)，可以減少數(shù)據(jù)傳輸開銷，進(jìn)一步優(yōu)化P2P 短程作用力的計算和傳輸總時間。

3.2 頁鎖定內(nèi)存優(yōu)化和CUDA 流優(yōu)化

對CUDA 架構(gòu)而言，傳統(tǒng)的內(nèi)存拷貝方式是由操作系統(tǒng)API malloc( )在主機(jī)上分配的同時用cudaMalloc( )分配設(shè)備內(nèi)存，這種內(nèi)存，方式稱為可分頁內(nèi)存（page-able memory）。可分頁內(nèi)存的內(nèi)存拷貝非常耗時，為了加速內(nèi)存拷貝，CUDA 提供了頁鎖定內(nèi)存（page-lock memory 或 pinned memory）機(jī)制來而分配主機(jī)內(nèi)存，即采用API cudaHostAlloc( )分配主機(jī)內(nèi)存。使用頁鎖定內(nèi)存方式分配的主機(jī)內(nèi)存會常駐物理內(nèi)存中，操作系統(tǒng)將不會對該內(nèi)存進(jìn)行分頁和交換到磁盤上，這塊內(nèi)存不會被破壞或者重新定位，所以操作系統(tǒng)能夠安全地使某個應(yīng)用程序訪問該內(nèi)存的物理地址。CPU與GPU 內(nèi)存交互方式如圖7所示，CPU和GPU之間的總線是PCIe，是雙向傳輸?shù)模珻PU和GPU之間的數(shù)據(jù)拷貝使用DMA 機(jī)制來實現(xiàn)。

圖6 memcpyHtoD、memcpyHtoD 以及p2p_kernel所占的時間比重Fig.6 Time proportion of memcpyHtoD,memcpyHtoD and p2p_kernel

圖7 CPU與GPU之間的內(nèi)存交互示意圖Fig.7 Memory interaction diagram between CPU and GPU

在進(jìn)行內(nèi)存復(fù)制操作時，若采用傳統(tǒng)的可分頁方式，CUDA 驅(qū)動程序仍會通過DRAM 把數(shù)據(jù)傳給GPU，這樣復(fù)制操作會執(zhí)行兩遍：先從可分頁內(nèi)存復(fù)制一塊到臨時的頁鎖定內(nèi)存，再從這個臨時的頁鎖定內(nèi)存復(fù)制到GPU 上。當(dāng)從可分頁內(nèi)存中執(zhí)行復(fù)制時，復(fù)制速度將受限制于PCIE 總線的傳輸速度和系統(tǒng)前段速度相對較低的一方。使用了頁鎖定內(nèi)存時，GPU可以直接確定內(nèi)存的物理地址，從而能夠使用DMA 技術(shù)在GPU和CPU之間拷貝數(shù)據(jù)。顯然，使用頁鎖定內(nèi)存大大減少了CPU與GPU之間數(shù)據(jù)傳輸時間，最大化內(nèi)存吞吐量，提高內(nèi)存帶寬。

CUDA 流包括一系列異步的CUDA 操作（主機(jī)與設(shè)備間的數(shù)據(jù)傳輸，內(nèi)核啟動以及由主機(jī)發(fā)起但由設(shè)備處理的一些命令）。CUDA 流按照主機(jī)程序確定的順序在設(shè)備端執(zhí)行，相對于主機(jī)是異步的?？梢允褂肅UDA的API 來確保一個異步操作在運(yùn)行結(jié)果被使用前完成。在同一個CUDA 流中的操作有嚴(yán)格的執(zhí)行順序，而在不同CUDA 流中的操作在執(zhí)行順序上不受限制。

使用多個CUDA 流同時啟動多個kernel可實現(xiàn)網(wǎng)格級并發(fā)，因為所有排隊的操作都是異步的，所以在主機(jī)與設(shè)備系統(tǒng)中可重疊執(zhí)行部分操作。在同一時間內(nèi)將流中排隊的操作與其他并發(fā)操作一起執(zhí)行，可以隱藏部分操作的時間開銷，使設(shè)備利用率最大化。圖8為使用2個CUDA流的性能提升示意圖，可以看出，在一個流中設(shè)備端的計算kernel與另一個流中設(shè)備到主機(jī)的數(shù)據(jù)傳輸是重疊的。

3.3 混合精度和快速數(shù)學(xué)庫優(yōu)化

單精度和雙精度浮點運(yùn)算在通信和計算上的性能差異是不可忽略的。在我們的程序中，使用雙精度數(shù)值能夠使程序運(yùn)行總時間增加近一倍(這個結(jié)果取決于程序是計算密集型還是I/O 密集型)，在設(shè)備端進(jìn)行數(shù)據(jù)通信的時間也是使用單精度的兩倍，這是由于傳輸數(shù)組的字節(jié)長度相差兩倍。隨著全局內(nèi)存輸入/輸出數(shù)量和每條指令執(zhí)行的位操作數(shù)量的增加，設(shè)備上的計算時間也會增加。為了最大化指令吞吐率，我們使用單精度數(shù)(float)常量、變量和單精度計算函數(shù)，僅在有必要的時候使用雙精度運(yùn)算。

內(nèi)部函數(shù)__fdividef與除法運(yùn)算相比，在執(zhí)行浮點除法時速度更快但數(shù)值精度相對較低。用功能上等價的__fdividef 來替換除法操作，可以加速程序運(yùn)行。另外，快速數(shù)學(xué)庫（Fast Math Library）可以將核函數(shù)中使用的單精度計算函數(shù)替換為CUDA 內(nèi)部實現(xiàn)的高速版本。我們對核函數(shù)中單精度浮點類型的超越函數(shù)（erfc）調(diào)用快速數(shù)學(xué)庫，程序性能大大提升。

圖8 使用2個流的CUDA 操作性能提升示意圖Fig.8 Performance improvement diagram of CUDA operation using 2 streams

4 數(shù)值結(jié)果

在浪潮NF5280M5 服務(wù)器上進(jìn)行測試。系統(tǒng)環(huán)境配置：2個Intel Xeon Gold 6148（20 核心*2），主頻為2.4GHz，12個內(nèi)存條 DDR4 2400 DDR4 （32GB *12）；4 張NVIDIA TITAN V GPU卡，顯存容量12GBS，CUDA 版本為9.1，采用Intel 編譯器和“-O3”優(yōu)化。

對優(yōu)化后的CUDA 版本程序進(jìn)行加速效果測試，粒子規(guī)模為兩百萬（1283），使用4個進(jìn)程，在浪潮服務(wù)器上單個時間步的模擬。圖9為不同MaxPackage 三個實現(xiàn)版本程序的數(shù)據(jù)傳輸時間對比。其中，版本1是未進(jìn)行優(yōu)化的CUDA 程序，版本2是表示使用頁鎖定內(nèi)存和CUDA 流優(yōu)化后的程序，版本3為使用混合精度調(diào)用快速數(shù)學(xué)庫優(yōu)化的程序?？梢姡琈axPackage 越大，數(shù)據(jù)傳輸越快，系統(tǒng)性能越好。

圖9 不同MaxPackage 三個優(yōu)化版本的數(shù)據(jù)傳輸時間對比Fig.9 Data transfer time comparison between three optimized versions of different MaxPackage

圖10為PHoToNs-2.0的純MPI程序和三個CUDA 實現(xiàn)版本程序的短程力計算函數(shù)P2P 運(yùn)行時間的對比，圖11為三個CUDA 實現(xiàn)版本采用不同MaxPackage值時P2P的加速倍數(shù)。可以看出，當(dāng)MaxPackage為128 時，核心函數(shù)P2P 執(zhí)行時間從475.9s 降到了1.16秒，加速了410.3倍。

圖10 純MPI 程序與三個CUDA 實現(xiàn)版本程序的P2P 函數(shù)執(zhí)行時間對比Fig.10 Comparison of P2P function execution time between pure MPI program and three CUDA implementation version programs

圖11 不同MaxPackage 三個CUDA 實現(xiàn)版本程序的P2P 函數(shù)加速倍數(shù)Fig.11 P2P acceleration multiples of three CUDA versions of different MaxPackages

5 結(jié)論與展望

本文在CPU+GPU 混合異構(gòu)平臺的，對基于快速多極子方法和粒子網(wǎng)格方法的N體模擬軟件PHoToN-2.0 進(jìn)行CUDA 實現(xiàn)和性能優(yōu)化，包括優(yōu)化粒子包的參數(shù)以減少數(shù)據(jù)傳輸開銷、使用頁鎖定內(nèi)存提高內(nèi)存帶寬和最大化內(nèi)存吞吐量、使用CUDA流優(yōu)化最大化設(shè)備利用率、采用混合精度優(yōu)化最大化指令吞吐率、有效利用快速數(shù)學(xué)庫等手段。最終，F(xiàn)MM-PM 耗時90%以上的核心函數(shù)P2P 在在Titan V的GPU 平臺上采用4 張GPU 卡的計算速度相對采用4個Intel Xeon CPU 核提高了約410倍。這些優(yōu)化方法為在大規(guī)模GPU 異構(gòu)平臺上進(jìn)一步優(yōu)化和超大規(guī)模天文N體模擬奠定了基礎(chǔ)。

利益沖突聲明

所有作者聲明不存在利益沖突關(guān)系。