VLPL-S在Knights Landing上的優(yōu)化與性能評(píng)估*

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

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

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

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

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

1 引言

隨著高性能計(jì)算領(lǐng)域的快速發(fā)展，激光等離子體模擬程序可以在更大維度以及更高的精度下進(jìn)行?！疤摂M激光等離子體實(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)化了碰撞電離、場(chǎng)致電離以及兩體碰撞效應(yīng)和粒子跟蹤，并提供了對(duì)眾核平臺(tái)的支持。在代碼中，電子的離化是使用虛函數(shù)開發(fā)實(shí)現(xiàn)的，體現(xiàn)了面向?qū)ο缶幊陶Z言的特性；移動(dòng)窗口方法也應(yīng)用到了程序中以研究加速能量達(dá)GeV的長距離激光尾波場(chǎng)加速問題；通過加入以上的模擬方法，程序可以對(duì)Terahertz輻射、中性氣體與激光相互作用中的光脈沖及離化波前演化、激光尾波場(chǎng)加速及表面電子加速等物理過程進(jìn)行研究[5]。

在本文中，VLPL-S使用多種方法優(yōu)化，并在Intel?第二代Xeon PhiTM處理器（代號(hào)Knights Landing）上進(jìn)行性能評(píng)估。由于原程序的熱點(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章對(duì)優(yōu)化后的性能進(jìn)行評(píng)估；第7章給出結(jié)論以及后續(xù)優(yōu)化思路。

2 相關(guān)工作

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

此外，Surmin等人[7]開發(fā)了一套KNC平臺(tái)上的高性能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é)者還評(píng)估了PICADOR在KNL上的性能[8]，僅重新針對(duì)KNL平臺(tái)編譯代碼就獲得了相對(duì)于KNC上2.43倍的加速比。在優(yōu)化后，PIDADOR又獲得了1.89倍的加速比，并且實(shí)現(xiàn)了100 GFLOPS的雙精度峰值性能。

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

3 VLPL-S介紹

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

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

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

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

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

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

4 Intel Knights Landing

Knights Landing是Intel?Xeon PhiTM眾核處理器的第二代產(chǎn)品，它具有多個(gè)新的特性。首先KNL支持自啟動(dòng)，也就是說它可以作為一個(gè)獨(dú)立的處理器來使用而不僅僅是協(xié)處理器。其次，藉由其片上整合的MCDRAM，KNL提供了400+GB/s的訪存帶寬，同時(shí)還包含傳統(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)帶寬的對(duì)比可見表1。

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

KNL的芯片由36個(gè)Tile組成。這些Tile由二維片上網(wǎng)格所連接。每個(gè)Tile中有兩個(gè)核以及1 MB的共享L2緩存。每核中有兩個(gè)向量處理單元。詳細(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)

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

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

5.2 訪存優(yōu)化

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

預(yù)取是一個(gè)用于減少訪存延遲影響的常用優(yōu)化方法。硬件預(yù)取與軟件預(yù)取共同工作以給處理器的計(jì)算單元提供數(shù)據(jù)[15]。然而，目前的硬件預(yù)取很難對(duì)遍歷鏈表時(shí)的隨機(jī)訪存進(jìn)行預(yù)測(cè)，即便使用Intel ICC編譯器也不能對(duì)其做自動(dòng)軟件預(yù)取優(yōu)化。為了取得更高的訪存性能，使用Intrinsic庫將鏈表中下一個(gè)要處理的粒子數(shù)據(jù)預(yù)取到L1緩存。由于鏈表的訪存特性，在讀取到當(dāng)前節(jié)點(diǎn)的數(shù)據(jù)之前并不知道下一個(gè)節(jié)點(diǎn)的數(shù)據(jù)位置，因此單純的預(yù)取也不能完全隱藏鏈表訪存延遲。為此，在每個(gè)網(wǎng)格中維護(hù)了一個(gè)前向查找表，以便在計(jì)算每個(gè)網(wǎng)格時(shí)預(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算法中的粒子推動(dòng)函數(shù)，它占用了程序執(zhí)行時(shí)間的80%以上。其中由粒子運(yùn)動(dòng)產(chǎn)生的電流數(shù)據(jù)需要更新到網(wǎng)格當(dāng)中，網(wǎng)格中的不同粒子可能會(huì)同時(shí)更新同一個(gè)位置的數(shù)據(jù)。因此，在多線程并行化的同時(shí)需要處理數(shù)據(jù)沖突的問題。

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

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

（1）原子操作

OpenMP庫提供了原子操作語句，其指定的內(nèi)存地址同時(shí)只能由一個(gè)線程操作，這就可以避免數(shù)據(jù)寫沖突。具體來說OpenMP提供了一個(gè)最小化的關(guān)鍵區(qū)塊以實(shí)現(xiàn)原子操作[16]。從PIC算法上來說，每個(gè)粒子只需要對(duì)其周圍的網(wǎng)格做出貢獻(xiàn)，因此只需在每個(gè)線程的邊界區(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）被動(dòng)貢獻(xiàn)方法

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

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

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

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

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

5.4 向量化

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

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

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

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

6 結(jié)果討論

6.1 算例及測(cè)試平臺(tái)說明

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

Table 2 Details about workloads表2 測(cè)試所用算例說明

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

Table 3 Machine configuration表3 測(cè)試平臺(tái)配置說明

6.2 在KNL上對(duì)原版VLPL-S代碼的性能評(píng)估

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

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

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

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平臺(tái)上的加速比

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

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

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

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

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

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

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

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

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

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

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

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

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

（4）向量化

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

6.4 總結(jié)

總體來說，對(duì)于主要用于性能測(cè)試的算例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個(gè)進(jìn)程及每個(gè)進(jìn)程4個(gè)線程的運(yùn)行配置在KNL平臺(tái)上取得了最好的性能，18個(gè)進(jìn)程及每個(gè)進(jìn)程兩個(gè)線程的配置在CPU平臺(tái)取得了最好的成績。對(duì)于算例B而言，VLPL-S程序在Xeon Phi 7210上的性能要比其在雙路E5-2697v4上的低，這是由算例B的負(fù)載不均衡性造成的。負(fù)載不均衡性將兩個(gè)平臺(tái)的性能對(duì)比轉(zhuǎn)變?yōu)閱魏诵阅艿膶?duì)比，而KNL的核心主頻較低，因此不占優(yōu)勢(shì)。

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

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

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

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

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

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

致謝林新華特別致謝日本學(xué)術(shù)振興會(huì)JSPS的RONPAKU項(xiàng)目資助。陳民感謝Pukhov教授提供了VLPL程序的源代碼以及他對(duì)本文工作的支持。感謝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/.