王家潤(rùn)，謝海峰

(華北計(jì)算技術(shù)研究所基礎(chǔ)三部，北京 100083)

0 引 言

在國(guó)家“核高基”重大專項(xiàng)支持下，自主可控的軟硬件發(fā)展迅速。但國(guó)產(chǎn)基礎(chǔ)軟硬件技術(shù)成熟度不高、技術(shù)儲(chǔ)備不足、信息化生態(tài)環(huán)境不完善、體系化能力弱、整體安全性尚待提高等系列問題，還需要重點(diǎn)解決應(yīng)用遷移、系統(tǒng)適配、性能優(yōu)化等問題[1]。文獻(xiàn)[1]針對(duì)典型基礎(chǔ)軟硬件，優(yōu)化參數(shù)配置、語句性能、界面渲染、數(shù)據(jù)管理等；文獻(xiàn)[2]通過綜合運(yùn)用并發(fā)處理、數(shù)據(jù)緩存、分片處理和分層顯示等技術(shù)，對(duì)信息的產(chǎn)生、傳輸、處理、查詢及顯示的全過程進(jìn)行性能優(yōu)化；文獻(xiàn)[3]研究了嵌入式環(huán)境中的國(guó)產(chǎn)化GPU JM5400對(duì)漢字的顯示及大數(shù)據(jù)量圖形的顯示優(yōu)化，但與商用相比，性能差距比較大。

基于CPU多核的OpenMP[4,5]及基于GPU的CUDA[5]的并行優(yōu)化技術(shù)已在各領(lǐng)域中大量應(yīng)用。文獻(xiàn)[5]對(duì)現(xiàn)代處理器特性、串行代碼優(yōu)化、常用并行編程模型及設(shè)計(jì)方法等進(jìn)行了深入分析；文獻(xiàn)[6]研究了OpenMP在Kriging插值中的并行優(yōu)化；文獻(xiàn)[7]中將CUDA應(yīng)用在運(yùn)動(dòng)目標(biāo)檢測(cè)處理；文獻(xiàn)[8]研究了OpenMP在矢量空間數(shù)據(jù)拓?fù)渌惴ㄖ械牟⑿袃?yōu)化；文獻(xiàn)[9]研究了CUDA在電磁可視化中的性能加速；文獻(xiàn)[10]研究了GPU性能優(yōu)化的主要途徑并提出了一個(gè)可視化的性能指導(dǎo)模型。上述研究針對(duì)不同的實(shí)際問題進(jìn)行了并行優(yōu)化，但是處理過程相對(duì)復(fù)雜、多樣。

本文以離散數(shù)據(jù)格網(wǎng)化處理中的并行優(yōu)化研究為例，研究國(guó)產(chǎn)應(yīng)用性能提升技術(shù)。借鑒文獻(xiàn)[10]的性能指導(dǎo)模型思想，針對(duì)并行編程中的線程任務(wù)劃分這一共性難點(diǎn)，建立了一個(gè)用于指導(dǎo)任務(wù)劃分的線程任務(wù)分解通用處理模型。

1 軟硬件結(jié)合并行優(yōu)化策略及處理流程

1.1 軟硬件結(jié)合并行優(yōu)化策略

國(guó)產(chǎn)CPU與商用相比雖然還存在差距，但因?yàn)槠毡椴捎昧硕嗪思夹g(shù)，在一定程度上可彌補(bǔ)主頻不足的問題；國(guó)產(chǎn)GPU目前沒有面向桌面的產(chǎn)品，主要使用商用GPU，由于缺乏適配的驅(qū)動(dòng)程序，只能采用開源的驅(qū)動(dòng)程序替代，與商用性能差距較大，目前國(guó)產(chǎn)下GPU通用計(jì)算還無法實(shí)施。因此國(guó)產(chǎn)下性能優(yōu)化應(yīng)以國(guó)產(chǎn)CPU處理器提供的多核處理技術(shù)為主，以商用GPU為輔(主要用于對(duì)比驗(yàn)證及技術(shù)儲(chǔ)備)。

1.2 軟硬件結(jié)合并行優(yōu)化處理流程

(1)性能分析

通過性能測(cè)試工具Intel VTune等找出軟件中的性能瓶頸。

(2)并行處理

并行編程的思想是分治策略，即將大問題劃分為一些小問題，再把這些小問題交給相應(yīng)的處理單元并行處理。分解后的線程處理任務(wù)是否不相關(guān)(或相關(guān)性很低)，這是并行優(yōu)化的前提。

(3)對(duì)比驗(yàn)證

采用并行處理的軟件，復(fù)雜性及調(diào)試難度較高，建議先采用串行編程，再并行優(yōu)化，串行編程用于驗(yàn)證并行優(yōu)化算法的正確性及性能優(yōu)化對(duì)比。

