丁丹迪，文敏華，周 姍，陳 民，林新華

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

2.英特爾亞太研發(fā)中心，上海 200240

3.上海交通大學(xué) 物理與天文系，上海 200240

4.東京工業(yè)大學(xué)，日本 東京

1 引言

隨著高性能計(jì)算領(lǐng)域的快速發(fā)展，激光等離子體模擬程序可以在更大維度以及更高的精度下進(jìn)行。“虛擬激光等離子體實(shí)驗(yàn)室”程序（virtual laser plasma lab，VLPL）[1]是一份基于C++開發(fā)的三維激光等離子體模擬程序，被廣泛應(yīng)用于激光與等離子體相互作用的研究中[2-4]。VLPL-S由上海交通大學(xué)陳民研究員在VLPL基礎(chǔ)上進(jìn)一步開發(fā)，優(yōu)化了碰撞電離、場致電離以及兩體碰撞效應(yīng)和粒子跟蹤，并提供了對眾核平臺的支持。在代碼中，電子的離化是使用虛函數(shù)開發(fā)實(shí)現(xiàn)的，體現(xiàn)了面向?qū)ο缶幊陶Z言的特性；移動窗口方法也應(yīng)用到了程序中以研究加速能量達(dá)GeV的長距離激光尾波場加速問題；通過加入以上的模擬方法，程序可以對Terahertz輻射、中性氣體與激光相互作用中的光脈沖及離化波前演化、激光尾波場加速及表面電子加速等物理過程進(jìn)行研究[5]。

在本文中，VLPL-S使用多種方法優(yōu)化，并在Intel?第二代Xeon PhiTM處理器（代號Knights Landing）上進(jìn)行性能評估。由于原程序的熱點(diǎn)受限于訪存延遲，預(yù)取優(yōu)化方法顯著提升了程序在Knights Landing（KNL）上的性能。通過應(yīng)用多線程優(yōu)化，VLPL-S在KNL集群上的性能也得以提升。使用這些優(yōu)化方法后，VLPL-S相比于原版程序獲得了1.68倍的加速比。

本文組織結(jié)構(gòu)如下：第2章介紹相關(guān)工作；第3章和第4章簡要介紹Particle-in-Cell（PIC）算法和Intel KNL架構(gòu)；第5章講述詳細(xì)的優(yōu)化方法；第6章對優(yōu)化后的性能進(jìn)行評估；第7章給出結(jié)論以及后續(xù)優(yōu)化思路。

2 相關(guān)工作

在最近高性能計(jì)算快速發(fā)展的三十年內(nèi)，許多學(xué)者提出了各自PIC算法的實(shí)現(xiàn)。這些實(shí)現(xiàn)從側(cè)面說明了PIC算法的并行可能性。Bastrakov等人[6]最近講述了對其實(shí)現(xiàn)的PIC代碼在多核平臺、KNC（Knights Corner，Intel第一代Xeon Phi眾核處理器）以及GPGPU（general purpose GPU）上的性能評估。通過他們的優(yōu)化，程序的數(shù)據(jù)局部性得以提升，并能更充分利用處理器SIMD（single instruction multiple data）指令，最終在一個8核的Xeon E5-2690平臺上實(shí)現(xiàn)了7倍的加速比，以及在共享內(nèi)存集群上實(shí)現(xiàn)了99%的強(qiáng)擴(kuò)展性。對于其代碼在KNC上的實(shí)現(xiàn)，Bastrakov等人僅在7100X協(xié)處理器上實(shí)現(xiàn)了相對于8核Xeon E5-2690 CPU 1.5倍的加速比。其程序在KNC上的性能僅相當(dāng)于Xeon Phi 7100X峰值性能的20%。與此不同，在基于Fermi架構(gòu)的GPGPU平臺上他們通過原子操作以減少場部分?jǐn)?shù)據(jù)讀寫，最終單精度版本的程序?qū)崿F(xiàn)了相對于8核Xeon E5-2690 CPU 10倍以上的加速比。

此外，Surmin等人[7]開發(fā)了一套KNC平臺上的高性能PIC代碼，代碼名稱是PICADOR。優(yōu)化后，程序的數(shù)據(jù)局部性和向量化程度得以提高，最終在Intel Xeon E5-2600 CPU上實(shí)現(xiàn)了3.75倍的加速比，在Xeon Phi 5110P（KNC）協(xié)處理器上實(shí)現(xiàn)7.5倍的加速比。PICADOR在Xeon Phi 5110P上的性能是其在單Xeon CPU上性能的1.6倍。Surmin等學(xué)者還評估了PICADOR在KNL上的性能[8]，僅重新針對KNL平臺編譯代碼就獲得了相對于KNC上2.43倍的加速比。在優(yōu)化后，PIDADOR又獲得了1.89倍的加速比，并且實(shí)現(xiàn)了100 GFLOPS的雙精度峰值性能。

