基于OpenACC的高性能計算并行優化研究與應用

顧文靜，孫 晨，王 彬

(國家氣象信息中心 高性能計算室，北京 100081)

0 引 言

近年來，面對浩如煙海的數據，CPU已力不從心，并行計算技術日趨成熟，世界領先的超級計算機都裝備大量的GPU加速器或眾核處理器[1]。然而，要在GPU硬件上實現加速需要通過對底層的API進行編程來實現，程序編寫復雜、難度大，且容易形成高度依賴特定設備的代碼。因此，一種基于指令制導計算的編程模型OpenACC應運而生[2]。OpenACC的編程機制是在源程序中添加少量編譯標識，編譯器根據作者的意圖自動產生低級語言代碼，無需學習新的編程語言和加速器硬件知識。只需添加少量編譯標識，不破壞原碼，開發容易，既可并行執行又可恢復串行執行；硬件更新時，不需要修改代碼，重新編譯即可。OpenACC標準制定時就考慮了當前及將來多種加速器產品，同一份代碼可以在多種加速器設備上編譯、運行、無成本切換硬件平臺[1]。

GRAPES(global/regional assimilation and prediction system)是中國氣象局自主研發的靜力/非靜力多尺度通用數值預報系統,該系統是氣象與氣候研究的基礎和核心。對于一個大型、先進的數值預報系統而言，并行計算已成為其必需的結構特征[3]。

GRAPES系統的核心部分包括模式動力框架和物理過程，動力框架中計算時間最長的為半拉格朗日計算和求解Helmholtz大型線性代數方程組，物理過程中耗時較大的為輻射和微物理過程，其中輻射過程計算耗時占到GRAPES模式的40%，所以提升輻射過程的性能和計算效率，對整個GRAPES模式的性能提升具有重要科研價值和實際意義[4-5]。

以GRAPES模式物理過程中的長波輻射為加速研究對象，依托國家超級計算無錫中心的“神威·太湖之光”高性能計算機系統，設計一種基于OpenACC的加速算法，以實現利用更精簡的指令來優化代碼。

1 “神威·太湖之光”概況

2016年6月20日，新一期全球超級計算機500強榜單在德國法蘭克福舉辦的國際超算大會(ISC)上公布，“神威·太湖之光”超級計算機榮登榜首?！吧裢ぬ狻庇蓢也⑿杏嬎銠C工程技術研究中心研制，是全球首臺運行速度超過10億億次/秒的超級計算機，峰值性能高達12.54億億次/秒，持續性能達到9.3億億次/秒。該系統部署在國家超級計算無錫中心，系統由40個運算機柜和8個網絡機柜組成。打開柜門，4個由32個節點板組成的超節點分布其中，每個節點板上包含4個節點卡，每個節點卡由兩個裝有“申威26010”眾核處理器和32 GB的DDR3內存節點組成，一臺機柜就有1 024塊處理器。節點之間通過基于PCI-E3.0(外設部件互聯標準擴展3.0版)的神威高速網絡進行互聯。

“神威·太湖之光”系統全部采用“申威26010”眾核處理器，架構如圖1所示。

圖1 “申威26010”處理器架構

該芯片主頻145 GHz，由4個核組組成，每個核組包含一個主核和一個從核簇，每個從核簇由64個從核組成，共260個處理器核心。同一個核組內的主核和64個從核共享主存，每個從核有獨立的、容量為64 KB的高速緩存(SPM)，支持DMA(direct memory access)的方式在主存和SPM之間傳輸數據。主核作為處理器核心，可以進行通信、IO、計算等操作，同一核組內的64個從核作為加速計算部件，用來加速主核代碼中的計算密集部分。

“神威·太湖之光”的冷卻系統采用計算節點板上全封閉式循環水冷技術和定制化的液體水冷單元，功耗比達到每瓦特60.51億次運算。

2 OpenACC*編程模型

OpenACC是由OpenACC組織(www.openacc-standard.org)于2011年推出的眾核加速編程語言，以編譯指示的方式提供眾核編程所需的語言功能，其主要目的是降低眾核編程的難度。用于C/C++和Fortran的OpenACC API和指令負責把底層的GPU任務交給編譯的同時，又提供跨系統、跨主機CPU和加速設備的支持[6]。

