“神威·太湖之光”上Tend_lin應(yīng)用的并行優(yōu)化研究

姜尚志，唐生林，高希然，花 嶸，陳 莉，劉 穎

(1.山東科技大學(xué)計算機(jī)科學(xué)與工程學(xué)院，山東 青島 266590;2.中國科學(xué)院計算技術(shù)研究所計算機(jī)體系結(jié)構(gòu)國家重點實驗室,北京 100190)

1 引言

隨著全球變暖的不斷加劇，氣候問題成為關(guān)系到人類社會生存與發(fā)展的重大戰(zhàn)略議題，受到各國政府、機(jī)構(gòu)和學(xué)術(shù)界的高度重視。西方發(fā)達(dá)國家紛紛開發(fā)具有更高精度的氣候模擬程序，如CAM5(Community Atmosphere Model version 5)[1]、ECHAM5(European Center HAmburg Model version 5)[2]、MRI-AGCM(Meteorological Research Institute Atmospheric General Circulation Models)[3]、BCC_AGCM(Beijing Climate Center Atmospheric General Circulation Model)[4]。我國也聯(lián)合多家研究機(jī)構(gòu)研制了中國科學(xué)院地球系統(tǒng)模式CAS-ESM(Earth System Model of Chinese Academy of Sciences)[5]。本文關(guān)心的IAP AGCM[6]是中國科學(xué)院大氣物理研究所IAP(Institute of Atmospheric Physics)研發(fā)的大氣環(huán)流模式AGCM(Atmospheric General Circulation Model)，是CAS-ESM的5個主要分量模式之一[7]，也是整個CAS-ESM中運行時間占比最大的分量模式。IAP AGCM自20世紀(jì)80年代發(fā)布第1代[8]以來，已歷經(jīng)4代，氣候模擬能力得到了顯著提升[6,9,10]，并亟需移植到新一代超級計算機(jī)上以進(jìn)一步提高其模擬速度。

E級超級計算系統(tǒng)將于2020年出現(xiàn)，各國都將全球氣候變化的模擬列為重點應(yīng)用領(lǐng)域。近年來，由于訪存墻、指令墻、并行墻和功耗墻等因素的限制，異構(gòu)眾核已成為超級計算機(jī)的主流體系結(jié)構(gòu)構(gòu)架。如何把復(fù)雜的氣候模擬程序移植到新一代的超級計算機(jī)上并取得更好的并行性能具有重要的現(xiàn)實意義?！吧裢ぬ狻笔俏覈灾餮邪l(fā)的超級計算機(jī)，也是全世界第1臺峰值性能超過100 PFLOPS的超級計算機(jī)，曾4次蟬聯(lián)TOP500榜首，其計算結(jié)點采用國產(chǎn)異構(gòu)眾核處理器申威26010[11]，是典型的異構(gòu)眾核加速結(jié)構(gòu)。研究氣候模式應(yīng)用在“神威·太湖之光”上的并行優(yōu)化，對于將來應(yīng)用到國產(chǎn)E級計算機(jī)的并行優(yōu)化具有借鑒意義。

在大規(guī)模異構(gòu)眾核的并行計算系統(tǒng)上，程序員主要使用MPI+X的混合編程模式進(jìn)行并行編程，其中MPI負(fù)責(zé)進(jìn)程級的并行和通信，而模型X對應(yīng)進(jìn)程內(nèi)異構(gòu)加速的并行編程模型，X的選擇有OpenMP[12]、OpenACC[13]、OpenCL和CUDA等?；谥茖?dǎo)的編程模型，比如OpenACC和OpenMP，因為把異構(gòu)加速代碼的生成交給編譯器，降低了程序員負(fù)擔(dān)而廣泛用于MPI應(yīng)用到大規(guī)模異構(gòu)眾核并行計算系統(tǒng)的遷移。但是，現(xiàn)有主流的異構(gòu)眾核編程模型對用戶標(biāo)記的待加速循環(huán)，主要以fork-join的鎖步推進(jìn)模式進(jìn)行加速，而難以發(fā)掘不同循環(huán)間的并行性。

近年來，數(shù)據(jù)驅(qū)動的任務(wù)編程模型受到國際學(xué)術(shù)界的廣泛關(guān)注，它因為能自然地表達(dá)并行性而有利于應(yīng)對當(dāng)代多核、眾核系統(tǒng)上的可編程性挑戰(zhàn)和性能挑戰(zhàn)。數(shù)據(jù)驅(qū)動的任務(wù)編程模型將應(yīng)用表示為以任務(wù)為節(jié)點、任務(wù)間數(shù)據(jù)依賴關(guān)系為有向依賴邊的有向無環(huán)圖DAG(Directed Acyclic Graph)，用戶只需提供每個任務(wù)的數(shù)據(jù)訪問信息，而依賴關(guān)系的構(gòu)建和任務(wù)調(diào)度執(zhí)行則交由底層運行時庫完成。這類編程模型有OpenMP (4.0以上)、StarPU[23]、OmpSs[24]、Cilk和AceMesh[25]等。