VLPL-S常用于離化相關(guān)的模擬流程，因此以相較于其他的PIC實(shí)現(xiàn)，VLPL-S在模擬過程中會產(chǎn)生更多的新粒子，從而此程序需要消耗更多的內(nèi)存，并且粒子移動也會產(chǎn)生更多的通信開銷。盡管現(xiàn)在已經(jīng)處于P級計(jì)算能力的時代，想要在眾核平臺上實(shí)現(xiàn)極致的性能也是不容易的。Nakashima[9]也在其文章中展示了在眾核平臺上實(shí)現(xiàn)高度并行及向量化的PIC程序的困難性。好消息是Intel的新一代Xeon Phi處理器提供了高達(dá)6 TFLOPS的峰值計(jì)算能力與由MCDRAM（multi-channel dynamic random access memory）帶來的極高訪存帶寬，并且也同時提供了與CPU平臺統(tǒng)一的編程環(huán)境，因此在KNL這個新平臺上移植并優(yōu)化VLPL-S程序?qū)且豁?xiàng)有意義的研究工作。

3 VLPL-S介紹

本部分將簡要介紹VLPL及VLPL-S程序中所使用的算法，PIC算法的介紹可參見文獻(xiàn)[10]。

程序中的模擬空間是一個由均勻的網(wǎng)格所劃分的三維空間。麥克斯韋方程組所定義的電場強(qiáng)度E和磁場強(qiáng)度B儲存在網(wǎng)格當(dāng)中，并用有限時域差分法（finite difference time domain，F(xiàn)DTD）求解[7]。粒子以有限大小的粒子云模型表示，眾多粒子分布由初始等離子體條件給定，單個粒子云在算法中具體體現(xiàn)為動量p、位置r、質(zhì)量m以及電荷q。粒子的位置r和速度v根據(jù)洛倫茲力以及Boris數(shù)值方程求解，進(jìn)而得出粒子的移動量。帶電粒子移動會產(chǎn)生電流j，是麥克斯韋方程組中的源項(xiàng)。

Fig.1 Stages and data dependency in PIC algorithm圖1 PIC算法的流程及數(shù)據(jù)依賴

PIC算法中包含4個主要的處理過程，每個過程所涉及到的方程以及數(shù)據(jù)依賴可參見圖1。主循環(huán)的一次執(zhí)行表示一個時間步的模擬。每個時間步中包含了上述的4個過程。場計(jì)算（field solver）部分會更新網(wǎng)格中儲存的電磁場強(qiáng)度；場插值及力計(jì)算部分會計(jì)算帶電粒子在電磁場中所受的洛倫茲力；粒子推動部分會根據(jù)粒子的動力學(xué)數(shù)據(jù)更新粒子位置并在網(wǎng)格中移動粒子；最后的部分是統(tǒng)計(jì)各個網(wǎng)格中由粒子運(yùn)動所產(chǎn)生的電流。從代碼實(shí)現(xiàn)角度來看，將部分場強(qiáng)計(jì)算、粒子受力及推動部分整合到一起可以提高數(shù)據(jù)局部性。整合后的部分在VLPL-S代碼move_particle()函數(shù)中實(shí)現(xiàn)。VLPL-S代碼中還考慮了激光與等離子體的相互作用，因此需要一個額外的流程來處理激光在電磁場中所受的影響。在PIC算法的所有處理流程中，粒子推動（particle push）部分耗時最多，因此move_particle()也是整個VLPLS程序的熱點(diǎn)。

關(guān)于空間劃分方法，VLPL-S在三維空間的每個維度都進(jìn)行了網(wǎng)格均分。VLPL-S中每個網(wǎng)格的大小都是相同的，每個 MPI（message passing interface）進(jìn)程包含相同數(shù)量的網(wǎng)格。由于采用了這樣的空間劃分方法，只有在模擬粒子的空間分布十分均勻的算例時才能取得良好的負(fù)載均衡性。

PIC方法中包含兩種不同的數(shù)據(jù)集。首先是空間中的粒子數(shù)據(jù)，其次是網(wǎng)格上的電磁場數(shù)據(jù)及電流、粒子密度等數(shù)據(jù)[7]。由于粒子推動及電流計(jì)算部分有復(fù)雜的計(jì)算過程以及對兩個數(shù)據(jù)集的讀寫，它們占去了大部分的運(yùn)行時間。VLPL-S程序中粒子由鏈表數(shù)據(jù)結(jié)構(gòu)表示。每個網(wǎng)格維護(hù)一條單向鏈表，其中儲存了此網(wǎng)格中的所有粒子。由鏈表所產(chǎn)生的隨機(jī)訪存對性能的影響將在之后討論。

4 Intel Knights Landing

