顧文靜，孫 晨，王 彬

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

0 引 言

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

GRAPES(global/regional assimilation and prediction system)是中國(guó)氣象局自主研發(fā)的靜力/非靜力多尺度通用數(shù)值預(yù)報(bào)系統(tǒng),該系統(tǒng)是氣象與氣候研究的基礎(chǔ)和核心。對(duì)于一個(gè)大型、先進(jìn)的數(shù)值預(yù)報(bào)系統(tǒng)而言，并行計(jì)算已成為其必需的結(jié)構(gòu)特征[3]。

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

以GRAPES模式物理過(guò)程中的長(zhǎng)波輻射為加速研究對(duì)象，依托國(guó)家超級(jí)計(jì)算無(wú)錫中心的“神威·太湖之光”高性能計(jì)算機(jī)系統(tǒng)，設(shè)計(jì)一種基于OpenACC的加速算法，以實(shí)現(xiàn)利用更精簡(jiǎn)的指令來(lái)優(yōu)化代碼。

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

2016年6月20日，新一期全球超級(jí)計(jì)算機(jī)500強(qiáng)榜單在德國(guó)法蘭克福舉辦的國(guó)際超算大會(huì)(ISC)上公布，“神威·太湖之光”超級(jí)計(jì)算機(jī)榮登榜首?！吧裢ぬ狻庇蓢?guó)家并行計(jì)算機(jī)工程技術(shù)研究中心研制，是全球首臺(tái)運(yùn)行速度超過(guò)10億億次/秒的超級(jí)計(jì)算機(jī)，峰值性能高達(dá)12.54億億次/秒，持續(xù)性能達(dá)到9.3億億次/秒。該系統(tǒng)部署在國(guó)家超級(jí)計(jì)算無(wú)錫中心，系統(tǒng)由40個(gè)運(yùn)算機(jī)柜和8個(gè)網(wǎng)絡(luò)機(jī)柜組成。打開柜門，4個(gè)由32個(gè)節(jié)點(diǎn)板組成的超節(jié)點(diǎn)分布其中，每個(gè)節(jié)點(diǎn)板上包含4個(gè)節(jié)點(diǎn)卡，每個(gè)節(jié)點(diǎn)卡由兩個(gè)裝有“申威26010”眾核處理器和32 GB的DDR3內(nèi)存節(jié)點(diǎn)組成，一臺(tái)機(jī)柜就有1 024塊處理器。節(jié)點(diǎn)之間通過(guò)基于PCI-E3.0(外設(shè)部件互聯(lián)標(biāo)準(zhǔn)擴(kuò)展3.0版)的神威高速網(wǎng)絡(luò)進(jìn)行互聯(lián)。

“神威·太湖之光”系統(tǒng)全部采用“申威26010”眾核處理器，架構(gòu)如圖1所示。

圖1 “申威26010”處理器架構(gòu)

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

“神威·太湖之光”的冷卻系統(tǒng)采用計(jì)算節(jié)點(diǎn)板上全封閉式循環(huán)水冷技術(shù)和定制化的液體水冷單元，功耗比達(dá)到每瓦特60.51億次運(yùn)算。

2 OpenACC*編程模型

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

“神威·太湖之光”計(jì)算機(jī)系統(tǒng)中的OpenACC*語(yǔ)言是在OpenACC2.0文本的基礎(chǔ)上，針對(duì)申威26010眾核處理器結(jié)構(gòu)特點(diǎn)進(jìn)行適當(dāng)?shù)木?jiǎn)和擴(kuò)充而來(lái)的。

2.1 執(zhí)行模型

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

圖2 OpenACC*執(zhí)行模型

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

OpenACC*支持三級(jí)并行機(jī)制：gang、worker、vector。gang是粗粒度并行，在加速設(shè)備上可以啟動(dòng)一定數(shù)量的gang。worker是細(xì)粒度并行，每個(gè)gang內(nèi)包含有一定數(shù)量的worker。vector并行是在worker內(nèi)通過(guò)SIMD或向量操作的指令級(jí)并行[1]。

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

2.2 存儲(chǔ)模型

在CUDA和OpenCL等低層次加速器編程語(yǔ)言中，主機(jī)和加速器存儲(chǔ)器分離的概念非常明確，內(nèi)存間移動(dòng)數(shù)據(jù)的語(yǔ)句占據(jù)用戶大部分代碼[7]。

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

圖3 OpenACC*存儲(chǔ)模型