地球系統(tǒng)模式到大規(guī)模異構(gòu)眾核平臺的移植和優(yōu)化得到了廣泛的研究。Yang等人[14]利用神威眾核線程加速庫Athread實現(xiàn)了MPI+Athread的高效并行；Fu等人[15,16]在“神威·太湖之光”上借助OpenACC、自研的循環(huán)變換工具、內(nèi)存足跡分析優(yōu)化工具對CAM的動力框架CAS-SE進(jìn)行了代碼重構(gòu)和并行，加速效果顯著，部分熱點取得了22倍的加速；李鑌洋等人[17]和吳琦等人[18]分別實現(xiàn)了海冰模式和區(qū)域海洋模式在“神威·太湖之光”上的并行化。肖灑等人[4]和魏敏等人[19]先后對國家氣候中心大氣環(huán)流模式BCC_AGCM的部分核心段在“神威·太湖之光”上使用OpenACC進(jìn)行異構(gòu)眾核加速；傅游等人[20]提取IAP AGCM-4的動力框架熱點模塊Tend_lin構(gòu)建了獨立的測試程序，在“神威·太湖之光”上使用OpenACC實現(xiàn)了最高1.4°分辨率下的異構(gòu)眾核加速，單核組多線程獲得6倍以上加速。

在并行算法優(yōu)化方面，Xiao等人[21]針對IAP AGCM的動力框架，為計算流推導(dǎo)出新的基本操作，通過引入近似項等方式，提出了通信避免算法，大規(guī)模降低了集合通信和點到點通信的頻率，同時也最大化了通信隱藏的機(jī)會，最終把整個模式的執(zhí)行時間降低了54%。針對二維經(jīng)緯網(wǎng)格的分解并行擴(kuò)展性低的問題，Wu等人[22]為IAP AGCM的動力框架引入三維區(qū)域分解及自適應(yīng)高斯濾波方案，使得0.5°分辨率的模擬在“天河二號”上可擴(kuò)展到32 768個CPU核。本文的工作集中在進(jìn)程內(nèi)的多線程并行，以及自動地實現(xiàn)通信-計算重疊，與以上2個進(jìn)程間的并行算法優(yōu)化是互補(bǔ)的。

本文面向“神威·太湖之光”并行計算平臺，選擇IAP AGCM-4中動力框架過程的適應(yīng)過程Tend_lin，分別利用OpenACC和中國科學(xué)院計算技術(shù)研究所研發(fā)的數(shù)據(jù)驅(qū)動的任務(wù)并行編程接口AceMesh對其進(jìn)行并行優(yōu)化，獲得了比較好的性能加速。本文貢獻(xiàn)包括：

(1)針對高分辨率模擬的并行特點，利用OpenACC對Tend_lin進(jìn)行眾核并行的優(yōu)化，指出了OpenACC的局限性。

(2)利用數(shù)據(jù)驅(qū)動的任務(wù)并行編程接口AceMesh對Tend_lin進(jìn)行異步并行的優(yōu)化。具體包括AceMesh并行區(qū)中全局?jǐn)?shù)據(jù)分塊的選擇，計算循環(huán)和通信代碼的任務(wù)并行方法，并針對Tend_lin應(yīng)用中計算/通信比低、進(jìn)程內(nèi)并行度不足的特點，討論了放松通信資源共享的通信去串行化，以及單層任務(wù)圖和嵌套任務(wù)圖下的任務(wù)映射優(yōu)化。

(3)在“神威·太湖之光”平臺上對Tend_lin的并行優(yōu)化進(jìn)行了評估和性能分析。測試表明,相比OpenACC并行版本，AceMesh在16～1 024進(jìn)程的不同并行配置下均得到了平均2倍以上的性能提升。本文還詳細(xì)分析了性能收益的來源。

2 Tend_lin過程和Tend_lin程序

大氣環(huán)流模式(AGCM)是由最初的數(shù)值天氣預(yù)報模式演變而來的，是地球系統(tǒng)模式中最為復(fù)雜的分量模式，計算量占地球系統(tǒng)模式中全部計算量的50%以上。IAP AGCM-4[6]是中國科學(xué)院大氣物理研究所研發(fā)的第4代大氣環(huán)流模式，采用均勻經(jīng)緯格點的有限差分?jǐn)?shù)值方法，目前最高支持到0.25°分辨率，主要包括物理過程和動態(tài)框架2個過程。其中，物理過程在IAP AGCM-4中的時間占比隨著CPU核數(shù)的增加而降低，在2進(jìn)程時時間占比23.8%，在1 024進(jìn)程時時間占比僅12.0%。自主設(shè)計的動力框架過程是IAP AGCM-4模式中最耗時的模塊之一，是用來求解關(guān)于時間的偏微分方程組。

動力框架的計算在空間上是三維的，每一格點的計算都可能與周圍的很多格點有關(guān)。動力框架方程組[6]如下所示：

(1)

其中，各參數(shù)的含義請見參考文獻(xiàn)[6]。

方程組中涉及平流過程(Tend_adv)和適應(yīng)過程(Tend_lin)2個主要的差分計算，其中適應(yīng)過程的調(diào)用頻率很高。

本文關(guān)注適應(yīng)過程(Tend_lin)，該過程主要是對經(jīng)緯圈和垂直層上的溫度、風(fēng)向模擬變量等進(jìn)行數(shù)值差分計算，得到U(緯向)、V(經(jīng)向)、T(溫度)、P(氣壓)的傾向值。動力框架采用緯度、高度的二維區(qū)域分解進(jìn)行MPI多進(jìn)程并行。為了保證數(shù)值計算的正確性和有效性，在其中插入濾波和平滑過程(Filt、Smoother)，剔除計算中偏差大的數(shù)值。