Knights Landing是Intel?Xeon PhiTM眾核處理器的第二代產(chǎn)品，它具有多個新的特性。首先KNL支持自啟動，也就是說它可以作為一個獨(dú)立的處理器來使用而不僅僅是協(xié)處理器。其次，藉由其片上整合的MCDRAM，KNL提供了400+GB/s的訪存帶寬，同時還包含傳統(tǒng)的DDR4內(nèi)存以提供大內(nèi)存容量支持。MCDRAM有不同的使用方式：cache、flat和hybrid。不同的方式有不同的地址空間映射以及物理內(nèi)存劃分方法。此外KNL還有不同的訪存路由處理模擬：All-to-All、Quadrant、SNC（sub NUMA clustering）。關(guān)于KNL的詳細(xì)介紹可參見文獻(xiàn)[11-12]。詳細(xì)的帶寬數(shù)據(jù)[13]以及和雙路Xeon節(jié)點(diǎn)帶寬的對比可見表1。

Table 1 Measured stream bandwidth of Xeon and Xeon Phi表1 Xeon和Xeon Phi的實(shí)測帶寬 GB/s

KNL的芯片由36個Tile組成。這些Tile由二維片上網(wǎng)格所連接。每個Tile中有兩個核以及1 MB的共享L2緩存。每核中有兩個向量處理單元。詳細(xì)的KNL架構(gòu)可見圖2（b）。OPA（Intel?Omini-Path architecture）也集成到KNL上以提供高速互連。

5 優(yōu)化方法

本部分首先介紹VLPL-S程序的關(guān)鍵特性以及性能瓶頸，然后詳細(xì)講述本文所采用的優(yōu)化方法。

5.1 基礎(chǔ)優(yōu)化

Fig.2 Knights Landing overview圖2 Knights Landing架構(gòu)

通常來說，平方根操作相比其他操作需要更多的時鐘周期，并且對512位寬數(shù)據(jù)的平方根也會產(chǎn)生更大的開銷。原版的VLPL-S程序基本使用的是單精度浮點(diǎn)數(shù)，但是由于編譯器在混合精度表達(dá)式中的自動精度擴(kuò)展，發(fā)現(xiàn)由Intel編譯產(chǎn)生的匯編程序中，部分單精度浮點(diǎn)數(shù)在平方根操作時被擴(kuò)展成為雙精度浮點(diǎn)數(shù)。精度擴(kuò)展的初衷是保證數(shù)據(jù)的精確性，因此單精度與雙精度的復(fù)合表達(dá)式中的單精度數(shù)據(jù)被自動擴(kuò)展成為雙精度數(shù)據(jù)。這在VLPL-S程序中不是必須的，因此將部分雙精度數(shù)據(jù)強(qiáng)制轉(zhuǎn)換成為單精度數(shù)據(jù)，從而可以減少部分平方根操作產(chǎn)生的數(shù)據(jù)等待。

虛函數(shù)在C++中是由虛表所實(shí)現(xiàn)的，它需要在運(yùn)行時根據(jù)內(nèi)存地址來決定函數(shù)調(diào)用的入口，而函數(shù)調(diào)用所產(chǎn)生的寄存器讀寫又是十分不利于性能的[14]。另外，自動內(nèi)聯(lián)是編譯器提供的一個對高性能計(jì)算程序十分有用的特性。自動內(nèi)聯(lián)可以減少核心計(jì)算時上下文切換所帶來的開銷。然而自動內(nèi)聯(lián)并不能作用于虛函數(shù)，因此將程序中所使用到的虛函數(shù)都進(jìn)行工程重寫，改為可以被自動內(nèi)聯(lián)的普通函數(shù)。

5.2 訪存優(yōu)化

在對VLPL-S的介紹中提到，程序中的粒子數(shù)據(jù)是使用鏈表的數(shù)據(jù)結(jié)構(gòu)進(jìn)行儲存的。PIC算法中每個時間步都需要處理空間中所有的粒子，因此模擬過程中對鏈表的多次遍歷就是不可避免的。鏈表數(shù)據(jù)結(jié)構(gòu)中每個節(jié)點(diǎn)的內(nèi)存是單獨(dú)分配的，因此這些節(jié)點(diǎn)處于內(nèi)存中的不同位置。由此造成遍歷鏈表時的低效隨機(jī)訪存。

預(yù)取是一個用于減少訪存延遲影響的常用優(yōu)化方法。硬件預(yù)取與軟件預(yù)取共同工作以給處理器的計(jì)算單元提供數(shù)據(jù)[15]。然而，目前的硬件預(yù)取很難對遍歷鏈表時的隨機(jī)訪存進(jìn)行預(yù)測，即便使用Intel ICC編譯器也不能對其做自動軟件預(yù)取優(yōu)化。為了取得更高的訪存性能，使用Intrinsic庫將鏈表中下一個要處理的粒子數(shù)據(jù)預(yù)取到L1緩存。由于鏈表的訪存特性，在讀取到當(dāng)前節(jié)點(diǎn)的數(shù)據(jù)之前并不知道下一個節(jié)點(diǎn)的數(shù)據(jù)位置，因此單純的預(yù)取也不能完全隱藏鏈表訪存延遲。為此，在每個網(wǎng)格中維護(hù)了一個前向查找表，以便在計(jì)算每個網(wǎng)格時預(yù)取所需要用到的粒子數(shù)據(jù)。這樣就可以提前預(yù)取數(shù)據(jù)以更好地隱藏隨機(jī)訪存延遲。