OpenACC標(biāo)準(zhǔn)中的數(shù)據(jù)管理功能將不產(chǎn)生實(shí)際的作用，但直接訪問(wèn)主存往往帶來(lái)性能損失，因此需要充分利用從核中的高速緩沖SPM，提升數(shù)據(jù)訪問(wèn)效率。在“神威·太湖之光”中，OpenACC*對(duì)Open-ACC標(biāo)準(zhǔn)所做的主要功能延伸和語(yǔ)法擴(kuò)展就是為了解決共享內(nèi)存架構(gòu)下片內(nèi)高速存儲(chǔ)空間的使用問(wèn)題。

“申威26010”中，從一個(gè)加速線程的角度來(lái)看，其可見的數(shù)據(jù)空間有三種。

(1)主線程數(shù)據(jù)空間：位于主存，主線程的內(nèi)存空間對(duì)其創(chuàng)建的加速線程直接可見，加速線程可以直接訪問(wèn)相應(yīng)的數(shù)據(jù)；對(duì)于主線程創(chuàng)建的多個(gè)加速線程而言，這部分空間是共享的；程序中在加速區(qū)外定義的變量，均位于該空間內(nèi)。

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

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

3 長(zhǎng)波輻射過(guò)程OpenACC優(yōu)化方案

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

3.1 輻射過(guò)程的構(gòu)成

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

輻射過(guò)程的函數(shù)包結(jié)構(gòu)如圖4所示。

圖4 長(zhǎng)波輻射過(guò)程的函數(shù)包結(jié)構(gòu)

3.2 并行優(yōu)化方案

并行化的唯一目的是充分利用硬件資源來(lái)提高程序運(yùn)行速度，縮短運(yùn)行時(shí)間。程序中最耗時(shí)間的是循環(huán)，計(jì)算并行化的目標(biāo)是將循環(huán)迭代步分散到多個(gè)不同的線程上執(zhí)行，這些線程運(yùn)行在多個(gè)加速器核心，從而將計(jì)算任務(wù)由CPU轉(zhuǎn)移到加速器上，減輕CPU的負(fù)擔(dān)[1]。

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

3.2.1 RRTM模塊優(yōu)化方案

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

圖5 循環(huán)內(nèi)移代碼結(jié)構(gòu)示意圖

3.2.2 mcica_subcol_gen_lw函數(shù)優(yōu)化方案

經(jīng)測(cè)試分析，函數(shù)中各進(jìn)程負(fù)載不均衡，通過(guò)插樁確定耗時(shí)主要集中在兩個(gè)分塊，一塊是隨機(jī)數(shù)生成、另一塊是數(shù)組計(jì)算。隨機(jī)數(shù)部分調(diào)配到從核運(yùn)行，并行生成0～1之間的隨機(jī)數(shù)，效率大幅提升。數(shù)組計(jì)算部分包含6個(gè)循環(huán)，考慮到各循環(huán)結(jié)構(gòu)和各層閾值均不一致，且沒有緊嵌套關(guān)系，不適合進(jìn)行循環(huán)合并，考慮直接添加OpenACC加速指示語(yǔ)句。下面取其中三個(gè)做實(shí)例分析。

(1)標(biāo)量、數(shù)組拷貝。

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

本段核心段代碼中循環(huán)體結(jié)構(gòu)為二維循環(huán)，外層循環(huán)閾值為[1,61]，內(nèi)層循環(huán)閾值為[1,90]，計(jì)算數(shù)組大小，LDM空間可以滿足，考慮直接進(jìn)行眾核加速，OpenACC編譯指示語(yǔ)句如下：

!$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兩個(gè)標(biāo)量放在每個(gè)加速線程的LDM中，其性質(zhì)是加速線程私有的，運(yùn)算過(guò)程中不做數(shù)據(jù)傳輸。程序中訪問(wèn)的數(shù)組有6個(gè)，rei、rel、play、reicmcl、relqmcl、pmid，其中rei、rel、play是只讀的，reicmcl、relqmcl、pmid是只寫的，對(duì)于只讀的數(shù)據(jù)只需要拷入(copyin)，只寫的數(shù)據(jù)只需要拷出(copyout)。當(dāng)動(dòng)態(tài)數(shù)組在加速計(jì)算區(qū)中以數(shù)組的形式使用時(shí)，編譯器往往無(wú)法判斷出該動(dòng)態(tài)數(shù)組所指向數(shù)組的維度信息，進(jìn)而無(wú)法在設(shè)備存儲(chǔ)器中申請(qǐng)適當(dāng)大小的空間并拷貝數(shù)據(jù)。dimension暗示中所指出的數(shù)組維度信息幫助編譯器在加速計(jì)算區(qū)中申請(qǐng)?jiān)O(shè)備空間并拷貝數(shù)據(jù)。collapse(2)表明nm和km兩層循環(huán)合并，有利于提高并行度。

(2)數(shù)組轉(zhuǎn)置。

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

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