“神威·太湖之光”計算機系統中的OpenACC*語言是在OpenACC2.0文本的基礎上，針對申威26010眾核處理器結構特點進行適當的精簡和擴充而來的。

2.1 執行模型

OpenACC*程序的執行模型是在host(主處理器)的指導下，host和device(加速設備)協作的加速執行模型，如圖2所示。

圖2 OpenACC*執行模型

程序首先在host上啟動運行，以一個主線程串行執行，或者通過使用OpenMP或MPI等編程接口以多個主線程并行執行，計算密集區域則在主線程的控制下作為kernel(加速任務)被加載到加速設備上執行。kernel的執行過程包括：在設備內存上分配所需數據空間；加載kernel代碼(包含kernel參數)至device；kernel將所需數據從主存傳輸至設備內存；等待數據傳輸完成；device進行計算并將結果傳送回主存；釋放設備上的數據空間等。大多數情況下，host可以加載一系列kernels，并在加速設備上逐個執行[1]。

OpenACC*支持三級并行機制：gang、worker、vector。gang是粗粒度并行，在加速設備上可以啟動一定數量的gang。worker是細粒度并行，每個gang內包含有一定數量的worker。vector并行是在worker內通過SIMD或向量操作的指令級并行[1]。

在“申威26010”中，一個運算控制核心(簡稱主核)僅控制一個運算核心陣列(加速設備，簡稱從核陣列)的運行，每個運算核心陣列內有64個運算核心(簡稱從核)，每個運算核心可以運行一個加速線程。默認情況下，64個加速線程被組織成64個gang、每個gang內一個worker、worker內可以vector并行的邏輯視圖。通常情況下，盡量用滿64個加速線程可以獲得較好的運行性能。

2.2 存儲模型

在CUDA和OpenCL等低層次加速器編程語言中，主機和加速器存儲器分離的概念非常明確，內存間移動數據的語句占據用戶大部分代碼[7]。

“申威26010”中從核和主核共享內存，從核可直接訪問主存空間，并在從核內提供加速線程私有的高速緩沖(local data memory,LDM)，加速計算需要存放到LDM的數據由從核控制傳輸。存儲模型如圖3所示。

圖3 OpenACC*存儲模型

OpenACC標準中的數據管理功能將不產生實際的作用，但直接訪問主存往往帶來性能損失，因此需要充分利用從核中的高速緩沖SPM，提升數據訪問效率。在“神威·太湖之光”中，OpenACC*對Open-ACC標準所做的主要功能延伸和語法擴展就是為了解決共享內存架構下片內高速存儲空間的使用問題。

“申威26010”中，從一個加速線程的角度來看，其可見的數據空間有三種。

(1)主線程數據空間：位于主存，主線程的內存空間對其創建的加速線程直接可見，加速線程可以直接訪問相應的數據；對于主線程創建的多個加速線程而言，這部分空間是共享的；程序中在加速區外定義的變量，均位于該空間內。

(2)加速線程私有空間：位于主存，每個加速線程有獨立的私有空間，程序中使用private、firstprivate子句修飾的變量將存放于該空間內。

(3)加速線程本地空間：位于設備內存，每個加速線程有獨立的本地空間(LDM)，本地空間的訪問性能是三種空間中最高的。程序中使用local、copy等數據子句修飾的變量將由編譯系統控制全部或局部存放于LDM內。主存與本地空間的數據交互由加速線程控制。

3 長波輻射過程OpenACC優化方案

GRAPES模式的輻射模型采用歐洲中期天氣預報中心(ECMWF)的長短波輻射方案，該方案將整個大氣層劃分成水平方向和垂直方向的三維網格，水平方向為橫跨地球表面的二維網格，垂直方向則表示大氣的層數[8-9]。數值天氣預報模式是一種非線性化離散計算模式，其計算量巨大[10]，最初的都是通過CPU計算實現，GPU出現之后，雖然可以把該算法移植到GPU上并行執行，但編程難度大且易出錯。因此，文中設計一種基于OPENACC加速方法，以實現通過簡化的編程算法和簡潔的Fortran語言代碼完成對長波輻射過程并行處理過程。

3.1 輻射過程的構成

