扶月月，王武，王喬

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

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

3.中國(guó)科學(xué)院國(guó)家天文臺(tái)，北京 100101

引言

N體問題計(jì)算相互作用的N個(gè)粒子的相互作用和運(yùn)動(dòng)，是高性能計(jì)算中具有代表性和挑戰(zhàn)性的問題之一，在天體物理、分子動(dòng)力學(xué)、磁流體動(dòng)力學(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)，且計(jì)算精度可控；FMM的復(fù)雜度為O(N)，兼具精度和速度優(yōu)勢(shì)，被IEEE 計(jì)算機(jī)協(xié)會(huì)和美國(guó)計(jì)算物理學(xué)會(huì)列為20世紀(jì)十大算法之一[3]；基于快傅里葉變換（Fast Fourier Transform，F(xiàn)FT）的PM 方法[4]復(fù)雜度為O(NlogN)，雖然計(jì)算速度較快，但計(jì)算短程力時(shí)不精確，只適用于長(zhǎng)程力計(jì)算。為了彌補(bǔ)PM方法計(jì)算精度不足的缺點(diǎn)，PM 混合算法應(yīng)運(yùn)而生，PM 方法結(jié)合直接法和樹方法的混合算法（P3M[5]、TreePM[6]）應(yīng)用廣泛。

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

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

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

1 FMM-PM 算法原理

1.1 快速多極子方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 FMM-PM 耦合方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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表示每個(gè)進(jìn)程的粒子數(shù)，常數(shù)K表示每個(gè)盒子的近場(chǎng)相鄰盒子數(shù)，包括該盒子本身，取值與截?cái)喟霃胶团袛鄿?zhǔn)則有關(guān)。

MaxPackage 參數(shù)分別與數(shù)據(jù)傳輸量和數(shù)據(jù)傳輸時(shí)間關(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ù)傳輸時(shí)間的關(guān)系Fig.5 The relationship between MaxPackage and the Data transfer time

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論與展望

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

利益沖突聲明

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