Tend_lin過程要計算4個模擬變量，這些變量都會經(jīng)歷2個階段的計算，分別是差分方程的計算和模擬結(jié)果的濾波平滑。后文把2個階段分別簡稱為stencil階段和濾波階段。這4個模擬變量的計算過程中還存在一些共用的變量和通信，整個適應(yīng)過程可以用圖1的模塊依賴關(guān)系圖進(jìn)行抽象。

從圖1所示的依賴關(guān)系圖中可以看出，DU/DV2個變量總是一起計算，而DT和DPsa在計算過程中也存在一定的關(guān)聯(lián)性，比如它們的極點區(qū)域是在一起計算的，極點區(qū)域的DT/DPsa的計算涉及一些變量在高度方向的歸約求和。而TT1則要對DIV和TT2個變量在每個高度通信域上進(jìn)行mpi_allgatherv收集，涉及的通信量很大。在差分計算階段有很多緯度方向的近鄰?fù)ㄐ拧ｊP(guān)于平滑處理，DU/DV采用的是一維緯向平滑方法(SHAP1)，計算量較小，也不需要進(jìn)程間的通信。而DT/DPsa采用二維平滑方法(SHAP2)，需要緯度方向的近鄰?fù)ㄐ拧?/p>

Figure 1 Diagram of dependencies between modules computed with four variables圖1 4個變量計算模塊之間的依賴關(guān)系圖

3 數(shù)據(jù)驅(qū)動的AceMesh并行優(yōu)化

3.1 AceMesh任務(wù)并行編程接口

AceMesh是基于制導(dǎo)的、數(shù)據(jù)驅(qū)動的任務(wù)并行編程接口，其編譯器是一個基于rose框架的、源源變換的編譯器，變換后的源代碼將再經(jīng)過目標(biāo)平臺上的本地編譯器編譯、并鏈接AceMesh運行時庫得到可執(zhí)行代碼。AceMesh編譯系統(tǒng)目前支持X86的多核平臺和國產(chǎn)神威平臺。

AceMesh支持增量并行化，表1列出了它的主要制導(dǎo)。begin/end制導(dǎo)用于指定一個待并行化的程序代碼區(qū)域，這個單入單出的代碼塊也稱為任務(wù)并行區(qū)，因為并行化后的程序以有向無環(huán)圖的任務(wù)圖方式執(zhí)行，這個并行區(qū)也稱為DAG并行區(qū)。AceMesh提供了2種類型的任務(wù)制導(dǎo)(task制導(dǎo)和do制導(dǎo))，讓用戶描述相關(guān)任務(wù)的數(shù)據(jù)依賴信息。DAG區(qū)內(nèi)的源代碼分為任務(wù)代碼和任務(wù)構(gòu)建代碼2類。任務(wù)代碼是那些出現(xiàn)在任務(wù)制導(dǎo)動態(tài)范圍內(nèi)的代碼，它們將被外聯(lián)為任務(wù)函數(shù)，這些任務(wù)的執(zhí)行被推遲，推遲到什么時刻在哪個計算核心上執(zhí)行取決于運行時庫。任務(wù)構(gòu)建代碼由主線程執(zhí)行，負(fù)責(zé)構(gòu)建任務(wù)并行區(qū)內(nèi)的所有任務(wù)，包括構(gòu)建其任務(wù)對象和參數(shù)、任務(wù)之間的依賴關(guān)系、實現(xiàn)任務(wù)之間的同步等。由于并發(fā)任務(wù)相對于構(gòu)圖階段而言是延遲執(zhí)行的，用戶應(yīng)該保證不存在從任務(wù)代碼到任務(wù)構(gòu)建階段的數(shù)據(jù)依賴關(guān)系，否則需要追加taskwait制導(dǎo)實現(xiàn)同步。arrayTile以數(shù)據(jù)分塊的方式提供了DAG并行區(qū)中的并行模板。AceMesh允許出現(xiàn)嵌套的DAG并行區(qū)，也就是在一個并發(fā)任務(wù)內(nèi)部，還可以定義tasks，使用do制導(dǎo)，聲明arrayTile、taskwiat制導(dǎo)。task/do制導(dǎo)的map子句的參數(shù)可以是master或acc，分別表示并發(fā)任務(wù)將在主線程上執(zhí)行，或者在加速器中執(zhí)行。

Table 1 Summary of AceMesh’s directives表1 AceMesh主要制導(dǎo)

對于并行循環(huán)，我們使用的do制導(dǎo)如下所示：

!$acemesh do [map(master|acc)] tasktile(loop_tile_list) in() out() inout() nested [if()] [nested] [private()] [firstprivate()] [lastprivate()] [reduction()] [label〈string〉]