輻射過程是GRAPES系統中一個重要的物理過程，包含云結構的描述和輻射算法(RRTM模塊)兩部分。其中云對調解輻射平衡起著至關重要的作用[11]，該系統在云產生器的基礎上利用McICA進行輻射計算。輻射過程首先調用mcica_subcol_lw函數實現云結構的描述，之后調用RRTM模塊進行輻射計算。RRTM模塊初始化函數rrtmg_lw_init讀取數據文件中的輸入數據；然后調用rrtmg_lw_rad函數計算輻射傳輸的過程；最后調用輻射模型子函數rrtm_lw，rrtm_lw子函數一次計算一個垂直列的輻射傳輸值[9]。它包含以下幾個基本步驟：inatm：準備大氣剖面；cldprmc：計算云層光學深度；setcoef：計算這種大氣剖面下輻射傳輸算法所需的各種量；taumol：計算氣體光學深度和普朗克分數；rtrnmc：計算晴空和云層的輻射傳輸值。其中，taumol中的16個子函數taugb，計算了全部16個長波波域譜帶的氣態光學厚度和普朗克(Planck)函數。

輻射過程的函數包結構如圖4所示。

圖4 長波輻射過程的函數包結構

3.2 并行優化方案

并行化的唯一目的是充分利用硬件資源來提高程序運行速度，縮短運行時間。程序中最耗時間的是循環，計算并行化的目標是將循環迭代步分散到多個不同的線程上執行，這些線程運行在多個加速器核心，從而將計算任務由CPU轉移到加速器上，減輕CPU的負擔[1]。

結合程序代碼的分析，確定需要眾核加速的代碼段，并根據源代碼數據結構進行預處理，對相關數組和循環進行調整，以適合OpenACC加速。優化的基本思想是確定將程序中的標量和關鍵數組盡可能多地放到局存(LDM)中，并設法提高程序中數據的傳輸效率，將需要使用眾核加速的循環的前后使用OpenACC加速編譯指示進行標注，其中PARALLEL表明該部分代碼是需要加速執行的并行代碼；LOOP表示下面緊跟著的循環需要在加速線程間進行并行劃分。

3.2.1 RRTM模塊優化方案

rrtm_lw函數包含一層循環，涵蓋inatm、cldprmc、setcoef、taumol和rtrnmc五個子函數調用，每個子函數內部有多個do循環分塊，考慮最外層的iplon循環閾值為[1,90]，直接進行眾核并行，從核使用率僅為41%，沒有充分發揮核組計算能力。且函數間調用關系復雜，如果在該層循環外層添加OpenACC加速指示語句，循環內部的所有函數將全部被調配到從核運行，即使添加routine函數將部分從核函數的棧數組放到主存，但由于子函數中變量眾多，對從核空間需求較高，最終仍難避免棧的容量超出局存(LDM)，造成模式運行出錯。結合輻射過程的計算模型為單柱計算模型，核心迭代針對每個柱進行多個物理過程計算，分析到各子函數間存在數據依賴關系，為保證數據正確性且增加從核運算并行度，考慮對rrtmg_lw函數進行拆分并進行循環內移。內移工作從rrtmg_lw函數調用的函數inatm開始，依次對cldprmc、setcoef、taumol和rtrnmc子函數進行循環內移，代碼結構示意如圖5所示。

圖5 循環內移代碼結構示意圖

3.2.2 mcica_subcol_gen_lw函數優化方案

經測試分析，函數中各進程負載不均衡，通過插樁確定耗時主要集中在兩個分塊，一塊是隨機數生成、另一塊是數組計算。隨機數部分調配到從核運行，并行生成0～1之間的隨機數，效率大幅提升。數組計算部分包含6個循環，考慮到各循環結構和各層閾值均不一致，且沒有緊嵌套關系，不適合進行循環合并，考慮直接添加OpenACC加速指示語句。下面取其中三個做實例分析。

(1)標量、數組拷貝。

do nm=1, nlay

do km=1, ncol

reicmcl(km,nm)=rei(km,nm)

relqmcl(km,nm)=rel(km,nm)

pmid(km,nm)=play(km,nm)*1.e2_rb

enddo

enddo

本段核心段代碼中循環體結構為二維循環，外層循環閾值為[1,61]，內層循環閾值為[1,90]，計算數組大小，LDM空間可以滿足，考慮直接進行眾核加速，OpenACC編譯指示語句如下：