5.3 多線程優(yōu)化

與Xeon相比較，盡管KNL的主頻要低一些，但KNL擁有更多的核心。更多核心需要運(yùn)行更多的進(jìn)程才能全部利用其計(jì)算能力，因此要想在KNL上實(shí)現(xiàn)極致的性能，就必須實(shí)現(xiàn)極致的并行性。

VLPL-S程序的熱點(diǎn)是PIC算法中的粒子推動函數(shù)，它占用了程序執(zhí)行時間的80%以上。其中由粒子運(yùn)動產(chǎn)生的電流數(shù)據(jù)需要更新到網(wǎng)格當(dāng)中，網(wǎng)格中的不同粒子可能會同時更新同一個位置的數(shù)據(jù)。因此，在多線程并行化的同時需要處理數(shù)據(jù)沖突的問題。

在每個時間步的迭代之后，粒子可能會移動到屬于其他進(jìn)程的網(wǎng)格中。每個進(jìn)程邊界上的電磁場數(shù)據(jù)也需要和其周圍的進(jìn)程進(jìn)行同步更新，因此每個時間步之后需要同步一次數(shù)據(jù)。由此，整個程序運(yùn)行過程中會產(chǎn)生大量的MPI通信。在Intel ITAC工具的幫助下，發(fā)現(xiàn)當(dāng)MPI進(jìn)程增加時MPI通信開銷所占的時間比也在增加。此外，負(fù)載均衡性較差的算例會產(chǎn)生更多的MPI通信開銷。此問題的原因是進(jìn)程之間計(jì)算量不同而每個時間步都需要同步，因此計(jì)算量較少的進(jìn)程會阻塞在MPI通信函數(shù)等待其他進(jìn)程。為了降低通信帶來的開銷以及在MPI進(jìn)程內(nèi)提升負(fù)載均衡性，本文使用OpenMP對VLPL-S程序進(jìn)行了多線程優(yōu)化。

如上文所述，一個粒子的移動會對其周圍多個網(wǎng)格的電流產(chǎn)生影響。不同線程中的粒子可能會同時對一個網(wǎng)格的電流做出貢獻(xiàn)，因而產(chǎn)生了數(shù)據(jù)寫沖突問題。本文提出3種方法來處理數(shù)據(jù)寫沖突問題。

（1）原子操作

OpenMP庫提供了原子操作語句，其指定的內(nèi)存地址同時只能由一個線程操作，這就可以避免數(shù)據(jù)寫沖突。具體來說OpenMP提供了一個最小化的關(guān)鍵區(qū)塊以實(shí)現(xiàn)原子操作[16]。從PIC算法上來說，每個粒子只需要對其周圍的網(wǎng)格做出貢獻(xiàn)，因此只需在每個線程的邊界區(qū)域應(yīng)用Atomic操作即可。應(yīng)用此方法后的線程網(wǎng)格示意圖如圖3（a）所示。

Fig.3 Scheme of methods avoiding data hazard of thread level parallelism圖3 避免數(shù)據(jù)寫沖突的處理方法示意圖

（2）被動貢獻(xiàn)方法

在原版的VLPL-S代碼當(dāng)中，粒子移動所產(chǎn)生的電流數(shù)據(jù)被直接寫入到對應(yīng)的網(wǎng)格當(dāng)中。對于粒子而言，這種寫數(shù)據(jù)方法是主動貢獻(xiàn)法。在被動貢獻(xiàn)方法中，電流在產(chǎn)生之后先會保存在粒子自己的數(shù)據(jù)結(jié)構(gòu)體之中，當(dāng)所有粒子的電流計(jì)算完成后再從網(wǎng)格的粒度遍歷數(shù)據(jù)，以將電流從粒子中讀取并統(tǒng)計(jì)寫入到網(wǎng)格當(dāng)中。如上而言，每個粒子需要一片緩存空間來保存產(chǎn)生的電流數(shù)據(jù)。被動貢獻(xiàn)方法的示意圖如圖3（b）所示。

（3）緩存規(guī)約方法

這個方法是在被動貢獻(xiàn)方法的基礎(chǔ)上提出的。為了保存粒子移動所產(chǎn)生電流的中間數(shù)據(jù)，需要在粒子結(jié)構(gòu)體中增加3個32位的變量。VLPL-S模擬中通常包含了百萬以上個粒子，這些新增的數(shù)據(jù)空間將會使整個程序所需要的內(nèi)存空間急劇增加，同時也會產(chǎn)生更多的訪存需求。緩存規(guī)約方法在被動貢獻(xiàn)方法基礎(chǔ)上做出了改進(jìn)，將電流數(shù)據(jù)暫存到網(wǎng)格數(shù)據(jù)中，這樣只需要在每個網(wǎng)格中新增9個變量的數(shù)據(jù)空間即可保存電流數(shù)據(jù)。相較于被動貢獻(xiàn)法，緩存規(guī)約方法所需的內(nèi)存空間更低。圖4為緩存規(guī)約方法的示意圖。