tasktile子句描述了如何將相關(guān)循環(huán)的迭代空間劃分為塊，并將每個塊封裝到一個數(shù)據(jù)流任務(wù)中。do/task制導(dǎo)的nested子句指示該任務(wù)內(nèi)部存在嵌套的任務(wù)依賴關(guān)系圖。in/out/inout子句描述每個任務(wù)的數(shù)據(jù)訪問區(qū)域，其列表項可以使用變量和數(shù)組區(qū)域。數(shù)組區(qū)域的每個下標(biāo)可以是一個三元組' [lo_expr]:[up_expr]:[step] '，其中l(wèi)o_expr和up_expr是tasktile子句中循環(huán)索引量的仿射表達(dá)式，我們還使用一個特殊的*表示整個維展區(qū)間。實際上，很多時候數(shù)組訪問區(qū)域可以簡化為一個數(shù)組名，因為具體的數(shù)組訪問區(qū)域通?？梢杂葾ceMesh編譯器自動推導(dǎo)出來。label是圖形化工具展示用的，提供task_graph和loop_graph上的循環(huán)標(biāo)記。

對于do制導(dǎo)，AceMesh編譯器會根據(jù)tasktile子句提供的循環(huán)分塊信息對并行循環(huán)進(jìn)行分段變換，每個循環(huán)分塊對應(yīng)一個并發(fā)任務(wù)，并按照map子句的指示，生成主核任務(wù)代碼或者從核任務(wù)代碼。

3.2 全局?jǐn)?shù)據(jù)分塊和計算循環(huán)的任務(wù)并行方法

用戶在對任務(wù)并行區(qū)進(jìn)行并行的時候，要先對并行區(qū)進(jìn)行全局分析，確定主流的數(shù)據(jù)分塊模式，盡量讓每個并行循環(huán)套的并行策略與數(shù)據(jù)分塊模板對齊。然后根據(jù)數(shù)據(jù)分塊的模式，具體決定每個并行循環(huán)的任務(wù)分塊方法。

(1)stencil階段的并行化。

stencil階段(也就是差分計算階段)的多數(shù)計算模式是對緯度-高度-經(jīng)度組成的三維網(wǎng)格空間進(jìn)行掃描計算，大多數(shù)循環(huán)套在3個軸向都是可并行的，只有少數(shù)幾個循環(huán)套的高度維不能并行而是要做復(fù)雜的規(guī)約計算。對整個階段而言，本文選擇對緯度空間進(jìn)行一維數(shù)據(jù)分塊，采用default子句為大多數(shù)數(shù)組聲明其分塊方式。然后對每個循環(huán)套進(jìn)行并行時，讓循環(huán)的任務(wù)分塊盡量與數(shù)據(jù)分塊對齊。算法1是其中一個典型循環(huán)套的并行化算法。DT數(shù)組的allocate語句是從其他文件中提取出來的，用來展示數(shù)組維展。第1～2行給出的是數(shù)據(jù)分塊制導(dǎo)，其中緯度分塊尺寸是tilej。第3行給出的是循環(huán)的任務(wù)并行方法，因為DT和ST的數(shù)組下標(biāo)表達(dá)式的仿射系數(shù)都是1，所以循環(huán)分塊的尺寸與對應(yīng)數(shù)組維的分塊尺寸相同，這樣每個任務(wù)正好為一個數(shù)據(jù)分塊進(jìn)行定值，而in/out子句的列表項，本文采用簡化寫法，只給出數(shù)組名，由編譯器推導(dǎo)任務(wù)訪問的數(shù)據(jù)區(qū)域。

算法1數(shù)據(jù)分塊制導(dǎo)和stencil階段第14循環(huán)的任務(wù)并行算法

allocate (DT(NX,beglev:endlev,beglatdynex:endlatdynex) )

… …

1 !$acemesh arrayTile dimtile(*,tilej) default

2 !$acemesh arrayTile dimtile(*,*,tilej) default

3 !$acemesh do tasktile(j:tilej) inout(DT) in(ST)

4doJ=beglatdyn,endlatdyn

5doK=beglev,endlev

6doI= 1,NX

7DT(I,K,J) =DT(I,K,J)+ST(I,K,J)

8enddo

9enddo

10enddo

11 …

(2)濾波階段的并行化。

本階段調(diào)用一個FILT2D函數(shù)對每一個經(jīng)緯平面進(jìn)行濾波和平滑處理。算法2展示的是DT數(shù)組的濾波過程的任務(wù)并行算法，對于每個高度，先把待濾波的經(jīng)緯面拷貝到一個臨時數(shù)組WW，濾波處理后再拷貝回原始數(shù)組。由于這里沒有使用nested制導(dǎo)，我們稱這種任務(wù)并行方法為單層任務(wù)圖。在FILT2D函數(shù)內(nèi)部，有多級函數(shù)調(diào)用，其中的多數(shù)循環(huán)只有緯度方向可以并行分塊。算法3給出的是其中一個典型循環(huán)——低緯度濾波循環(huán)的任務(wù)并行情況。

算法2DT變量濾波的任務(wù)并行算法

1dok=beglev,endlev

2 $acemesh do tasktile (j:tilej) in(DT) out(WW2)

3doj=beglatdyn,endlatdyn

4WW(:,j)=DT(:,K,j)

5enddo

6 CallFILT2D(WW,0,1,IBCFFT)

7 !$acemesh do tasktile(j:tilej) out(DT) in(WW2)

8doj=beglatdyn,endlatdyn

9DT(:,K,j)=WW(:,j)

10enddo

11enddo

算法3濾波階段一個循環(huán)的AceMesh并行算法

1 real(r8),intent(inout) ::CH(NX,NY)

2 !$acemesh arrayTile dimtile(*,tilej) default