(4)性能調(diào)優(yōu)

根據(jù)CPU及GPU處理器硬件的結(jié)構(gòu)特征，對(duì)并行優(yōu)化過程中的細(xì)節(jié)調(diào)優(yōu)(編譯器選項(xiàng)、數(shù)據(jù)分塊、內(nèi)存使用、負(fù)載平衡、調(diào)度優(yōu)化、緩存優(yōu)化等[5,11,12])。

2 離散數(shù)據(jù)格網(wǎng)化IDW算法及并行優(yōu)化設(shè)計(jì)

2.1 離散數(shù)據(jù)格網(wǎng)化IDW算法

離散數(shù)據(jù)格網(wǎng)化有廣泛的應(yīng)用，例如：氣象觀測(cè)站的分布是離散不規(guī)則的，需要先根據(jù)這些觀測(cè)站數(shù)據(jù)來插值出規(guī)則網(wǎng)格點(diǎn)上的屬性值(溫度、氣壓等)，基于規(guī)則網(wǎng)格，再采用等值線、填色圖等進(jìn)行可視化等。基于空間相近相似原理的反距離權(quán)重IDW插值算法是一種常用的方法。

反距離權(quán)重IDW插值算法[13]：

(1)輸入

離散采樣點(diǎn)數(shù)據(jù)(位置、屬性等)集合；待插值規(guī)則網(wǎng)格的范圍、行列數(shù)目等。

(2)過程

A：鄰近搜索(Searching)

針對(duì)需要構(gòu)建的規(guī)則網(wǎng)格中的每個(gè)格點(diǎn)，依據(jù)某種搜索規(guī)則，從離散采樣點(diǎn)數(shù)據(jù)集合中搜索與該格點(diǎn)相距較近的若干個(gè)離散采樣點(diǎn)，使用這些離散采樣點(diǎn)插值計(jì)算該格點(diǎn)的屬性值。目前有多種搜索規(guī)則：圓域、四方向、八方向、基于Delaunay三角網(wǎng)一階鄰近點(diǎn)、基于遮蔽程度的搜索等[13-15]，本文采用四方向搜索法。

B：插值計(jì)算(Computing)

針對(duì)規(guī)則網(wǎng)格中的每個(gè)格點(diǎn)搜索出的若干個(gè)離散采樣點(diǎn)，采用反距離權(quán)重IDW插值公式對(duì)這若干個(gè)離散采樣點(diǎn)屬性數(shù)據(jù)進(jìn)行計(jì)算，得到該網(wǎng)格格點(diǎn)的屬性值。IDW插值計(jì)算公式

其中，權(quán)重系數(shù)

其中，Z(S0)為S0處的預(yù)測(cè)值，N為計(jì)算過程中使用的鄰近采樣點(diǎn)數(shù)量，Zi為第i個(gè)采樣點(diǎn)的值，di為插值點(diǎn)S0到各采樣點(diǎn)的距離，2是距離的冪。

(3)輸出

插值計(jì)算出的各網(wǎng)格的格點(diǎn)屬性值。根據(jù)需要可進(jìn)一步進(jìn)行交叉驗(yàn)證，計(jì)算插值誤差等。

由IDW插值算法過程可知，規(guī)則格網(wǎng)的格點(diǎn)數(shù)目一般比較大、計(jì)算密集，每個(gè)格點(diǎn)的搜索及插值與鄰近的格點(diǎn)沒有依賴，相互獨(dú)立，非常符合并行化處理。下面遵循軟硬件結(jié)合的并行優(yōu)化策略，進(jìn)行并行優(yōu)化。

2.2 基于OpenMP的并行優(yōu)化PIDW算法

OpenMP使用Fork-Join并行編程模型，程序創(chuàng)建一個(gè)獨(dú)立的串行執(zhí)行的主線程，執(zhí)行到OpenMP標(biāo)識(shí)的并行域時(shí)，即進(jìn)入Fork過程：主線程會(huì)根據(jù)任務(wù)的劃分創(chuàng)建一組并行的子線程，然后并行域中的程序代碼在不同的線程中并行執(zhí)行；Join過程：當(dāng)主線程創(chuàng)建的并行子線程在并行域中執(zhí)行完任務(wù)，退出并行域時(shí)，OpenMP內(nèi)部同步各子線程，結(jié)束各子線程，返回到只有主線程運(yùn)行狀態(tài)。

基于Fork-Join并行編程模型的并行優(yōu)化PIDW算法：

(1)輸入

待插值網(wǎng)格格點(diǎn)數(shù)據(jù)、離散采樣點(diǎn)數(shù)據(jù)、網(wǎng)格范圍及行列數(shù)目等。

(2)過程