!$ACC PARALLEL LOOP local(nm,km)

!$ACC copyin(rei,rel,play) copyout(reicmcl,relqmcl,pmid)

!$ACC annotate(dimension(rei(ncol,nlay),rel(ncol,nlay),play(ncol,nlay),reicmcl(ncol,nlay),relqmcl(ncol,nlay),pmid(ncol,nlay))) collapse(2)

………//源代碼

!$ACC END PARALLEL LOOP

其中，local(nm,km)表明將nm、km兩個標量放在每個加速線程的LDM中，其性質是加速線程私有的，運算過程中不做數據傳輸。程序中訪問的數組有6個，rei、rel、play、reicmcl、relqmcl、pmid，其中rei、rel、play是只讀的，reicmcl、relqmcl、pmid是只寫的，對于只讀的數據只需要拷入(copyin)，只寫的數據只需要拷出(copyout)。當動態數組在加速計算區中以數組的形式使用時，編譯器往往無法判斷出該動態數組所指向數組的維度信息，進而無法在設備存儲器中申請適當大小的空間并拷貝數據。dimension暗示中所指出的數組維度信息幫助編譯器在加速計算區中申請設備空間并拷貝數據。collapse(2)表明nm和km兩層循環合并，有利于提高并行度。

(2)數組轉置。

do i=1, ncol

do isubcol=1, nsubcol

do ilev=2,nlay

if(CDF(isubcol, i, ilev-1)>1._rb-cldf(i,ilev-1) )

then

CDF(isubcol,i,ilev)=CDF(isubcol,i,ilev-1)

else

CDF(isubcol,i,ilev)=CDF(isubcol,i,ilev) *

(1._rb-cldf(i,ilev-1))

endif

enddo

enddo

enddo

本段源碼核心段是矩陣乘運算代碼段，由于數組CDF的訪問與ilev相關，因此ilev-循環不能分塊(tile)處理，因此將ilev循環調到最內層，循環的調整導致CDF數組的訪問不連續。不連續的數據訪問會導致數據傳輸效率較低，而swapout的作用是對CDF數組的原始數據按照指定的維度順序進行轉置，使用轉置后的數據代替原始CDF數組的訪問，這樣可以使不連續的數據訪問變成連續的數據訪問。

調整后，ilev-循環對應CDF數組的第三維，isubcol-循環對應CDF數組的第一維，i-循環對應CDF數組的第二維，OpenACC編譯指示語句如下：

!$ACC PARALLEL LOOP swap(CDF(dimension order:3,1,2)) copyin(cldf) annotate(dimension(cldf(ncol,nlay))) collapse(2)

……….//源代碼

!$ACC END PARALLEL LOOP

(3)循環分塊。

本段核心代碼中，數組拷貝需要的LDM空間過大，處理器無法滿足需求，考慮對循環進行分塊處理，求解過程如下：

clwp、ciwp和tauc數組的數據類型均為雙精度浮點，數據大小為8字節，計算該循環使用三個數組的大小，其值為(90*61+90*61+140*90*61)*8=6 236 640字節，遠超出局存大小(65 536字節)，須分塊后才能裝入從核。

①計算最外層循環(ilev-循環)的分塊大小。令該層的分塊大小為1，這時三個數組數據塊的總大小為(90*1+90*1+140*90*1)*8=102 240字節，仍超過了局存大小。表明對數組該維的最小分塊仍會使局存溢出，需要對下一維進行劃分，取該層的分塊值為1[12]。

②計算次外層循環(i-循環)的分塊大小。先令該層的分塊大小為l，這時三個數組的第2維按塊大小1劃分，三個數組數據塊總大小為(1*1+1*1+140*1*1)*8=1 136字節，遠小于局存大小。為充分利用存儲空間，應增大分塊值。當該層的分塊值增大到58時，三個數組數據塊總大小為(58*1+58*1+140*58*1)*8=65 888，大于局存空間，達到臨界值。所以該層分塊大小確定為57，算法結束[12]。

該算法得到的循環分塊方案為(1，57，0)，即最外層按塊大小1劃分，次外層按塊大小57劃分，最內層不劃分[12]，OpenACC代碼如下：