Fig.4 Scheme of buffering and reduction for avoiding data hazard圖4 應(yīng)對數(shù)據(jù)寫沖突的緩存規(guī)約方法

在MPI層面，不同的進(jìn)程擁有相同數(shù)量的網(wǎng)格。由于粒子在網(wǎng)格之間的頻繁移動，粒子在網(wǎng)格中的分布可能逐步變得不均勻。在OpenMP層面，水平和垂直方向的嵌套循環(huán)可以展開成為一個循環(huán)，因而可以使用OpenMP動態(tài)調(diào)度方法來處理粒子在進(jìn)程內(nèi)的分布不均衡問題。

5.4 向量化

KNL提供的一個重要特性是每個核心擁有兩個512位的向量處理單元。本節(jié)將講述一個對于鏈表數(shù)據(jù)結(jié)構(gòu)的中間向量化方法。對于此方法的評估及討論將在第6章進(jìn)行。

鏈表數(shù)據(jù)結(jié)構(gòu)是針對有頻繁插入與刪除操作的數(shù)據(jù)集所設(shè)計(jì)的，它的優(yōu)點(diǎn)是任意節(jié)點(diǎn)的插入與刪除的時間復(fù)雜度都是O(1)。然而，遍歷鏈表時的隨機(jī)內(nèi)存訪問卻是高性能計(jì)算應(yīng)用所不希望看到的。

為了充分利用512位寬的向量化處理單元以及對向量化所帶來的性能提升有一個初步的認(rèn)識，本文采用了工作量較小但有一些額外開銷的臨時向量化方法。此方法的簡要示意圖見圖5。在原程序中，盡管每個粒子的質(zhì)量、電荷及動量等數(shù)據(jù)儲存在連續(xù)的內(nèi)存空間，但是4個相鄰粒子的質(zhì)量數(shù)據(jù)卻是離散的，這樣的內(nèi)存布局不能直接向量化。此臨時向量化方法會一次讀取數(shù)個粒子的數(shù)據(jù)，并將同類數(shù)據(jù)拼裝打包到連續(xù)的內(nèi)存當(dāng)中，這樣就可以組成512位的向量化數(shù)據(jù)。VLPL-S程序使用的是32位單精度浮點(diǎn)數(shù)，為了組成512位的向量化數(shù)據(jù)，每次計(jì)算之前需要讀取的粒子數(shù)量由1個變?yōu)榱?6個。

Fig.5 An example of packing 4 linked particles圖5 拼裝組合4個粒子數(shù)據(jù)的示意圖

6 結(jié)果討論

6.1 算例及測試平臺說明

本文使用兩個不同的算例來進(jìn)行性能測試。表2展示了詳細(xì)的算例配置。其中，算例A是一個主要用戶性能測試的算例，因?yàn)榱Ｗ拥姆植际窍鄬鶆虻?，所以算例A有較好的負(fù)載均衡性。此外，算例A模擬的空間在水平和垂直方向上的網(wǎng)格數(shù)量是相當(dāng)?shù)?，這對于大多數(shù)的空間劃分方法都是適用的。算例B是一個用于實(shí)際研究的算例。初始化時，粒子在空間中的分布是完全均衡的，但隨著模擬的進(jìn)行，粒子的分布變得不均衡。此外，算例B的模擬空間是一個高度較低的空間，這導(dǎo)致了在空間劃分時垂直方向最多只能劃分為60個。

Table 2 Details about workloads表2 測試所用算例說明

在3個平臺上對VLPL-S的性能進(jìn)行了測試，表3展示了3個平臺的詳細(xì)配置。VLPL-S在KNL上使用了多種的Memory mode和Cluster mode進(jìn)行了測試，詳細(xì)的配置將隨性能數(shù)據(jù)一起給出。測試時KNL平臺的代碼編譯選項(xiàng)是“-O3-xMIC-AVX512”，Xeon平臺的編譯選項(xiàng)是“-O3-xCORE-AVX2”。

Table 3 Machine configuration表3 測試平臺配置說明

6.2 在KNL上對原版VLPL-S代碼的性能評估

KNL一個重要的新特性是對CPU平臺程序的二進(jìn)制兼容性。為了了解KNL平臺對原生CPU程序的支持及其性能，參照Intel編譯器團(tuán)隊(duì)的指導(dǎo)，僅使用“-xMIC-AVX512”對原版程序進(jìn)行了重新編譯。圖6展示了VLPL-S在KNL的不同Memory mode及Cluster mode上的性能。