3 !$acemesh do tasktile(j:tilej) inout(CH)

4doJ= max(JBL,beglatdyn),min(JEL,endlatdyn)

5X0=ZERO

6doI= 1,IM,2

7X0 = (CH(I,J)-CH(I+1,J))+X0

8enddo

9X0=X0/FIM

10doI= 1,IM,2

11II=I+1

12CH(I,J) =CH(I,J)-X0

13CH(II,J)=CH(II,J)+X0

14enddo

15 callperiodp(CH(1,J))

16enddo

3.3 利用嵌套任務(wù)圖優(yōu)化計算任務(wù)的從核映射

本小節(jié)主要優(yōu)化濾波階段的并行性。濾波階段的代碼隱含了變量之間、高度之間、緯度之間3個層次的并行性，AceMesh編程接口可以自然地發(fā)掘這3個層次的并行性。

單層任務(wù)圖的任務(wù)映射缺陷如算法2所示，每個高度上都會進(jìn)行濾波計算，由于AceMesh根據(jù)任務(wù)的數(shù)據(jù)訪問區(qū)域決定親和性，而不同高度上的濾波都是在局部數(shù)組CH上進(jìn)行，在AceMesh看來它們的緯度區(qū)間是重疊的，這就導(dǎo)致這些并行循環(huán)上的并行任務(wù)總是映射到同樣的從核子集，而沒有均勻分散到從核陣列，于是丟失了第1層和第2層的濾波并行性。同時，AceMesh編譯器的運行時系統(tǒng)還沒有在“神威·太湖之光”上支持任務(wù)竊取，從而濾波階段的并行度沒有得到充分的發(fā)揮。

利用嵌套任務(wù)圖優(yōu)化濾波階段：算法4利用AceMesh的嵌套任務(wù)圖(后面簡稱嵌套圖)對濾波階段進(jìn)行并行化。相比算法2的單層任務(wù)圖(后面簡稱單層圖)風(fēng)格的代碼，這里新增加了第1行的數(shù)據(jù)分塊——對DT數(shù)組進(jìn)行緯度為tilej、高度為1的二維分塊，第2行的do制導(dǎo)對k循環(huán)進(jìn)行任務(wù)分塊，分塊尺寸是1——與數(shù)據(jù)分塊的尺寸相同。需要指出的是k循環(huán)內(nèi)的并行化制導(dǎo)以及FILT2D函數(shù)內(nèi)的制導(dǎo)都不需要修改(與算法3的版本一致)。第2行的并行循環(huán)將對應(yīng)很多并發(fā)任務(wù)，這里每個任務(wù)本身又是一個任務(wù)圖，從而形成嵌套的任務(wù)圖。AceMesh編譯系統(tǒng)對第2行進(jìn)行處理時，根據(jù)每個任務(wù)的DT訪問區(qū)域?qū)γ總€任務(wù)分配對應(yīng)的從核集合。編譯系統(tǒng)遇到算法2的二級任務(wù)時，將在一級任務(wù)所確定的從核集合中進(jìn)行從核映射。采用嵌套任務(wù)圖的并行方法后，濾波階段3個層次的并行性得到充分挖掘。

算法4用嵌套任務(wù)圖并行化濾波階段

1 !$acemesh arrayTile dimtile(*,1,tilej) dim(NX,beglev:endlev,beglatdyn:endlatdyn) arrvar(DT)

2 !$acemesh do tasktile(k:1) inout(DT(*,k,*)) private(WW) nested

3dok=beglev,endlev

!original loop body

4enddo

3.4 用task制導(dǎo)實現(xiàn)通信和計算的動態(tài)重疊

原程序?qū)c到點通信實現(xiàn)了一定程度的通信-計算的重疊優(yōu)化，但是沒有對集合通信做通信隱藏的優(yōu)化(實際上靜態(tài)的代碼調(diào)度很難隱藏這些優(yōu)化)。stencil階段對點到點通信的發(fā)起盡量前提，而把通信等待延遲到使用時。而filter階段的通信調(diào)度更加局限，沒有與計算進(jìn)行重疊，而僅僅是對每一對通信的收、發(fā)進(jìn)行了調(diào)度。

AceMesh版本對所有點到點通信和集合通信都進(jìn)行了任務(wù)并行的優(yōu)化。算法5給出的是通信發(fā)起操作的任務(wù)并行制導(dǎo)，因為消息接收是對消息緩沖的寫操作，而同時建立一個通信句柄，也就是通信句柄也是寫操作。算法6給出的是wait操作的任務(wù)并行制導(dǎo)，根據(jù)MPI規(guī)范的語義通信句柄是inout參數(shù)，所以該句柄的數(shù)據(jù)流屬性是既讀又寫。

對于MPI集合通信，AceMesh任務(wù)并行方法與此類似，不再贅述。需要指出的是，AceMesh自動把連續(xù)出現(xiàn)的通信發(fā)起代碼和通信等待代碼在調(diào)度時分開，盡可能早地提前發(fā)起通信操作，而對通信等待任務(wù)進(jìn)行不斷輪詢，避免線程在mpi_wait上等待而浪費處理器的時間，也避免了由此帶來的MPI死鎖。

算法5非阻塞通信的post操作的AceMesh并行算法

1if(myid_y.ne.0)then

2 !$acemesh task map (master) private(ierr) out(rbuf(*))out(rreq(2))

