趙文龍，王武

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

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

引言

N 體問(wèn)題描述了多個(gè)粒子在經(jīng)典力學(xué)下的運(yùn)動(dòng)規(guī)律，是廣泛運(yùn)用于天體物理學(xué)和分子動(dòng)力學(xué)等多個(gè)重要領(lǐng)域中的基本問(wèn)題之一。為了能夠更加精準(zhǔn)地預(yù)測(cè)宇宙中物質(zhì)的分布以及預(yù)想模型的正確性，往往需要進(jìn)行大規(guī)模的模擬，這對(duì)于傳統(tǒng)計(jì)算機(jī)集群來(lái)說(shuō)有極大的壓力。而隨著硬件與算法的快速發(fā)展，用來(lái)求解N 體問(wèn)題的軟件也在不斷地更新迭代，如GreeM[1]、CUBE[2]和PHoToNs[3]。

目前模擬N 體問(wèn)題的主要數(shù)值方法有粒子-粒子法（Particle-Particle,P2P），粒子-網(wǎng)格法[4]（Particle-Mesh Method,PM），Barnes-Hut 樹(shù)算法[5]（Barnes-Hut Tree Method,BH-Tree）以及快速多極子算法[6]（Fast Multipole Method,FMM）。P2P 方法的時(shí)間復(fù)雜度為，計(jì)算效率較低，但是精度最高；PM 方法將空間網(wǎng)格離散，通過(guò)快速傅里葉變換（Fast Fourier Transform,FFT）求解泊松方程得到作用力，但近距離時(shí)誤差較大；而樹(shù)方法，包括BH-Tree 算法與FMM 算法，可以控制計(jì)算精度。為了兼顧精度與效率，P3M[7]、TreePM[8]與FMM-PM[9]這樣的混合算法成為了主流。

得益于目前高性能計(jì)算的快速發(fā)展，尤其以CPU+GPU 為主流的并行異構(gòu)平臺(tái)的強(qiáng)大的計(jì)算能力，N 體問(wèn)題成為了加速卡異構(gòu)平臺(tái)的熱點(diǎn)應(yīng)用之一。

國(guó)際上已有不少針對(duì)GPU 這種SIMT（Single Instruction Multiple Thread）并行機(jī)制的加速卡的N體模擬異構(gòu)算法與實(shí)現(xiàn)。Nyland[10]針對(duì)P2P 算法的優(yōu)化提出了“i-parallel”方法以及“walk”的概念，在GPU 上加速了50 倍。而Yokota[11]在Nyland的基礎(chǔ)上，實(shí)現(xiàn)了對(duì)107個(gè)粒子的150 倍加速。由Hamada 提出的“Chamomile Scheme”算法被完全移植到了GPU 上，最終模擬131,072 個(gè)粒子性能達(dá)到了256GFlops[12]，之后通過(guò)對(duì)樹(shù)算法進(jìn)行一系列優(yōu)化，用576 個(gè)GPU 模擬32 億個(gè)粒子，性能達(dá)到了190.5TFlops[13-15]。2011 年，波士頓大學(xué)針對(duì)FMM算法在天文問(wèn)題上的應(yīng)用進(jìn)行GPU 加速[16]，相較于CPU 端加速了70 倍。2017 年，Potter 將其開(kāi)發(fā)的PKDGRAV3[17]中的FMM 算法成功移植到GPU 上，并用4,000 個(gè)GPU 節(jié)點(diǎn)實(shí)現(xiàn)了2 萬(wàn)億粒子的天文模擬。在國(guó)內(nèi)，針對(duì)國(guó)家天文臺(tái)研發(fā)的PHoToNs-2 有相關(guān)的多GPU 的移植與優(yōu)化[18-19]。上海交通大學(xué)以及廈門大學(xué)對(duì)CUBE 的PP 部分進(jìn)行GPU 移植，單步運(yùn)行速度提升了5.16 倍，并在64 個(gè)GPU 上有著83.8%的弱擴(kuò)展效率。