!$ACC PARALLEL LOOPcopyin(clwp,ciwp,tauc) copyout(clwp_stoch(*,*,i),ciwp_stoch(*,*,i),tauc_stoch(*,*,i))

do ilev=1,nlay

!$ACC LOOP tile(57)

do i=1, ncol

do isubcol=1, nsubcol

..........

clwp_stoch(isubcol,i,ilev)=clwp(i,ilev)

ciwp_stoch(isubcol,i,ilev)=ciwp(i,ilev)

tauc_stoch(isubcol,i,ilev)=tauc(isubcol,i,ilev)

..........

enddo

enddo

enddo

!$ACC END PARALLEL LOOP

3.2.3 使用ldm數學函數庫

在OpenACC加速指導語句的基礎上，使用ldm數學函數庫。經測試，ldm數學函數精度與默認庫保持一致，計算性能提升30%～80%不等，進一步提升了GRAPES模式整體運算效率。

3.2.4 結果與分析

實驗采用分辨率為0.5°×0.5°的GRAPES全球模式長波輻射方案，積分步長為36，優化選項使用-O2。由于cldprmc、setcoef耗時較少，沒有進行優化。加速后，對比輸出的全球長波輻射通量與模式輸出值，經驗證誤差，相對誤差的量級在10-5～10-4之間，在允許范圍內。優化結果對比如表1所示。

表1 優化結果對比

4 結束語

隨著編譯器的進一步優化和硬件技術的發展，OpenACC加速與CUDA等在底層技術實現上的差距將越來越小，而且支持OpenACC加速指令轉換將不僅僅針對CUDA設備，還包括其他更多廠商的硬件加速設備，從而能極大地提高OpenACC加速指令的普適性[13]。

文中以GRAPES模式長波輻射模塊為加速對象，目前只應用了為數不多的OpenACC指令，就能得到較高的加速比，優化還有很大的提升空間。下一步將繼續深入研究源代碼數據結構，加強數據預處理，減少從核訪主存操作，充分利用從核中高速緩存SPM，提升數據訪問效率；同時嘗試多種OpenACC加速指令組和使用，進一步提升加速效果。

參考文獻:

[1] 何滄平.OpenACC并行編程實戰[M]．北京:機械工業出版社，2017．

[2] 曾文權，胡玉貴，何擁軍，等．一種基于OPENACC的GPU加速實現高斯模糊算法[J]．計算機技術與發展，2013，23(7)：147-150．

[3] 伍湘君，金之雁，陳德輝，等．新一代數值預報模式GR-APES的并行計算方案設計與實現[J]．計算機研究與發展，2007，44(3)：510-515．

[4] 伍湘君，金之雁，黃麗萍，等．GRAPES模式軟件框架與實現[J]．應用氣象學報，2005，16(4)：539-546．

[5] 陳德輝，沈學順．新一代數值預報系統GRAPES研究進展[J]．應用氣象學報，2006，17(6)：773-777．

[6] 郭 妙，金之雁，周 斌.基于通用圖形處理器的GRAPES長波輻射并行方案[J].應用氣象學報，2012，23(3)：348-354．

[7] 劉文志.并行編程方法與優化實踐[M]．北京:機械工業出版社，2015．

[8] LI L,XUE W,RANJAN R,et al.A scalable Helmholtz solver in GRAPES over large-scale multicore cluster[J].Concurrency and Computation:Practice and Experience,2013，25(12)：1722-1737．

[9] 鄭 芳，許先斌，向冬冬,等.基于GPU的GRAPES數值預報系統中RRTM模塊的并行化研究[J].計算機科學，2012，39(6A)：370-374.

[10] 王 為，張悠慧，姚 駿,等.基于線性陣列處理器的GR-APES核心代碼優化[J].計算機學報,2013，36(10)：2053-2061.

[11] 荊現文,張 華.McICA云-輻射方案在國家氣候中心全球氣候模式中的應用與評估[J]．大氣科學，2012，36(5)：945-958.

[12] 李雁冰，趙榮彩，趙 博,等.面向異構多核處理器的循環分塊[J]．計算機工程與設計，2015，36(1)：168-173．

[13] FARBER R.Parallel programming with OpenACC[M]．Massachusetts:Morgan Kaufmann Publishers,2015.