步驟1 設(shè)置OpenMP運(yùn)行參數(shù)：獲取CPU核數(shù)，根據(jù)核數(shù)設(shè)置要啟動(dòng)的子線程數(shù)。設(shè)置子線程數(shù)目與核數(shù)相同，可最大程度上利用各核提供的硬件能力，使各核具有較好的負(fù)載均衡。

步驟2 進(jìn)入OpenMP并行域(Fork)：根據(jù)子線程數(shù)目進(jìn)行任務(wù)劃分，確定每個(gè)線程負(fù)責(zé)的任務(wù)。子線程計(jì)算網(wǎng)格格點(diǎn)的屬性值(搜索鄰近采樣點(diǎn)，并采用IDW插值公式計(jì)算)。各子線程間的計(jì)算相對(duì)獨(dú)立，互不干擾。詳細(xì)可參見基于OpenMP的線程任務(wù)分解模型部分。

步驟3 退出OpenMP并行域(Join)：OpenMP并行域執(zhí)行完畢后，即結(jié)束并行計(jì)算模式，此時(shí)網(wǎng)格格點(diǎn)數(shù)據(jù)插值處理完畢。

(3)輸出

經(jīng)過插值處理的網(wǎng)格格點(diǎn)數(shù)據(jù)。

基于OpenMP實(shí)現(xiàn)的代碼示例：

//gridData：待插值網(wǎng)格格點(diǎn)數(shù)據(jù)，待插值網(wǎng)格列

//數(shù)、行數(shù)：gridNumX，girdNumY

//離散采樣點(diǎn)數(shù)據(jù)：sampleData

//獲取CPU處理器的核數(shù)目

intthreads=omp_get_Num_threads()；

//設(shè)置并行啟用threads個(gè)子線程，啟用CPU所有//核

omp_set_threads(threads)；

//進(jìn)入并行域(Fork)，啟用threads個(gè)子線程

//各線程相互獨(dú)立，并行處理計(jì)算任務(wù)

#pragmaompparallel//OpenMP并行域標(biāo)識(shí)