Gadget-2[20]是天體物理領(lǐng)域最有影響力的并行模擬軟件之一，它由Springel 開(kāi)發(fā)，采用TreePM 算法，是國(guó)際上第一個(gè)可計(jì)算粒子規(guī)模達(dá)到百億量級(jí)的開(kāi)源宇宙模擬軟件，在天文領(lǐng)域獲得廣泛認(rèn)可，相關(guān)論文被引用超過(guò)6,500 次。針對(duì)不同的計(jì)算平臺(tái)，Gadget 系列已經(jīng)衍生出了眾多分支版本來(lái)適應(yīng)不同的并行架構(gòu)和編程模式，如OpenAcc-Gadget3[21]、MP-Gadget[22]、AX-Gadget[23]。在2021 年推出的Gadget-4[24]采用了FMM-PM 算法。

本文針對(duì)國(guó)產(chǎn)高性能異構(gòu)加速卡平臺(tái)，基于HIP 完成了Gadget-2 的移植，并對(duì)其在核函數(shù)的訪存效率進(jìn)行了優(yōu)化，最終整體性能加速13.27 倍，短程力計(jì)算加速35.67 倍，并行效率達(dá)到57.29%。第1 節(jié)主要介紹TreePM 以及相關(guān)算法的基本原理；第2 節(jié)介紹加速卡平臺(tái)的計(jì)算環(huán)境，以及最耗時(shí)的短程力計(jì)算在加速卡上的實(shí)現(xiàn)方法；第3 節(jié)介紹短程力計(jì)算在加速卡上的性能優(yōu)化方法；第4 節(jié)驗(yàn)證移植優(yōu)化版本的結(jié)果正確性，并給出性能測(cè)試結(jié)果。

1 TreePM 算法原理

1.1 BH-Tree 算法

BH-Tree 算法是由Barnes 和Hut 提出的，其主要思想是將模擬空間看作是一個(gè)立方體，每一個(gè)維度遞歸地劃分子空間，將空間中的粒子劃分到不同的子節(jié)點(diǎn)，直到每個(gè)子空間只包含了一個(gè)粒子，以此構(gòu)造一個(gè)八叉樹(shù)，并將模擬的整個(gè)空間看作是根節(jié)點(diǎn)。

構(gòu)造樹(shù)之后，從葉子節(jié)點(diǎn)的父節(jié)點(diǎn)開(kāi)始，由下往上依次計(jì)算每個(gè)節(jié)點(diǎn)的質(zhì)點(diǎn)與質(zhì)量，直到根節(jié)點(diǎn)。若節(jié)點(diǎn)只包含一個(gè)粒子，則該節(jié)點(diǎn)的質(zhì)量與質(zhì)心由該粒子的質(zhì)量與位置決定，否則由子節(jié)點(diǎn)決定。在計(jì)算作用力的時(shí)候，自根節(jié)點(diǎn)向下遍歷，粒子之間的相互作用轉(zhuǎn)換為包含多個(gè)粒子的樹(shù)節(jié)點(diǎn)與粒子之間的相互作用，這使得計(jì)算復(fù)雜度降低為，是否能夠?qū)鄠€(gè)粒子的節(jié)點(diǎn)看作一個(gè)整體來(lái)對(duì)目標(biāo)粒子進(jìn)行作用，是通過(guò)以下準(zhǔn)則來(lái)判定：

圖1 BH-Tree 算法中八叉樹(shù)的劃分及粒子受力Fig.1 Partitioning of octrees and the forces acting on the particles in the BH-Tree algorithm

1.2 粒子-網(wǎng)格算法

粒子網(wǎng)格算法是將空間離散網(wǎng)格化，先計(jì)算引力勢(shì)然后再通過(guò)引力勢(shì)計(jì)算其所受作用力。粒子的引力勢(shì)滿足泊松方程[25]：

1.3 TreePM 算法

TreePM[26]算法由Xu 等人提出，其主要思想是由于PM 的運(yùn)算效率比樹(shù)方法快得多，但同時(shí)在小尺度范圍內(nèi)計(jì)算精度較低，而樹(shù)方法可靈活控制精度但計(jì)算大尺度范圍相互作用的效率較低，模擬大規(guī)模N 體問(wèn)題依然存在性能瓶頸，因此綜合權(quán)衡兩個(gè)算法的優(yōu)點(diǎn)和不足之處，TreePM 算法將引力勢(shì)在譜空間中劃分為兩部分：長(zhǎng)程力與短程力，樹(shù)方法計(jì)算短程力，PM 方法計(jì)算長(zhǎng)程力[27]：