調(diào)整后，ilev-循環(huán)對(duì)應(yīng)CDF數(shù)組的第三維，isubcol-循環(huán)對(duì)應(yīng)CDF數(shù)組的第一維，i-循環(huán)對(duì)應(yīng)CDF數(shù)組的第二維，OpenACC編譯指示語(yǔ)句如下：

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

……….//源代碼

!$ACC END PARALLEL LOOP

(3)循環(huán)分塊。

本段核心代碼中，數(shù)組拷貝需要的LDM空間過(guò)大，處理器無(wú)法滿足需求，考慮對(duì)循環(huán)進(jìn)行分塊處理，求解過(guò)程如下：

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

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

②計(jì)算次外層循環(huán)(i-循環(huán))的分塊大小。先令該層的分塊大小為l，這時(shí)三個(gè)數(shù)組的第2維按塊大小1劃分，三個(gè)數(shù)組數(shù)據(jù)塊總大小為(1*1+1*1+140*1*1)*8=1 136字節(jié)，遠(yuǎn)小于局存大小。為充分利用存儲(chǔ)空間，應(yīng)增大分塊值。當(dāng)該層的分塊值增大到58時(shí)，三個(gè)數(shù)組數(shù)據(jù)塊總大小為(58*1+58*1+140*58*1)*8=65 888，大于局存空間，達(dá)到臨界值。所以該層分塊大小確定為57，算法結(jié)束[12]。

該算法得到的循環(huán)分塊方案為(1，57，0)，即最外層按塊大小1劃分，次外層按塊大小57劃分，最內(nèi)層不劃分[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數(shù)學(xué)函數(shù)庫(kù)

在OpenACC加速指導(dǎo)語(yǔ)句的基礎(chǔ)上，使用ldm數(shù)學(xué)函數(shù)庫(kù)。經(jīng)測(cè)試，ldm數(shù)學(xué)函數(shù)精度與默認(rèn)庫(kù)保持一致，計(jì)算性能提升30%～80%不等，進(jìn)一步提升了GRAPES模式整體運(yùn)算效率。

3.2.4 結(jié)果與分析

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

表1 優(yōu)化結(jié)果對(duì)比

4 結(jié)束語(yǔ)

隨著編譯器的進(jìn)一步優(yōu)化和硬件技術(shù)的發(fā)展，OpenACC加速與CUDA等在底層技術(shù)實(shí)現(xiàn)上的差距將越來(lái)越小，而且支持OpenACC加速指令轉(zhuǎn)換將不僅僅針對(duì)CUDA設(shè)備，還包括其他更多廠商的硬件加速設(shè)備，從而能極大地提高OpenACC加速指令的普適性[13]。

文中以GRAPES模式長(zhǎng)波輻射模塊為加速對(duì)象，目前只應(yīng)用了為數(shù)不多的OpenACC指令，就能得到較高的加速比，優(yōu)化還有很大的提升空間。下一步將繼續(xù)深入研究源代碼數(shù)據(jù)結(jié)構(gòu)，加強(qiáng)數(shù)據(jù)預(yù)處理，減少?gòu)暮嗽L主存操作，充分利用從核中高速緩存SPM，提升數(shù)據(jù)訪問(wèn)效率；同時(shí)嘗試多種OpenACC加速指令組和使用，進(jìn)一步提升加速效果。

參考文獻(xiàn):

[1] 何滄平.OpenACC并行編程實(shí)戰(zhàn)[M]．北京:機(jī)械工業(yè)出版社，2017．

[2] 曾文權(quán)，胡玉貴，何擁軍，等．一種基于OPENACC的GPU加速實(shí)現(xiàn)高斯模糊算法[J]．計(jì)算機(jī)技術(shù)與發(fā)展，2013，23(7)：147-150．

[3] 伍湘君，金之雁，陳德輝，等．新一代數(shù)值預(yù)報(bào)模式GR-APES的并行計(jì)算方案設(shè)計(jì)與實(shí)現(xiàn)[J]．計(jì)算機(jī)研究與發(fā)展，2007，44(3)：510-515．

[4] 伍湘君，金之雁，黃麗萍，等．GRAPES模式軟件框架與實(shí)現(xiàn)[J]．應(yīng)用氣象學(xué)報(bào)，2005，16(4)：539-546．

[5] 陳德輝，沈?qū)W順．新一代數(shù)值預(yù)報(bào)系統(tǒng)GRAPES研究進(jìn)展[J]．應(yīng)用氣象學(xué)報(bào)，2006，17(6)：773-777．

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

[7] 劉文志.并行編程方法與優(yōu)化實(shí)踐[M]．北京:機(jī)械工業(yè)出版社，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數(shù)值預(yù)報(bào)系統(tǒng)中RRTM模塊的并行化研究[J].計(jì)算機(jī)科學(xué)，2012，39(6A)：370-374.

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

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

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

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