{

//獲取當(dāng)前子線程索引，索引范圍：0～threads-1

intthreadIdx=omp_get_thread_num();

//設(shè)置該子線程的循環(huán)計(jì)算任務(wù)：不越界

//(i

//(i+=threads)，不重復(fù)且保證所有的線程計(jì)算

//覆蓋所有的網(wǎng)格格點(diǎn)計(jì)算

for(inti=threadIdx;i

for(intj=0;j

//搜索與插值計(jì)算

SearchCompute(i,j,gridData,

gridNumX,gridNumY，sampleData);

}//退出并行域,結(jié)束threads個(gè)子線程，返回主線

//程(Join)，此時(shí)網(wǎng)格gridData的各格點(diǎn)數(shù)據(jù)已

//經(jīng)過插值處理，可用等值線跟蹤算法等進(jìn)一步處//理，可視化顯示等。

2.3 基于OpenMP的線程任務(wù)分解模型

OpenMP總共分配有N個(gè)子線程(在整個(gè)并行處理域內(nèi)，OpenMP會(huì)為每個(gè)子線程分配唯一的索引threadIdx，索引號(hào)范圍：0～N-1)，需要處理columnNum列、rowNum行的二維網(wǎng)格格點(diǎn)的計(jì)算任務(wù)(網(wǎng)格格點(diǎn)從左下角(0,0)開始，沿X、Y方向編號(hào))。

在一般情況下，子線程數(shù)目遠(yuǎn)遠(yuǎn)小于網(wǎng)格格點(diǎn)數(shù)目，因此，每個(gè)子線程需要負(fù)責(zé)多個(gè)格點(diǎn)的計(jì)算。基本劃分思想是將網(wǎng)格格點(diǎn)按列(或行)分解到N個(gè)組中，每個(gè)組由一個(gè)子線程負(fù)責(zé)計(jì)算，重點(diǎn)是建立線程索引與網(wǎng)格格點(diǎn)索引的映射關(guān)系，下面通過二重循環(huán)來完成每個(gè)子線程的任務(wù)分解：

intthreadIdx=omp_get_thread_num()；

for(inti=threadIdx;i

for(intj=0;j

GridNodeCompute(…);

即針對(duì)并行域中的當(dāng)前線程threadIdx，由該線程負(fù)責(zé)處理網(wǎng)格中從第i=threadIdx列(即將線程索引threadIdx映射到網(wǎng)格的第i列)開始到columnNum-1列范圍內(nèi)的各列格點(diǎn)的所有計(jì)算(j:0～rowNum)，其中列i的遞增間隔為線程總數(shù)N，原因是已分配有N個(gè)線程，前0～N-1列已分配到這N個(gè)子線程中并行處理，所以前N列無需重復(fù)計(jì)算，這是線程并行任務(wù)劃分的關(guān)鍵。i=i+N表示沿X方向繼續(xù)搜尋需要該線程計(jì)算的列，保證計(jì)算不遺漏。i

線程任務(wù)分解模型如圖1所示，共使用4個(gè)子線程：0～3，將網(wǎng)格中的前4列直接映射到4個(gè)線程(圖中左側(cè)細(xì)線矩形范圍內(nèi)的4列網(wǎng)格格點(diǎn))，其中線程0負(fù)責(zé)網(wǎng)格的列0與列4格點(diǎn)計(jì)算(參見圓符號(hào)標(biāo)記)，線程1、2與線程0處理類似，線程3只處理網(wǎng)格列3(參見矩形符號(hào)標(biāo)記)。

圖1 基于OpenMP的線程任務(wù)分解模型

線程任務(wù)分解模型的計(jì)算均衡(負(fù)載平衡)：一般情況下，首先每個(gè)線程參與計(jì)算，對(duì)應(yīng)使用了CPU的各核，對(duì)計(jì)算任務(wù)進(jìn)行了基于全部核的分解；其次，通過任務(wù)分解模型示意圖可看出，各線程之間最多相差網(wǎng)格的一列或一行格點(diǎn)，各線程的計(jì)算量相差不太大，整體來說，計(jì)算比較均衡，即線程任務(wù)分解模型具有較好的負(fù)載均衡。

2.4 基于CUDA的并行優(yōu)化PIDW算法

根據(jù)CUDA架構(gòu)的編程模型(HostDevice)[5]，基于CUDA的并行計(jì)算主要分為以下3個(gè)過程：將待處理數(shù)據(jù)由Host端(CPU)讀入Device端(GPU)顯存中；在Device端(GPU)中，進(jìn)行并行計(jì)算：采用多線程并行處理數(shù)據(jù)，并將處理結(jié)果存儲(chǔ)在顯存中；把Device端(GPU)中的處理過的數(shù)據(jù)讀回Host端(CPU)的內(nèi)存中。CUDA內(nèi)核函數(shù)運(yùn)行在GPU上，啟動(dòng)后CUDA中的每一個(gè)線程都將會(huì)同時(shí)并行地執(zhí)行內(nèi)核函數(shù)中的代碼。GPU是由多個(gè)流多處理器構(gòu)成，以塊(Block)為基本調(diào)度單元，因此對(duì)于流處理器較多的GPU，一次可以處理的塊(Block)更多，從而運(yùn)算速度更快。

使用CUDA編程主要注意：①線程的設(shè)置。獲取GPU硬件參數(shù)：流處理器及線程塊最大線程個(gè)數(shù)等參數(shù)及處理問題規(guī)模等。②線程的索引標(biāo)識(shí)。每個(gè)線程都有線程索引threadIdx、所在的塊索引blockIdx、所在塊的維數(shù)blockDim、所在的線程網(wǎng)格維數(shù)gridDim等CUDA內(nèi)置的信息(圖2)，通過這些信息可計(jì)算出線程在整個(gè)線程網(wǎng)格Grid中的索引位置，如在線程網(wǎng)格的X索引方向的索引位置

threadIdx.x+blockIdx.x*blockDim.x

基于CUDA架構(gòu)編程模型的并行優(yōu)化PIDW算法：

(1)輸入

待插值網(wǎng)格格點(diǎn)數(shù)據(jù)、離散采樣點(diǎn)數(shù)據(jù)、網(wǎng)格范圍及行列數(shù)目等。

(2)過程

步驟1 獲取GPU設(shè)備的硬件參數(shù)：GPU的流處理器(SP)個(gè)數(shù)、線程塊的最大線程數(shù)量等。

步驟2 根據(jù)GPU設(shè)備參數(shù)及網(wǎng)格格點(diǎn)數(shù)據(jù)維數(shù)(columnNum，rowNum)，設(shè)置CUDA中線程網(wǎng)格的劃分gridDim；線程塊的劃分blockDim，這樣在CUDA中分配的總線程數(shù)：gridDim.x*gridDim.y*；blockDim.x*blockDim.y。

因?yàn)椴煌腉PU硬件的限制及網(wǎng)格數(shù)據(jù)規(guī)模等的限制，網(wǎng)格及塊劃分要根據(jù)實(shí)際環(huán)境調(diào)整。gridDim及blockDim的劃分有多種方式，本文采用二維形式的gridDim及blockDim為例。一種常用的劃分如下：gridDim(ceil(columnNum/16),ceil(rowNum/16));blockDim(16,16)，其中ceil取上限整數(shù)。

步驟3 在顯存(GPU)和內(nèi)存(CPU)中為待插值網(wǎng)格格點(diǎn)屬性數(shù)據(jù)開辟存儲(chǔ)空間，在顯存中為離散采樣點(diǎn)數(shù)據(jù)開辟存儲(chǔ)空間，將離散采樣點(diǎn)的數(shù)據(jù)由內(nèi)存拷貝到顯存。

步驟4 根據(jù)CUDA中設(shè)置的線程數(shù)和任務(wù)量進(jìn)行任務(wù)劃分，指定每個(gè)線程負(fù)責(zé)的任務(wù)。各線程計(jì)算網(wǎng)格格點(diǎn)的屬性值(搜索鄰近采樣點(diǎn)，并采用IDW插值公式計(jì)算)。各線程間的計(jì)算相對(duì)獨(dú)立，互不干擾。詳細(xì)參見基于CUDA的線程任務(wù)分解模型。

步驟5 啟動(dòng)并行計(jì)算模式(即CPU端啟動(dòng)核函數(shù))，在GPU端由GPU中的多個(gè)線程并行計(jì)算網(wǎng)格格點(diǎn)的屬性值。

步驟6 結(jié)束并行計(jì)算模式(即GPU端核函數(shù)結(jié)束，返回CPU端)，將網(wǎng)格格點(diǎn)數(shù)據(jù)由顯存(GPU)拷貝到內(nèi)存(CPU)[11,12]。

(3)輸出

輸出已插值計(jì)算的網(wǎng)格格點(diǎn)數(shù)據(jù)。

CUDA部分實(shí)現(xiàn)示例代碼：

//gridData：待插值網(wǎng)格格點(diǎn)數(shù)據(jù)，待插值網(wǎng)格行//數(shù)、列數(shù)：gridNumX，girdNumY

//離散采樣點(diǎn)數(shù)據(jù)：sampleData

//核函數(shù)

GridingParalleComputing<<>>()

//核函數(shù)類似OpenMP中并行域概念

{

//設(shè)置為二維線程網(wǎng)格，當(dāng)前線程索引CUDA內(nèi)部

//采用二維(threadIdx.x,threadIdx.y)標(biāo)識(shí)，

//所在的線程塊由(blockIdx.x,blockIdx.y)標(biāo)

//識(shí)

intI=threadIdx.x+blockIdx.x*blockDim.x;

intJ=threadIdx.y+blockIdx.y*blockDim.y;

//該線程的任務(wù)分解：不越界(i

//j

//i+=gridDim.x*blockDim.x及

//j+=gridDim.y*blockDim.y,類似OpenMP中的

//計(jì)算任務(wù)處理

for(inti=I;

i

i+=gridDim.x*blockDim.x)

for(intj=J;

j

j+=gridDim.y*blockDim.y)

//搜索與計(jì)算

SearchCompute(i,j,gridData,

gridNumX,gridNumY,

sampleData);

}

2.5 基于CUDA的線程任務(wù)分解模型

網(wǎng)格格點(diǎn)數(shù)據(jù)維數(shù)(columnNum，rowNum)，采用二維形式的gridDim及blockDim，CUDA總共分配有threads=gridDim.x*gridDim.y*blockDim.x*blockDim.y個(gè)線程。基本劃分思想是將網(wǎng)格格點(diǎn)分解到threads個(gè)組中，每個(gè)組由一個(gè)線程負(fù)責(zé)計(jì)算，重點(diǎn)是建立線程索引與網(wǎng)格格點(diǎn)索引的映射關(guān)系。在每個(gè)塊內(nèi)，CUDA會(huì)為每個(gè)線程分配唯一的索引threadIdx(threadIdx.x,threadIdx.y)，因?yàn)椴捎枚S形式的blockDim，所以threadIdx也取二維形式。需要處理columnNum列、rowNum行的網(wǎng)格格點(diǎn)的計(jì)算任務(wù)(網(wǎng)格格點(diǎn)從左下角(0,0)開始編號(hào))，通過二重循環(huán)來完成線程的任務(wù)分解：

intcol=threadIdx.x+

blockIdx.x*blockDim.x;

introw=threadIdx.y+

blockIdx.y*blockDim.y;

for(inti=col;

i

i+=gridDim.x*blockDim.x)

for(intj=threadIdx.y;

j

j+=gridDim.y*blockDim.y)

GridNodeCompute(…);

每個(gè)線程相對(duì)于所在的線程塊，有唯一索引，則在整個(gè)線程網(wǎng)格中的索引為：(col,row)=(threadIdx.x+blo-ckIdx.x*blockDim.x，threadIdx.y+blockIdx.y*blockDim.y);

即針對(duì)塊中的線程threadIdx，由該線程負(fù)責(zé)處理對(duì)應(yīng)的網(wǎng)格數(shù)據(jù)中的col列、row行索引格點(diǎn)的計(jì)算(如果col列、row行格點(diǎn)存在)。從第i=col列、第j=row行開始到columnNum-1列、rowNum-1行范圍內(nèi)的各格點(diǎn)的所有計(jì)算，其中列i的遞增間隔為gridDim.x*blockDim.x，其中行j的遞增間隔為gridDim.y*blockDim.y，原因是已分配有threads個(gè)線程，其中：x方向間隔gridDim.x*blockDim.x；y方向間隔gridDim.y*blockDim.y。網(wǎng)格數(shù)據(jù)的前0～gridDim.x*blockDim.x-1列及0～gridDim.y*blockDim.y-1行已映射到threads個(gè)線程中并行處理，所以無需重復(fù)計(jì)算，這是線程并行劃分的關(guān)鍵。i+=gridDim.x*blockDim.x與j+=gridDim.y*blockDim.y表示繼續(xù)沿x、y兩個(gè)方向搜索該線程負(fù)責(zé)計(jì)算的格點(diǎn)，主要保證計(jì)算不遺漏。i

線程任務(wù)分解模型如圖3所示，其中g(shù)ridDim(3,2),blockDim(2,2),共分配線程個(gè)數(shù)：24個(gè)，直接映射到網(wǎng)格格點(diǎn)(圖中左下方細(xì)線矩形范圍內(nèi)的24個(gè)網(wǎng)格格點(diǎn))，線程網(wǎng)格索引(0,0)(參見菱形符號(hào)標(biāo)記)，所在的塊索引block(0,0)及該塊中線程索引thread(0,0),分配的計(jì)算任務(wù)是計(jì)算網(wǎng)格格點(diǎn)(0,0)、(6，0)、(0,4)、(6,4)；線程網(wǎng)格索引(5,2)(參見三角符號(hào)標(biāo)記)，所在的塊索引block(2,1)及該塊中線程索引thread(1,0)，分配的計(jì)算任務(wù)是計(jì)算網(wǎng)格格點(diǎn)(5,2)、(5，6)。

圖3 基于CUDA的線程任務(wù)分解模型

線程任務(wù)分解模型的計(jì)算均衡(負(fù)載平衡)：一般情況下，首先GPU分配的所有線程參與計(jì)算，對(duì)計(jì)算任務(wù)進(jìn)行了基于全部線程的分解；其次，通過任務(wù)分解模型示意圖可看出，各線程處理的格點(diǎn)在二維網(wǎng)格上均勻間隔分布，整體來說，計(jì)算比較均衡，即線程任務(wù)分解模型具有較好的負(fù)載均衡。

3 算法實(shí)驗(yàn)

采用傳統(tǒng)的串行算法(CPU單核)、并行算法(CPU多核OpenMP)、并行算法(GPU多核CUDA)3種方式，統(tǒng)計(jì)總消耗時(shí)間、加速比(串行算法時(shí)間/并行算法時(shí)間之比)等。

實(shí)驗(yàn)數(shù)據(jù)選自全球氣象觀測(cè)709個(gè)站點(diǎn)，地理范圍：經(jīng)度-180°～180°，維度-90°～90°，數(shù)據(jù)網(wǎng)格以360×180為基準(zhǔn)(即經(jīng)緯度按1°劃分)，采用網(wǎng)格系數(shù)表達(dá)網(wǎng)格疏密程度，例如網(wǎng)格系數(shù)為2，表示360×2×180×2。

3.1 實(shí)驗(yàn)一、OpenMP并行優(yōu)化及分析

(1)實(shí)驗(yàn)環(huán)境：國(guó)產(chǎn)麒麟Kylin操作系統(tǒng)(32位)，CPU：申威(SW410)、龍芯(Loongson3A2000)、英特爾(Intel I5)，均為四核。

(2)實(shí)驗(yàn)結(jié)果：統(tǒng)計(jì)消耗時(shí)間及加速比，并繪制成曲線圖，結(jié)果參見表1、表2、圖4～圖7。

表1 網(wǎng)格系數(shù)2時(shí)多核并行加速性能測(cè)試統(tǒng)計(jì)

表2 網(wǎng)格系數(shù)4時(shí)多核并行加速性能測(cè)試統(tǒng)計(jì)

圖4 網(wǎng)格系數(shù)2時(shí)CPU消耗時(shí)間對(duì)比

圖5 網(wǎng)格系數(shù)2時(shí)CPU加速比曲線

圖6 網(wǎng)格系數(shù)4時(shí)CPU消耗時(shí)間對(duì)比

圖7 網(wǎng)格系數(shù)4時(shí)CPU加速比曲線

(3)實(shí)驗(yàn)分析：選取3種CPU及兩種網(wǎng)格劃分，對(duì)IDW算法并行實(shí)驗(yàn)，從表格及曲線圖可看出：

加速比：經(jīng)過并行加速后，加速比達(dá)到3.75～5.614。龍芯(最大加速比3.82)及申威CPU(最大加速比3.75)接近4倍，這二者基本相當(dāng)，英特爾CPU超過4倍(最大加速比5.614)，加速比最高。

總消耗時(shí)間：經(jīng)過并行加速后，消耗時(shí)間大量減少。龍芯與申威二者基本相當(dāng)，英特爾消耗最少：6.883。

基本結(jié)論：采用OpenMP并行優(yōu)化IDW算法是有效的，整體加速比基本符合與CPU的核數(shù)相當(dāng)這一經(jīng)驗(yàn)值，也說明了OpenMP加速受CPU核心數(shù)的限制。從加速比可看出，國(guó)產(chǎn)與商用CPU存在差距。

3.2 實(shí)驗(yàn)二、CUDA并行優(yōu)化及分析

(1)實(shí)驗(yàn)環(huán)境：CPU：Intel I2，雙核；GPU：NVIDIA Quadro 240；操作系統(tǒng)：windows XP 32位；CUDA Cores：16個(gè)；CUDA v6.5。

(2)實(shí)驗(yàn)結(jié)果：統(tǒng)計(jì)消耗時(shí)間及加速比，并繪制成曲線圖，結(jié)果參見表3、圖8、圖9。

(3)實(shí)驗(yàn)分析：通過采用不同疏密網(wǎng)格劃分，在Windows操作系統(tǒng)下對(duì)IDW算法并行優(yōu)化實(shí)驗(yàn)，從表格及曲線圖可看出：

加速比：加速比在1.9～5.9之間，性能提升顯著。隨著網(wǎng)格數(shù)據(jù)量的增加，加速比更高，可看出GPU比較適合大規(guī)模的數(shù)據(jù)并行計(jì)算。

表3 GPU并行加速性能測(cè)試統(tǒng)計(jì)

圖8 GPU加速前后消耗時(shí)間對(duì)比

圖9 GPU加速比曲線

總消耗時(shí)間：經(jīng)過GPU并行加速后，優(yōu)化后消耗時(shí)間比優(yōu)化前消耗的時(shí)間大量減少。針對(duì)網(wǎng)格系數(shù)4，比較表2與表3消耗時(shí)間(Intel I5：34.352， GPU：20.187)，可看出CUDA比OpenMP更有效。

基本結(jié)論：采用CUDA并行優(yōu)化IDW算法更加有效。

3.3 結(jié)論及應(yīng)用效果

基于OpenMP與CUDA的并行優(yōu)化，具有較好的性能提升。將本文并行優(yōu)化技術(shù)應(yīng)用到國(guó)產(chǎn)下某型號(hào)氣象業(yè)務(wù)系統(tǒng)項(xiàng)目中，該應(yīng)用系統(tǒng)性能得到較大的改善，效果如圖10所示。

圖10 局部區(qū)域等值線顯示效果

4 結(jié)束語

針對(duì)傳統(tǒng)的IDW算法采用并行優(yōu)化處理后，實(shí)現(xiàn)了較好的性能提升，體現(xiàn)了CPUGPU多核的硬件加速能力，驗(yàn)證了軟硬件結(jié)合并行優(yōu)化思路的正確性，線程任務(wù)分解模型也極大地減輕了并行編程任務(wù)劃分的難度，為國(guó)產(chǎn)軟件的性能提升提供了一定的指導(dǎo)。下一步可研究：①目前CUDA僅支持NVIDIA顯卡，支持多種顯卡的OpenCL通用計(jì)算技術(shù)；②基于Hadoop MapReduce、Spark等分布式大規(guī)模數(shù)據(jù)的并行計(jì)算技術(shù)；③針對(duì)OpenMP及CUDA的性能優(yōu)化。

參考文獻(xiàn)：

[1]ZHANG Xiaoqing,GONG Bo,TIAN Liyun,et al.Study on the performance optimization of application running on the platform built with Chinese owned technology[J].Software,2015,36(2):5-9(in Chinese).[張曉清,龔波,田麗韞,等.國(guó)產(chǎn)自主可控應(yīng)用性能優(yōu)化研究[J].軟件，2015,36(2):5-9.]

[2]XIE Zhouyu,ZANG Fei.Performance optimization technology for information system based on domestic software and hardware[J].Command Information System and Technology,2014,5(3):59-63(in Chinese).[謝宙宇,臧飛.基于國(guó)產(chǎn)軟硬件的信息系統(tǒng)性能優(yōu)化技術(shù)[J].指揮信息系統(tǒng)與技術(shù),2014,5(3):59-63.]

[3]FU He,XIE Yongfang.From design of a cockpit display system based on homemade GPU JM5400[J].Computer Engineering & Science,2016,38(10):2083-2090(in Chinese).[符鶴,謝永芳.基于國(guó)產(chǎn)化圖形芯片JM5400的座艙顯示系統(tǒng)設(shè)計(jì)[J].計(jì)算機(jī)工程與科學(xué),2016,38(10):2083-2090.]

[4]OpenMPI[EB/OL].[2017-04-07].http://www.open-mpi.org.

[5]LIU Wenzhi.Parallel computing and performance optimization[M].Beijing:China Machine Press,2015:2-174(in Chinese).[劉文志.并行算法設(shè)計(jì)與性能優(yōu)化[M].北京:機(jī)械工業(yè)出版社,2015:2-174.]

[6]CHEN Huan,XIE Jian.Kriging interpolation algorithm based on OpenMP[J].Computer Science,2012,39(6A):392-394(in Chinese).[陳歡,謝健.基于OpenMP的Kriging插值算法研究[J].計(jì)算機(jī)科學(xué),2012,39(6A):392-394.]

[7]LING Bin.DENG Yan.YU Shibo.Processing for accelerated sparse PCNN moving target detection algorithm with CUDA[J].Computer Engineering and Design,2016,37(12):3301-3305(in Chinese).[凌濱,鄧艷,于士博.CUDA并行加速的稀疏PCNN運(yùn)動(dòng)目標(biāo)檢測(cè)算法[J].計(jì)算機(jī)工程與設(shè)計(jì),2016,37(12):3301-3305.]

[8]GU Yuhang,ZHAO Wei,LI Li,et al.Design and realization of parallel topology algorithm of vector spatial data based on OpenMP[J].Engineering of Surveying and Mapping,2015,24(11):22-27(in Chinese).[谷宇航,趙偉,李力,等.基于OpenMP的矢量空間數(shù)據(jù)并行拓?fù)渌惴ㄔO(shè)計(jì)與實(shí)現(xiàn)[J].測(cè)繪工程,2015,24(11):22-27.]

[9]WU Lingda,HAO Liyun,FENG Xiaomeng,et al.Combining isosurface rendering and volume rendering for method of electromagnetic environment visualization[J].Journal of Beijing University of Aeronautics and Astronautics,2017,43(5):887-893(in Chinese).[吳玲達(dá),郝利云,馮曉萌,等.結(jié)合等值面繪制與體繪制的電磁環(huán)境可視化方法[J].北京航空航天大學(xué)學(xué)報(bào),2017,43(5):887-893.]

[10]JIA Haipeng.Research of parallel optimization technologies on GPU computing platforms[D].Qingdao:Ocean University of China,2012:1-105(in Chinese).[賈海鵬.面向GPU計(jì)算平臺(tái)的若干并行優(yōu)化關(guān)鍵技術(shù)研究[D].青島:中國(guó)海洋大學(xué),2012:1-105.]

[11]DAI Chen,CHEN Peng,YANG Donglei,et al.On parallel programming and optimization for multi-core[J].Computer Application and Software,2013,30(12):198-202(in Chinese).[戴晨,陳鵬,楊冬蕾,等.面向多核的并行編程和優(yōu)化研究[J].計(jì)算機(jī)應(yīng)用與軟件,2013,30(12):198-202.]

[12]SHAO Tianjiao.Research on GPU-based parallel computation performance optimization[D].Changchun:Jilin University,2014:1-48(in Chinese).[邵天驕.基于GPU并行計(jì)算的性能優(yōu)化研究[D].長(zhǎng)春:吉林大學(xué),2014:1-48.]

[13]ZHOU Hongwen,DU Yuxing,REN Gaofeng,et al.Analysis of ore reserves based on Arc-Engine and IDW[J].Journal of Wuhan University of Technology,2014,36(2):115-119(in Chinese).[周洪文,杜玉星,任高峰,等.基于Arc-Engine與IDW法的礦石儲(chǔ)量分析研究[J].武漢理工大學(xué)學(xué)報(bào),2014,36(2):115-119.]

[14]DUAN Ping,SHENG Yehua,LI Jia,et al.Adaptive IDW interpolation method and its application in the temperature field[J].Geographical Research,2014,33(8):1417-1426(in Chinese).[段平,盛業(yè)華,李佳,等.自適應(yīng)的IDW插值方法及其在氣溫場(chǎng)中的應(yīng)用[J].地理研究,2014,33(8):1417-1426.]

[15]LI Zhengquan,WU Yaoxiang.Inverse distance weighted interpolation involving position shading[J].Acta Geodaetica et Cartographica Sinica,2015,44(1):91-98 (in Chinese).[李正泉,吳堯祥.顧及方向遮蔽性的反距離權(quán)重插值法[J].測(cè)繪學(xué)報(bào),2015,44(1):91-98.]