實(shí)現(xiàn)時(shí)并不直接計(jì)算公式中的超越函數(shù)，而是通過(guò)樣點(diǎn)查表和插值多項(xiàng)式快速計(jì)算

TreePM 方法與P3M 方法相比誤差近似甚至更小，并且在樹(shù)方法中只需要遍歷目標(biāo)粒子周圍的一個(gè)小空間，兼顧了數(shù)值計(jì)算效率，使得該算法在多個(gè)領(lǐng)域中獲得了廣泛的應(yīng)用，帶來(lái)了巨大的性能提升。

但即便如此，TreePM 的性能瓶頸隨著模擬尺度與粒子規(guī)模的增大依然明顯。經(jīng)過(guò)測(cè)試，PM 部分通過(guò)FFT 計(jì)算，計(jì)算時(shí)間在整個(gè)算法的運(yùn)行過(guò)程中占用比例極小，而樹(shù)方法中計(jì)算短程作用力的部分是計(jì)算密集型的任務(wù)，密度分布較集中時(shí)短程力計(jì)算時(shí)間平均占比90%左右，因此通過(guò)高性能異構(gòu)平臺(tái)對(duì)其進(jìn)行硬件加速，利用多線程計(jì)算來(lái)提高這部分的計(jì)算速度，能夠有效地提升整個(gè)算法的求解效率。

2 TreePM 在加速卡上的實(shí)現(xiàn)

2.1 ROCm 計(jì)算平臺(tái)與HIP 編程模型

“中科先導(dǎo)1 號(hào)”是基于國(guó)產(chǎn)處理器和加速卡的高性能異構(gòu)平臺(tái)，其基本架構(gòu)如圖2 所示，其中每個(gè)節(jié)點(diǎn)有一個(gè)32 核的CPU 處理器以及4 張類GPU 加速卡，分別為Host 端與Device 端。Host 端的內(nèi)存為128GB，每張加速卡有64 個(gè)流多處理器和16GB 顯存，訪存帶寬1TB/s，與CPU 之間使用32 個(gè)PCIE 3.0 互聯(lián)。加速卡與CPU 之間的訪存帶寬是16GB/s，節(jié)點(diǎn)間的網(wǎng)絡(luò)帶寬是25GB/s。CPU 支持x86 指令集以及常用的并行編程模型如OpenMP、MPI 等，加速卡支持ROCm 計(jì)算環(huán)境和HIP 編程模型。

圖2 “中科先導(dǎo)1 號(hào)”架構(gòu)Fig.2 Architecture of Xiandao-1,CAS

ROCm[28-29]開(kāi)源計(jì)算平臺(tái)基于新一代的Radeon硬件，支持最新的數(shù)學(xué)函數(shù)庫(kù)和匯編語(yǔ)言環(huán)境，以及多加速卡的計(jì)算，為開(kāi)發(fā)大規(guī)模、高性能計(jì)算軟件提供了便利，并大大降低了功耗，同時(shí)該平臺(tái)還提供了機(jī)器學(xué)習(xí)框架等多種開(kāi)源技術(shù)，并在不斷地更新，有著優(yōu)秀的可擴(kuò)展性，其中用于可移植性的異構(gòu)計(jì)算接口HIP 是一種C++運(yùn)行時(shí)API 和內(nèi)核語(yǔ)言，它主要是為了簡(jiǎn)化CUDA 應(yīng)用向可移植C++代碼的轉(zhuǎn)換，接口可以跨主機(jī)/內(nèi)核邊界使用模板和類，以及其他C++功能集。移植工具HIPify 可以自動(dòng)完成大部分CUDA 到HIP 的轉(zhuǎn)換工作，代碼在AMD 的DCU 與NVIDIA 的GPU 上通過(guò)各自的編譯器進(jìn)行編譯后便可以正常運(yùn)行，同時(shí)沒(méi)有性能上的損失。

2.2 樹(shù)方法的短程力計(jì)算在加速卡上的實(shí)現(xiàn)