3rreq(2)=MPI_REQUEST_NULL

4 call mpi_irecv(rbuf,leng,mpir8,myid_y-1,10+myid_z,comm_y,rreq(2),ierr)

5 !$acemesh task end

6endif

算法6非阻塞通信的wait操作的AceMesh并行算法

1 !$acemesh task map (master) inout(req(*))

2 call mpi_waitall(2,req,status,ierr)

3 !$acemesh task end

3.5 利用tag私有化降低通信之間的串行化依賴

通信信封和通信保序：MPI為了實現(xiàn)通信行為的確定性，用通信信封對點到點通信進(jìn)行消息配對。通信信封是個四元組〈源進(jìn)程、目的進(jìn)程、通信標(biāo)簽、通信域〉。信封相同的消息在配對的時候要保持原來的發(fā)起順序，先發(fā)起的先配對，也就是說P1進(jìn)程發(fā)給P2進(jìn)程的信封相同的2個消息會按照確定的順序進(jìn)行，在P2進(jìn)程上接收順序與P1上的發(fā)送順序一致。對于非阻塞通信，序關(guān)系作用在通信發(fā)起操作(比如mpi_irecv)之間。

Tend_lin應(yīng)用中有4個變量需要進(jìn)行濾波處理，其中DT和DPsa2個變量的濾波過程中需要二維平滑處理，從而引入緯度方向的近鄰?fù)ㄐ?。FILT2D函數(shù)內(nèi)這些通信均采用了相同的通信標(biāo)簽和通信域，從而導(dǎo)致不同高度、不同變量的濾波通信之間必須保序，也就是所有向左的通信必須串行化發(fā)起，所有向右的通信也必須串行化發(fā)起，這就降低了應(yīng)用中的任務(wù)并行性。

通過使用更多的通信標(biāo)簽(tag)進(jìn)行去串行化：可以用2種方法來放松通信信封帶來的序關(guān)系：(1)復(fù)制出更多的通信域；(2)使用個性化的通信標(biāo)簽。因為通信域的復(fù)制會帶來很大的開銷，本文選擇對通信標(biāo)簽進(jìn)行重構(gòu)。對通信標(biāo)簽進(jìn)行基于計數(shù)器的資源池管理，提供一定數(shù)量的、不同的通信標(biāo)簽，降低通信發(fā)起操作之間的串行化趨勢，提高通信推進(jìn)的自由性。后文把這個優(yōu)化稱為tag私有，而原來的版本則稱為tag共享。目前最終測試中每個軸向通信的標(biāo)簽總數(shù)是20。

3.6 從核任務(wù)代碼的數(shù)據(jù)傳輸優(yōu)化

從核任務(wù)代碼的質(zhì)量在很大程度上影響到AceMesh程序的性能。AceMesh編譯器所生成的從核任務(wù)代碼已經(jīng)對變量進(jìn)行了LDM(Local Data Memory)局部存儲優(yōu)化。但是，在跨語言對比的時候，仍然發(fā)現(xiàn)了AceMesh編譯器從核代碼優(yōu)化的不足，本文對此進(jìn)行了手工補(bǔ)充。

非整除塊的reshape優(yōu)化：AceMesh的tasktile與OpenACC的循環(huán)映射的分塊語義不一樣，并行循環(huán)的第1個分塊和最后1個分塊都可能出現(xiàn)非完整分塊的情況。若并行分塊循環(huán)對應(yīng)到數(shù)組的非最高維，則非完整的循環(huán)分塊導(dǎo)致LDM數(shù)組中的目標(biāo)區(qū)域在地址上不連續(xù)，從而DMA數(shù)據(jù)傳輸被迫分成多次進(jìn)行，將會明顯降低傳輸性能。由于AceMesh編譯器采用了源源變換的方法，用戶可以查看變換后的從核任務(wù)代碼，并進(jìn)行進(jìn)一步的手工優(yōu)化。對于主存源數(shù)組區(qū)域的內(nèi)存地址連續(xù)的情況，本文采用一次DMA通信把整塊數(shù)據(jù)傳遞到LDM，然后在片上進(jìn)行數(shù)組轉(zhuǎn)置。類似地，這種非完整分塊的數(shù)據(jù)寫回也采用“先reshape-后DMA傳輸”的方式。

4 性能評估

實驗平臺：“神威·太湖之光”計算機(jī)系統(tǒng)的每個結(jié)點采用 申威26010異構(gòu)眾核處理器，該處理器集成4個運算核組，每個核組包含1個運算控制核心(主核)和64個運算核心(從核)，從核以8×8的Mesh結(jié)構(gòu)組成運算核心陣列(從核陣列)。每個從核擁有64 KB大小、SPM形式的高速局部存儲空間。

神威平臺的結(jié)點內(nèi)提供athread和OpenACC 2種異構(gòu)編程模型。其中OpenACC是一個面向多核、異構(gòu)眾核平臺的并行編程接口，支持C、C++、Fortran編程語言，支持多種硬件架構(gòu)。神威平臺的OpenACC*基于OpenACC2.0標(biāo)準(zhǔn)，針對申威26010處理器結(jié)構(gòu)特點進(jìn)行了適當(dāng)?shù)木喓蛿U(kuò)充。一些有特色的擴(kuò)展包括：(1)針對帶跨步的數(shù)組傳輸，提供pack/swap等子句，利用數(shù)據(jù)打包、數(shù)組轉(zhuǎn)置優(yōu)化DMA數(shù)據(jù)傳輸效率；(2)為data制導(dǎo)擴(kuò)展index子句，允許用戶設(shè)置數(shù)據(jù)打包和轉(zhuǎn)置的放置點，把數(shù)據(jù)打包和轉(zhuǎn)置動作進(jìn)行外提。