當(dāng)使用MCDRAM時，VLPL-S在不同的Memory mode及Cluster mode下的性能基本是一致的。而僅使用DDR內(nèi)存則性能會差很多。使用MCDRAM時，程序從128個進(jìn)程擴(kuò)展到256個進(jìn)程的加速比約為1.5，如果使用DDR內(nèi)存則基本沒有加速比。原因是MCDRAM能提供相比于DDR內(nèi)存更高的內(nèi)存帶寬。程序使用256個進(jìn)程運(yùn)行時，其性能已受限于訪存帶寬。

6.3 優(yōu)化方法的效果評估

Fig.6 Apreview of performance of VLPL-S on KNL圖6 原版VLPL-S程序在KNL上的性能

Fig.7 Performance gain of using different optimizations on Xeon and Xeon Phi platforms圖7 使用不同優(yōu)化方法在Xeon及Xeon Phi平臺上的加速比

第5章提到的優(yōu)化方法所取得的性能提升如圖7所示。KNL平臺使用了Flat memory mode以及Quadrant cluster mode。VLPL-S在Xeon Phi 7210平臺上運(yùn)行了256個進(jìn)程，在雙路Xeon E5-2699v3平臺上運(yùn)行了36個進(jìn)程。程序在KNL上使用256個進(jìn)程時性能是最佳的，其原因是MCDRAM能夠提供足夠的訪存帶寬。KNL平臺每個核心有兩個VPU，單精度的標(biāo)量計(jì)算程序?qū)PU不會產(chǎn)生極高的壓力，因此VLPL-S可以在超線程上獲得性能提升。對于CPU平臺而言，其內(nèi)存帶寬是相對有限的，同時每個核心僅有一個VPU，因此CPU上運(yùn)行36個進(jìn)程時性能是最佳的。

（1）基礎(chǔ)優(yōu)化

如圖7所示，基礎(chǔ)優(yōu)化方法在KNL平臺上的加速比比CPU平臺上的高。原因是KNL主頻較低，上下文切換的開銷更大，因此消除虛函數(shù)的優(yōu)化在KNL上更加有效。此外，算例A在基礎(chǔ)優(yōu)化后的性能提升高于算例B的性能提升。原因是粒子在算例A中移動得較少，因此通信所占的時間更少，使得主要作用于計(jì)算部分的優(yōu)化更加有效。

（2）訪存優(yōu)化

為了評估預(yù)取的效果，使用Intel VTuneTM收集了運(yùn)行時的硬件計(jì)數(shù)器數(shù)據(jù)。如圖8所示（測試平臺為雙路E5-2699v3），緩存缺失率下降了。預(yù)取優(yōu)化對CPU和KNL都十分有效。對算例A而言，其他KNL平臺上的加速比是1.21倍，這比它在Xeon平臺上的要高。原因是KNL平臺上運(yùn)行了8倍于CPU平臺的進(jìn)程數(shù)量，MCDRAM提供了更高的訪存帶寬。對于算例B而言，由于其負(fù)載不均衡性，預(yù)取主要作用于相對任務(wù)較重的進(jìn)程，從而不能有效提升整體程序的帶寬。因此算例A在KNL平臺上獲得了更高的加速比。

Fig.8 Cache miss rate reduced圖8 緩存缺失率下降

（3）多線程優(yōu)化

應(yīng)用多線程優(yōu)化后，VLPL-S在單節(jié)點(diǎn)上取得了相當(dāng)?shù)男阅芴嵘?。對算例A而言，程序在KNL平臺取得了1.14倍的加速比，在CPU平臺上取得了1.05倍的加速比。對于算例B而言，程序在KNL平臺上實(shí)現(xiàn)了1.27倍的加速比，在CPU平臺上實(shí)現(xiàn)了1.11倍的加速比。兩個算例在KNL平臺上所獲得的加速比都更高。在未使用多線程優(yōu)化的情況下，VLPL-S最佳的運(yùn)行配置是KNL平臺256個進(jìn)程，CPU平臺36個進(jìn)程。多線程優(yōu)化僅在程序的熱點(diǎn)應(yīng)用，在超線程的幫助下它的硬件過度使用開銷相對更低。此外，算例B在兩個平臺上所取得的加速比都比算例A要高。其原因是多線程優(yōu)化后的線程動態(tài)調(diào)度策略可以在線程層面上緩解算例B的負(fù)載不均衡問題。

MPI+OpenMP的混合使用模式極大提升了VLPLS在多個KNL節(jié)點(diǎn)上的擴(kuò)展性。擴(kuò)展性在使用OPA互連的Intel Endeavor集群上測試得出，KNL和BDW節(jié)點(diǎn)都采用OPA進(jìn)行節(jié)點(diǎn)互連。如圖9所示，VLPLS在BDW集群上的強(qiáng)擴(kuò)展性略高于線性。當(dāng)運(yùn)行于更多的節(jié)點(diǎn)時，每個節(jié)點(diǎn)所需要的內(nèi)存變得更少，因此可以緩存的數(shù)據(jù)比例更多。這使得程序的數(shù)據(jù)局部性隨節(jié)點(diǎn)數(shù)增加而提高。