TreePM 方法中，短程力部分采用樹(shù)方法計(jì)算。加速度的計(jì)算是整個(gè)N 體模擬計(jì)算流程最主要的開(kāi)銷，占比可以達(dá)到90%以上，而這一部分中計(jì)算短程力的部分占了極大比例。以表1 為例，粒子數(shù)為1283，采用4 個(gè)進(jìn)程，迭代200 個(gè)時(shí)間步。加速度的計(jì)算耗時(shí)905.1 秒，占總時(shí)間的93.6%，而其中計(jì)算短程力的部分可以占總耗時(shí)的89%，因此對(duì)這一部分進(jìn)行加速，可以大幅度地提升整個(gè)流程的效率。

表1 Gadget-2 的各部分耗時(shí)Table 1 Time consumption for each part of Gadget-2

在Gadget-2 中，樹(shù)方法的實(shí)現(xiàn)是一個(gè)“邊走邊算”的過(guò)程，即獲得目標(biāo)粒子后，從根節(jié)點(diǎn)開(kāi)始，對(duì)目標(biāo)粒子附近的空間節(jié)點(diǎn)進(jìn)行遍歷，獲得相關(guān)的作用節(jié)點(diǎn)，然后計(jì)算兩者的作用力，因此本文樹(shù)方法在加速卡上的實(shí)現(xiàn)中，將對(duì)樹(shù)遍歷的過(guò)程與粒子作用力計(jì)算同時(shí)轉(zhuǎn)移到設(shè)備端進(jìn)行操作，相應(yīng)的我們需要將粒子信息與樹(shù)的相關(guān)信息共同傳入設(shè)備端，即如圖3 所示。

圖3 Tree 方法在加速卡的實(shí)現(xiàn)Fig.3 Implementation of BH-Tree method on accelerators

短程力計(jì)算中，首先在CPU 端將樹(shù)節(jié)點(diǎn)的相關(guān)信息傳入設(shè)備端，再通過(guò)函數(shù)CountParticlesInCurrTime()遍歷粒子列表，得到當(dāng)前時(shí)間步的目標(biāo)粒子數(shù)量以用于分配傳送粒子信息所需要的內(nèi)存。

開(kāi)始對(duì)粒子進(jìn)行引力計(jì)算。對(duì)短程力的計(jì)算包含兩個(gè)部分：local 和export，分別表示當(dāng)前進(jìn)程以及其他進(jìn)程的粒子對(duì)目標(biāo)粒子進(jìn)行的作用力。在計(jì)算local list 時(shí)，遍歷樹(shù)的過(guò)程中會(huì)判定并標(biāo)記當(dāng)前粒子是否會(huì)對(duì)其他進(jìn)程的粒子進(jìn)行作用，當(dāng)計(jì)算完本地進(jìn)程粒子之間的作用力之后，根據(jù)標(biāo)記從其他進(jìn)程獲取需要進(jìn)行相互作用的粒子列表，開(kāi)始計(jì)算export list。實(shí)現(xiàn)時(shí)在函數(shù)force_treeevaluate_shortrange_hip()中調(diào)用核函數(shù)進(jìn)行引力計(jì)算，并傳入?yún)?shù)0 或者1 來(lái)表示計(jì)算local list 或是export list，從而在核函數(shù)開(kāi)始前進(jìn)行對(duì)應(yīng)的數(shù)據(jù)預(yù)處理。最終計(jì)算完成之后對(duì)結(jié)果進(jìn)行匯總，傳回CPU 端后，更新粒子加速度以進(jìn)行下一次迭代。

在原本的Gadget-2 中，是將當(dāng)前時(shí)間步中參與引力計(jì)算的粒子劃分為多次任務(wù)后進(jìn)行計(jì)算，為了減少核函數(shù)的調(diào)用次數(shù)以及數(shù)據(jù)傳輸開(kāi)銷，在計(jì)算第一次任務(wù)時(shí)即可將當(dāng)前進(jìn)程中所有粒子的本地加速度在加速卡上計(jì)算完成，之后的多次任務(wù)只進(jìn)行export list 中的粒子的引力計(jì)算。