本節(jié)在“神威·太湖之光”平臺上分別對OpenACC和AceMesh實現(xiàn)的2個Tend_lin并行版進(jìn)行性能評估和分析。本文采用2種分辨率，其三維網(wǎng)格空間分別是720×361×30的0.5°×0.5°(50 km)和1440×721×30的0.25°×0.25°(25 km)。

對于OpenACC版本，本文檢查了源源變換后的每個從核代碼，確保變量都進(jìn)行了LDM優(yōu)化。

參數(shù)調(diào)優(yōu)：對于OpenACC和AceMesh 2個版本，本文對每個并行循環(huán)設(shè)置了不同的分塊參數(shù)，以方便調(diào)優(yōu)。MPI二維區(qū)域分解中，不同進(jìn)程布局對性能影響較大，本文對進(jìn)程數(shù)進(jìn)行因式分解的窮舉測試；對以上輸入采用調(diào)優(yōu)參數(shù)配置文件的形式完成批量測試，測試完成后自動在測試結(jié)果中選優(yōu)。

4.1 Tend_lin的OpenACC優(yōu)化

OpenACC并行化基于傅游等人[20]的前期工作，此前已經(jīng)進(jìn)行了循環(huán)分布的預(yù)處理，實現(xiàn)了OpenACC循環(huán)映射、數(shù)據(jù)傳輸?shù)膬?yōu)化以及函數(shù)調(diào)用的從核化等幾方面程序變換，對于1.4°分辨率得到了良好的性能加速。本文面向100、1 000進(jìn)程的結(jié)點并行，考慮0.5°和0.25°高分辨率下的大規(guī)模眾核并行。

4.2 濾波階段的優(yōu)化評估與分析

由前可知，Tend_in應(yīng)用包含了stencil階段和濾波(filter)階段2個階段。其中，濾波階段的并行度很高，其AceMesh版本與OpenACC版本的行為差異很大，AceMesh相對加速比很高，因此本小節(jié)單獨對這個階段代碼的并行優(yōu)化方法進(jìn)行對比分析。

濾波階段主要有2個AceMesh優(yōu)化：影響任務(wù)映射的是單層圖和嵌套圖2種并行方式，影響通信并行度的是tag共享和tag私有2個不同的通信標(biāo)簽處理方式。以上2種優(yōu)化組合起來總共形成4個獨立的AceMesh版本，其特征描述如表2所示。表2中，“子圖”特指單個高度的濾波計算對應(yīng)的任務(wù)圖。簡單分析后可以知道，這4個版本按照性能由低到高的順序是tag共享的嵌套圖、tag共享的單層圖、tag私有的單層圖、tag私有的嵌套圖。

Table 2 Description of four AceMesh versions of the filter phase表2 Filter階段4個AceMesh版本的介紹

下面對比4個AceMesh版本相對OpenACC的加速比，由于tag共享的嵌套圖性能太差，本文沒有對其進(jìn)行測試，在此只展現(xiàn)后面3個版本的性能提升的情況。每種進(jìn)程數(shù)下，OpenACC和AceMesh版本采用相同的進(jìn)程網(wǎng)格分解方式(都選擇AceMesh版整體性能最優(yōu)的進(jìn)程分解方式)，給出AceMeshi/OpenACC的相對加速比，其中AceMeshi表示AceMesh的第i個優(yōu)化版本，i=2,3,4。圖2所示為0.5°分辨率下，16～1 024進(jìn)程時，3種AceMesh并行版本的相對加速比的差異。

Figure 2 AceMeshi/OpenACC speedups in the filter phase,i∈{2,3,4}圖2 濾波階段3個AceMesh優(yōu)化版本的AceMeshi/OpenACC加速比,i∈{2,3,4}

首先，tag共享的單層圖也實現(xiàn)了子圖之間一定程度的并行性。由于多個高度間屬于同一層任務(wù)依賴圖，當(dāng)一個高度的任務(wù)圖推進(jìn)由于通信等待而被阻塞時，計算資源可以分配給其他高度濾波的就緒任務(wù)，從而不同高度間的任務(wù)也存在并行機(jī)會，也挖掘了通信和計算的重疊。相比OpenACC只能挖掘循環(huán)內(nèi)不同緯度間的并行性，tag共享的單層圖相比OpenACC平均加速1.18倍。

其次，作為一個典型的通信密集型應(yīng)用，Tend_lin在濾波階段尤其表現(xiàn)出計算量小、通信密集的特點，通信之間的并行性、通信-計算之間的并行性對性能影響很大。tag私有的單層圖因為消除了通信之間的序關(guān)系，在保證部分高度并行的同時，為多個通信的并行創(chuàng)造了條件，在各種進(jìn)程配置下都獲得了較為明顯的加速效果，把相對加速比進(jìn)一步提高，平均再提高了0.64倍。此時不同高度之間尚未完全并行，仍然存在性能提升空間。