算例A中，模擬的網(wǎng)格空間共有1 200×1 200個網(wǎng)格，每個網(wǎng)格中有16個粒子。總共的內(nèi)存開銷是1 200×1 200×（網(wǎng)格數(shù)據(jù)大小+16×粒子數(shù)據(jù)大?。?1 735.83 MB。當(dāng)擴(kuò)展到32個雙路E5-2697v4節(jié)點(diǎn)時，可用的L3緩存大小是2 880 MB，這已經(jīng)大于程序的總內(nèi)存需求了。因此VLPL-S程序在BDW集群上表現(xiàn)出了極好的擴(kuò)展性。

在Intel Trace Analyzer and Collector的幫助下，本文測試了每個MPI操作所消耗的時間。程序中唯一的一個MPI集體通信操作MPI_Bcast占據(jù)了大部分MPI通信時間，尤其是在使用4個及更多KNL節(jié)點(diǎn)的時候。程序中MPI_Bcast僅用于同步一個32位狀態(tài)變量，在假設(shè)程序狀態(tài)沒有異常的情況下消除了MPI_Bcast操作。如圖9所示，程序移除MPI_Bcast后的擴(kuò)展性獲得了較大的提升。對MPI_Bcast操作在KNL平臺和CPU平臺進(jìn)行了詳細(xì)的性能測試，結(jié)果如圖10所示。當(dāng)使用兩個及更多節(jié)點(diǎn)時，MPI_Bcast在KNL平臺上的延遲要比其在CPU平臺上的延遲高出3倍以上。同時，在單節(jié)點(diǎn)上運(yùn)行多于64個進(jìn)程之后，MPI_Bcast的性能也極度下降。測試結(jié)果見圖11。這可能是超線程帶來的硬件共享所導(dǎo)致的。

Fig.9 Strong scalability of VLPL-S on BDW and KNL clusters圖9 VLPL-S在BDW和KNL平臺上的強(qiáng)擴(kuò)展性

Fig.10 MPI_Bcast performance on BDW and KNL clusters圖10 MPI_Bcast在BDW及KNL集群上的性能

（4）向量化

本文采用的臨時向量化方法導(dǎo)致了程序性能下降，尤其是在KNL平臺。進(jìn)行數(shù)據(jù)拼裝打包后，訪存時間與計(jì)算時間不再重疊。本文測試了程序的緩存缺失率后發(fā)現(xiàn)，應(yīng)用此向量化方法缺失率上升。之前預(yù)取優(yōu)化后，算例A在KNL平臺上取得了較高的加速比，但向量化之后，每次計(jì)算前需要讀取16個粒子的數(shù)據(jù)，弱化了之前預(yù)取帶來的訪存隱藏效果。由此原因，在KNL平臺上原本受益于預(yù)取優(yōu)化最多的算例A在臨時的向量化優(yōu)化之后性能下降也最多。

顯然，當(dāng)前版本的VLPL-S程序還離極致的性能很遠(yuǎn)。將鏈表數(shù)據(jù)結(jié)構(gòu)改寫為Struct of Array（SOA）是優(yōu)化過程中關(guān)鍵的一步。

Fig.11 MPI_Bcast performance on single-node KNL圖11 MPI_Bcast在單個KNL節(jié)點(diǎn)上的性能

6.4 總結(jié)

總體來說，對于主要用于性能測試的算例A而言，優(yōu)化后的VLPL-S程序在Xeon Phi 7210上的性能是其在雙路E5-2699v3上的1.77倍，也是其在雙路E5-2697v4上性能的1.53倍。優(yōu)化后的程序性能如圖12所示。圖中HSW表示雙路E5-2699v3；BDW表示雙路E5-2697v4；KNL表示Xeon Phi 7210 Flat Quadrant模式。64個進(jìn)程及每個進(jìn)程4個線程的運(yùn)行配置在KNL平臺上取得了最好的性能，18個進(jìn)程及每個進(jìn)程兩個線程的配置在CPU平臺取得了最好的成績。對于算例B而言，VLPL-S程序在Xeon Phi 7210上的性能要比其在雙路E5-2697v4上的低，這是由算例B的負(fù)載不均衡性造成的。負(fù)載不均衡性將兩個平臺的性能對比轉(zhuǎn)變?yōu)閱魏诵阅艿膶Ρ龋鳮NL的核心主頻較低，因此不占優(yōu)勢。

Fig.12 Performance comparision of VLPL-S on different platforms圖12 VLPL-S在不同平臺上的性能比較

盡管不同的優(yōu)化方法在不同平臺上取得了不同性能提升，這些常用的優(yōu)化手段在KNL及CPU平臺都是有效的。盡管MPI的性能在使用OPA互連的KNL平臺上要低一些，使用多線程優(yōu)化后VLPL-S在KNL集群上的擴(kuò)展性也是可以接受的。由此可見，對于VLPL-S來說多線程優(yōu)化在KNL平臺上是十分必要的。