核函數(shù)中線程根據(jù)索引取得對(duì)應(yīng)的粒子信息，然后通過(guò)遍歷樹(shù)計(jì)算粒子與節(jié)點(diǎn)之間的作用力。遍歷從根節(jié)點(diǎn)開(kāi)始，其流程與原Gadget-2 基本相同，都是依靠NextNode 數(shù)組來(lái)獲得下一個(gè)需要遍歷的節(jié)點(diǎn)，然后進(jìn)行判定處理，計(jì)算加速度并更新。由于我們需要在遍歷樹(shù)的時(shí)候標(biāo)記粒子是否會(huì)對(duì)其他進(jìn)程中的粒子進(jìn)行作用，此時(shí)粒子在線程中的索引與其在CPU 端Particle 數(shù)組中的索引是不同的，因此還需要將粒子在CPU 端實(shí)際的索引傳到設(shè)備端。Gadget-2 在構(gòu)造樹(shù)時(shí)高效利用了節(jié)點(diǎn)之間的索引關(guān)系，相關(guān)內(nèi)存的消耗大大降低，并且提高了遍歷樹(shù)的效率，同時(shí)也有效地減少了需要傳入設(shè)備端的數(shù)據(jù)量。經(jīng)過(guò)測(cè)試，在模擬規(guī)模達(dá)到1283粒子的整個(gè)過(guò)程中，數(shù)據(jù)傳輸所占用的時(shí)間比例僅是2.5%。

3 核函數(shù)在加速卡上的優(yōu)化

對(duì)于加速卡來(lái)說(shuō)，核函數(shù)的內(nèi)存請(qǐng)求通常是在DRAM 設(shè)備和片上內(nèi)存間以內(nèi)存事務(wù)來(lái)實(shí)現(xiàn)的。連續(xù)對(duì)齊、合并的內(nèi)存訪問(wèn)對(duì)于設(shè)備內(nèi)存來(lái)說(shuō)是最佳訪問(wèn)模式，否則會(huì)因?yàn)閮?nèi)存事務(wù)中有大部分傳輸?shù)臄?shù)據(jù)不會(huì)被使用，而造成帶寬的浪費(fèi)，或者是事務(wù)請(qǐng)求的數(shù)量過(guò)于龐大，增加了訪存壓力，降低了最終的訪存效率。通常來(lái)說(shuō)，有兩種數(shù)據(jù)組織的方式，分別是Struct of Array (SOA)與Array of Struct(AOS)。相比于AOS，SOA 通常能更加有效地減少訪存事務(wù)的數(shù)量，提高訪存效率。

如圖4，采用AOS 模式在加速卡上讀取坐標(biāo)數(shù)據(jù)x 字段，由于隱式加載以及數(shù)據(jù)的交叉存儲(chǔ)訪問(wèn)，每次讀取事務(wù)中只有1/3 的數(shù)據(jù)是我們實(shí)際需要的，同時(shí)會(huì)由于加載了不需要的數(shù)據(jù)而浪費(fèi)緩存空間。如果線程0 到線程5 需要讀取x 軸坐標(biāo)，按照?qǐng)D示讀取方式需要三次事務(wù)請(qǐng)求，而在SOA 模式下讀取數(shù)據(jù)只需要一次事務(wù)請(qǐng)求，更好地利用了加速卡的內(nèi)存帶寬，以達(dá)到更高效的訪存效果。

圖4 SOA 與AOS 內(nèi)存布局Fig.4 The memory layout of SOA and AOS

如算法2 中的偽代碼所示，在核函數(shù)中，在遍歷樹(shù)的過(guò)程中會(huì)反復(fù)用到許多常數(shù)來(lái)判斷節(jié)點(diǎn)之間的關(guān)系，同時(shí)在更新加速度因子fac 時(shí)，會(huì)涉及到對(duì)表的查詢，該表是對(duì)修正系數(shù)的插值。為了提高線程的訪存效率，在不影響并行效率的情況下這些常數(shù)可以先預(yù)取保存在寄存器中，而對(duì)于查詢表可以將其保存在共享內(nèi)存。