最后，tag私有的嵌套圖由于實現(xiàn)了不同高度子圖間的完全并行，性能進(jìn)一步提高。在一些進(jìn)程配置下性能提升還比較顯著。該版本相比OpenACC得到了最高6.7倍的加速，以及超過3.9倍的平均加速。

4.3 全程序性能評估與分析

圖3給出了2種分辨率下Tend_lin全程序在“神威·太湖之光”平臺上的性能測試結(jié)果，無論在何種進(jìn)程下，Tend_lin測試程序的AceMesh版的總性能都優(yōu)于OpenACC的。在0.5°分辨率下，AceMesh得到最高3倍，平均2.17倍的加速；在0.25°分辨率下，AceMesh得到最高2.86倍，平均1.98倍的加速，獲得了較好的加速效果。為理解性能數(shù)據(jù)，本文分別測量stencil階段和濾波階段。

Figure 3 AceMesh/OpenACC performance comparison圖3 AceMesh/OpenACC性能對比

從分項測試數(shù)據(jù)來看，stencil階段的計算時間較長，濾波階段的計算時間較短，所以stencil階段的性能對程序總性能提升效果影響更大。stencil階段相比濾波階段數(shù)據(jù)足跡和計算量較大，當(dāng)進(jìn)程數(shù)較少(小于32進(jìn)程)時，進(jìn)程內(nèi)并行度足夠，OpenACC同步式的推進(jìn)方式容易產(chǎn)生集中的DMA訪問，導(dǎo)致帶寬受限，而AceMesh的異步調(diào)度方式能夠有效地將DMA和計算錯峰調(diào)度，獲得更好的DMA性能。隨著進(jìn)程數(shù)的增加，進(jìn)程內(nèi)并行度和帶寬壓力的降低，DMA調(diào)度帶來的性能提升隨之下降。當(dāng)緯向的進(jìn)程數(shù)超過32時，進(jìn)程內(nèi)的并行度成為主要矛盾，此時將任務(wù)分塊因子從2調(diào)整到1，可以增加任務(wù)個數(shù)，降低從核的空置率，stencil階段依然能保持性能的穩(wěn)定。OpenACC只能在本地循環(huán)嵌套中并行，隨著進(jìn)程數(shù)增加，面臨的并行度不足問題越來越嚴(yán)重，即便使用collapse合并多層循環(huán)，達(dá)到的加速效果仍然有限。尤其作為通信密集型應(yīng)用，OpenACC并不能提供良好的通信優(yōu)化手段。不足的并行度和密集的通信成為限制OpenACC性能提升的關(guān)鍵，而這也是AceMesh并行的主要收益來源。0.5°分辨率下，stencil階段在1 024進(jìn)程時加速效果上升，其實是因為OpenACC并行度極低狀態(tài)下，空轉(zhuǎn)的線程產(chǎn)生了負(fù)面影響，導(dǎo)致絕對性能變差。

需要注意的是，濾波階段的數(shù)據(jù)與4.2節(jié)的數(shù)據(jù)存在很大的不同。4.2節(jié)為了評估不同優(yōu)化效果的影響，OpenACC和AceMesh采用了相同的進(jìn)程網(wǎng)格分解方式。但是，圖3中濾波階段的進(jìn)程網(wǎng)格分解方式與并行編程接口整體最優(yōu)性能下的進(jìn)程網(wǎng)格分解方式一致，從而其數(shù)據(jù)與圖2的數(shù)據(jù)存在很大的不同?？梢钥闯?，AceMesh加速效果隨著進(jìn)程數(shù)量的增加而增加，這是因為隨著進(jìn)程數(shù)量的增加，每個循環(huán)中剩下的工作越來越少，較低的并行度成為OpenACC的主要性能瓶頸。AceMesh則可以利用不同數(shù)組、不同高度和不同緯度之間隱含的3層并行性，相對充分地利用眾核計算資源，得到更好的加速效果。當(dāng)進(jìn)程數(shù)超過256時，AceMesh的加速效果急劇提高。

5 結(jié)束語

本文基于國產(chǎn)超算平臺“神威·太湖之光”，研究了大氣環(huán)流模式熱點Tend_lin的異構(gòu)并行優(yōu)化方法，指出了OpenACC在并行優(yōu)化方面的局限性，利用數(shù)據(jù)驅(qū)動的任務(wù)并行模型AceMesh對其進(jìn)行并行優(yōu)化，并探討了通信去串行化、嵌套任務(wù)圖的并行優(yōu)化方法。最后對2個不同的并行編程接口進(jìn)行性能對比，結(jié)果表明MPI+AceMesh并行獲得平均2倍左右的性能提升，本文還詳細(xì)分析了AceMesh性能收益的來源。

未來將進(jìn)行以下2項工作：(1)利用AceMesh并行編程環(huán)境把最新的IAP AGCM 5.0熱點移植到“神威·太湖之光”平臺上，同時把數(shù)據(jù)流并行擴(kuò)大到該模式代碼的更大范圍，以進(jìn)一步提高模擬的并行性能；(2)在AceMesh編譯系統(tǒng)中自動實現(xiàn)本研究中提出的編譯優(yōu)化技術(shù)，比如非整除塊的reshape優(yōu)化，提高編譯系統(tǒng)的性能。