7 結(jié)論及未來工作

綜上所述，對于VLPL-S來說常規(guī)的優(yōu)化方法在Xeon PhiTM和Xeon平臺上都是有效的。KNL平臺的MCDRAM提供了極高的訪存帶寬，從而極大地提高了VLPL-S程序的性能。MPI與OpenMP混合的運(yùn)行模式在單節(jié)點(diǎn)及多節(jié)點(diǎn)上實(shí)現(xiàn)了最佳性能。

鑒于目前所采用的臨時向量化方法對性能沒有提升，VLPL-S程序必須使用SOA替換鏈表數(shù)據(jù)結(jié)構(gòu)才能實(shí)現(xiàn)極致的性能。這部分工作目前正在進(jìn)行中，對SOA版本程序進(jìn)行性能評估及分析后會做更進(jìn)一步的討論。在SOA版本的基礎(chǔ)上，可以做更多針對MCDRAM及512位向量處理單元的優(yōu)化。

Intel MPI在KNL平臺上的性能低于其在CPU平臺上的性能，尤其是MPI集體操作的性能。其原因是目前MPI在KNL平臺上的實(shí)現(xiàn)還不成熟，之后會對此問題進(jìn)行更詳細(xì)的研究與討論。

致謝林新華特別致謝日本學(xué)術(shù)振興會JSPS的RONPAKU項(xiàng)目資助。陳民感謝Pukhov教授提供了VLPL程序的源代碼以及他對本文工作的支持。感謝Intel亞太研發(fā)中心高性能計(jì)算組里王哲及其他同事的優(yōu)化建議。

[1]Pukhov A.Three-dimensional electromagnetic relativistic particle-in-cell code VLPL(Virtual Laser Plasma Lab)[J].Journal of Plasma Physics,1999,61(3):425-433.

[2]Faure J,Glinec Y,Pukhov A,et al.A laser-plasma accelerator producing monoenergetic electron beams[J].Nature,2004,431(7008):541-544.

[3]Pukhov A,Meyer-Ter-Vehn J.Laser wake field acceleration:the highly non-linear broken-wave regime[J].Applied Physics B,2002,74(4/5):355-361.

[4]Pukhov A.Strong field interaction of laser radiation[J].Reports on Progress in Physics,2002,66(1):R1-R55.

[5]Chen Min,Luo Ji,Li Feiyu,et al.Tunable synchrotron-like radiation from centimeter scale plasma channels[J].Light:Science&Applications,2016,5(1):e16015.

[6]Bastrakov S,Meyerov I,Surmin I,et al.Particle-in-cell plasma simulation on CPUs,GPUs and Xeon Phi coprocessors[C]//LNCS 8488:Proceedings of the 29th International Supercomputing Conference,Leipzig,Jun 22-26,2014.Berlin,Heidelberg:Springer,2014:513-514.

[7]Surmin I A,Bastrakov S I,Efimenko E S,et al.Particle-incell laser-plasma simulation on Xeon Phi coprocessors[J].Computer Physics Communications,2016,202:204-210.

[8]Surmin I,Bastrakov S,Matveev Z,et al.Co-design of a particle-in-cell plasma simulation code for Intel Xeon Phi:a first look at Knights Landing[C]//LNCS 10049:Proceedings of the 2016 International Conference on Algorithms and Architectures for Parallel Processing,Granada,Dec 14-16,2016.Berlin,Heidelberg:Springer,2016:319-329.

[9]Nakashima H.Manycore challenge in particle-in-cell simulation:how to exploit 1 TFlops peak performance for simulation codes with irregular computation[J].Computers&Electrical Engineering,2015,46:81-94.

[10]Grigoryev Y N,Vshivkov V A,Fedoruk M P.Numerical"particle-in-cell"methods:theory and applications[M].Berlin:Walter de Gruyter,2002.

[11]Reinders J.Knights Landing,an overview for developers[EB/OL].[2016-06-27].https://software.intel.com/sites/default/files/managed/f3/87/Slides--Knights-Landing--An-Overview-for-Developers.pdf.

[12]Sodani A.Knights landing(KNL):2nd generation Intel?Xeon Phi processor[C]//Proceedings of the Hot Chips 27 Symposium,Cupertino,Aug 22-25,2015.Washington:IEEE Computer Society,2015:1-24.

[13]McCalpin J D.A survey of memory bandwidth and machine balance in current high performance computers[J].IEEE Technical Committee on Computer Architecture Newsletter,1995:19-25.

[14]Driesen K,H?lzle U.The direct cost of virtual function calls in C++[J].ACM Sigplan Notices,1996,31(10):306-323.

[15]Vanderwiel S P,Lilja D J.Data prefetch mechanisms[J].ACM Computing Surveys,2000,32(2):174-199.

[16]Barney B.OpenMP tutorial[EB/OL].[2016-07-29].https://computing.llnl.gov/tutorials/openMP/.