優(yōu)化過(guò)程中也可以采用頁(yè)鎖定內(nèi)存以及異步多流的方式來(lái)提高數(shù)據(jù)在CPU 端與設(shè)備端之間的傳輸效率，但是如2.2 節(jié)所述，數(shù)據(jù)傳輸占據(jù)的整體運(yùn)行時(shí)間比例僅為2.5%，因此該類方法并不能大幅提升整體性能。

4 數(shù)值結(jié)果與正確性

4.1 正確性驗(yàn)證

下面將Gadget-2 的HIP 移植優(yōu)化版本記作G2H。實(shí)驗(yàn)平臺(tái)為“中科先導(dǎo)1 號(hào)”，軟硬件環(huán)境詳見(jiàn)2.1 節(jié)。宇宙模擬的初始條件由基于二階拉格朗日擾動(dòng)理論和Lambda Cold Dark Matter（LCDM）模型編寫的軟件2LPTic 生成。LCDM 模型是宇宙大爆炸理論的數(shù)學(xué)參數(shù)化描述，該模型假設(shè)宇宙中的大部分物質(zhì)由冷暗物質(zhì)組成，它們只存在引力相互作用。這是如今宇宙學(xué)中普遍認(rèn)同的標(biāo)準(zhǔn)化模型，也是目前可以為觀察到的宇宙微波背景輻射、宇宙的大尺度結(jié)構(gòu)等現(xiàn)象提供合理解釋的最簡(jiǎn)單的模型，同時(shí)在該宇宙模型中有一個(gè)與暗能量相關(guān)的常數(shù)Lambda，它能夠解釋哈勃定律，即宇宙隨時(shí)間不斷加速膨脹的過(guò)程。

文中用功率譜來(lái)驗(yàn)證G2H 的計(jì)算正確性。驗(yàn)證方法是通過(guò)計(jì)算統(tǒng)計(jì)暗物質(zhì)密度分布，然后根據(jù)模擬尺度得到密度分布對(duì)應(yīng)的功率譜。功率譜是隨機(jī)振幅的兩點(diǎn)相關(guān)函數(shù)的傅里葉變換，這個(gè)振幅會(huì)在引力的影響下隨著時(shí)間的推移而緩慢變化，因此我們可以通過(guò)比較程序模擬中相同紅移處的功率譜檢驗(yàn)正確性。

實(shí)驗(yàn)過(guò)程中模擬1283個(gè)粒子在尺寸為100h-1Mpc，以共動(dòng)坐標(biāo)系為參考的盒子中的運(yùn)動(dòng)狀態(tài)，并考慮周期邊界條件，初始紅移z=49，其他參數(shù)包括：同時(shí)采用PM 方法，網(wǎng)格參數(shù)設(shè)置為PMGRID=128。實(shí)驗(yàn)過(guò)程中選取紅移z=0、1、2、3、5 來(lái)進(jìn)行功率譜的比對(duì)。

Gadget-2 和G2H 在紅移z=0 時(shí)粒子分布的柱密度圖如圖5 所示，圖6 的功率譜中線條和符號(hào)分別表示G2H 和Gadget-2 的結(jié)果，兩者的殘差幾近為0，保證了G2H 的正確性。

圖5 Gadget-2（左）與G2H（右）的N 體模擬在紅移z=0 時(shí)的柱密度圖Fig.5 The density map for the N-body simulation of Gadget-2(left panel) and G2H (right panel) as red shift is 0

圖6 Gadget-2 與G2H 的功率譜比較Fig.6 Comparison of power spectrum of Gadget-2 and G2H

4.2 性能測(cè)試

下面分別測(cè)試加速比和并行效率。由于計(jì)算過(guò)程中每一步粒子的分布不同，因此每個(gè)時(shí)間步的運(yùn)行時(shí)間有波動(dòng)，我們進(jìn)行10 步迭代然后取平均值作為單步運(yùn)行時(shí)間進(jìn)行比較。粒子規(guī)模分別取1283、2563、5123，分別使用4、8、16 個(gè)進(jìn)程(記作Np)與相同數(shù)量的加速卡對(duì)Gadget-2、G2H 進(jìn)行性能測(cè)試，平均單步的總運(yùn)行時(shí)間和短程力運(yùn)行時(shí)間見(jiàn)表2。

由表2 可以得到，粒子規(guī)模越大，程序在加速卡上的加速效果越明顯，因?yàn)槎坛塘Φ挠?jì)算時(shí)間占比一般會(huì)和粒子規(guī)模一同增加。相比于Gadget-2，G2H 程序的整體加速比可以從11.47 倍提高到13.27倍，短程力部分計(jì)算加速比可以從24.25 倍提高到35.67 倍。

表2 Gadget-2、G2H 的單步平均耗時(shí)Table 2 Average step runtime of Gadget-2 and G2H

在以往PHoToNs 的移植優(yōu)化工作中，短程力計(jì)算效率可以提升上百倍，而Gadget-2 中熱點(diǎn)函數(shù)在加速卡上移植優(yōu)化后性能只能提升30 倍左右，這與兩者的計(jì)算流程有關(guān)。在PHoToNs 中提前進(jìn)行樹(shù)的遍歷，將相互作用的粒子放入包中，而加速卡提取對(duì)應(yīng)的包進(jìn)行計(jì)算，沒(méi)有過(guò)多的其他步驟，但在Gadget-2 中由于需要通過(guò)遍歷樹(shù)不斷判定來(lái)進(jìn)行引力對(duì)的計(jì)算，除計(jì)算外的冗余步驟以及邏輯判斷較多，這些都嚴(yán)重影響了G2H 中核函數(shù)的并行效率，導(dǎo)致兩者的優(yōu)化效果有著明顯的差異。

在并行效率測(cè)試中，粒子規(guī)模為5123，分別使用8、16、32、64、128 個(gè)進(jìn)程與相同數(shù)量的加速卡迭代10 步后取平均值作為單步耗時(shí)。圖7 顯示了G2H 在“中科先導(dǎo)一號(hào)”上的強(qiáng)可擴(kuò)展性，表3 是測(cè)試中對(duì)應(yīng)的具體耗時(shí)?？梢钥闯?，在加速卡上不僅有短程力的計(jì)算，還有樹(shù)的遍歷，這對(duì)短程力在加速卡上的并行效率有著較大的影響，而在主機(jī)端實(shí)現(xiàn)了非本地樹(shù)節(jié)點(diǎn)之間的通信以及PM 的通信，這些因素也成為了程序主要的性能瓶頸。最終G2H 的整體并行效率是57.29%，其中短程力并行效率是53.87%。

表3 G2H 強(qiáng)可擴(kuò)展性測(cè)試的單步平均耗時(shí)Table 3 Average step runtime of G2H in strong scalability test

圖7 G2H 的強(qiáng)可擴(kuò)展性測(cè)試Fig.7 The strong scalability test of G2H

5 結(jié)論與展望

本文將基于TreePM 算法的宇宙學(xué)N 體模擬軟件Gadget-2 中耗時(shí)占比90%以上的短程力計(jì)算模塊通過(guò)HIP 移植到加速卡異構(gòu)平臺(tái)“中科先導(dǎo)一號(hào)”。通過(guò)數(shù)據(jù)內(nèi)存重構(gòu)，并充分利用寄存器和共享內(nèi)存，減少了加速卡對(duì)內(nèi)存存取開(kāi)銷，提高了設(shè)備端訪存效率。通過(guò)功率譜驗(yàn)證了移植優(yōu)化的正確性，并進(jìn)行相關(guān)的性能測(cè)試，程序整體與短程力部分的計(jì)算速度最高分別提升了13.27 倍與35.67 倍，并行效率分別達(dá)到57.29%與53.87%。由于本文的方案是在加速卡上進(jìn)行樹(shù)的遍歷，導(dǎo)致較多的邏輯判斷，這對(duì)加速卡性能的充分挖掘是不利的，因此考慮將樹(shù)的遍歷盡可能地在CPU 上進(jìn)行，在加速卡端進(jìn)行主要的計(jì)算，同時(shí)這樣也能平衡主機(jī)與設(shè)備端之間的數(shù)據(jù)傳輸，這些將作為下一步的研究任務(wù)。本文的移植優(yōu)化方法可為類似的樹(shù)結(jié)構(gòu)算法在高性能異構(gòu)平臺(tái)上的實(shí)現(xiàn)與大規(guī)模異構(gòu)并行模擬提供借鑒。

利益沖突聲明

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