CPU-GPU協(xié)同加速K riging插值的負載均衡方法*

姜春雷，張樹清(.中國科學院東北地理與農(nóng)業(yè)生態(tài)研究所，吉林長春3002；2.中國科學院大學，北京00049)

CPU-GPU協(xié)同加速K riging插值的負載均衡方法*

姜春雷1，2，張樹清1
(1.中國科學院東北地理與農(nóng)業(yè)生態(tài)研究所，吉林長春130102；2.中國科學院大學，北京100049)

Kriging插值算法被廣泛應用于地學各領域，有著極其重要的現(xiàn)實意義，但在面對大規(guī)模輸出網(wǎng)格及大量輸入采樣點時，不可避免地遇到了性能瓶頸。利用OpenCL和OpenMP在異構平臺上實現(xiàn)了CPU與GPU協(xié)同加速普通Kriging插值。針對Kriging插值中采樣點的不規(guī)則分布及CPU和GPU由于體系結構差異對其的不同適應性，提出一種基于不同設備間計算性能的差異和數(shù)據(jù)分布特點的負載均衡方法。試驗結果表明，該方法能有效提高普通Kriging插值速度，同時還能節(jié)約存儲空間和提高訪存效率。

通用計算圖形處理器；開放運算語言；Kriging插值；負載均衡

Kriging插值算法是一種考慮區(qū)域變量空間變異結構(自相關結構)特點的線性無偏最優(yōu)估計插值方法［1］，它在地學的很多領域都有著重要的應用。但是，當Kriging插值的輸出網(wǎng)格增大及輸入采樣點增多時，其計算時間急劇增加［2］。當前，集成Kriging插值算法的主流地理資訊系統(tǒng)(Geographic Information System，GIS)仍以串行計算為基礎框架，不能充分利用多核處理器的優(yōu)勢［3］。隨著各種并行技術的不斷成熟，地理空間柵格數(shù)據(jù)算法的并行化研究成為新的熱點［4］。很多研究者試圖通過不同的并行技術來加速Kriging插值。目前，針對Kriging插值的加速方法主要基于的編程模型有消息傳遞接口(Message Passing Interface，MPI)［5］，OpenMP［6］和通用計算圖形處理器(General Purpose Graphics Processing Units，GPGPU)［7-8］?；谶@些不同編程模型的Kriging插值算法都取得了很好的加速效果，尤其是GPGPU編程模型的加速比更高。然而，GPGPU加速算法僅將CPU作為一個流程控制端，造成了GPU計算時CPU空閑等待，顯然，CPU與GPU共同承擔計算任務的方式是未來協(xié)同并行計算的發(fā)展方向［9］。

多設備協(xié)同加速需要解決的核心問題之一是設備間的負載均衡問題。目前關于計算設備間負載均衡的研究主要集中在根據(jù)不同設備的性能差異分配計算的靜態(tài)負載均衡［10-11］以及在計算過程中根據(jù)不同計算設備進度分配計算的動態(tài)負載均衡［12］。夏飛等［13］對于異構設備(例如CPU與GPU)間的負載均衡給出了根據(jù)計算能力分發(fā)計算任務的靜態(tài)負載均衡方法；趙斯思和周成虎［14］對于單個GPU上多個kernel同時執(zhí)行，提出一種基于動態(tài)規(guī)劃的動態(tài)負載均衡方法，有效地加速了算法的執(zhí)行。對于Kriging插值而言，由于受限于自然環(huán)境、經(jīng)濟和人力等條件的限制，采樣點通常不規(guī)則，表現(xiàn)為分布不均勻。由于體系結構的差異，對于數(shù)據(jù)量相同但分布密度不同的數(shù)據(jù)，不同體系結構的硬件通常表現(xiàn)出不同的計算速度。這表明負載均衡時不僅要考慮不同設備間性能的差異，還應考慮數(shù)據(jù)的特征。

本文提出了一種新的綜合考慮設備計算性能和數(shù)據(jù)特征的負載均衡方法(Load Balancing based on Computation Performance and Data Distribution，LBCPDD)，來實現(xiàn)設備間高效合理的負載均衡，以達到更快加速Kriging插值的目的。

1 背景知識

1.1 普通K riging

普通Kriging插值公式可以表示為:

a.無偏性條件。

其中，E表示數(shù)學期望，m是中間變量。

b.最優(yōu)性條件(即估計值誤差的方差最小)。

依據(jù)拉格朗日乘數(shù)法原理，建立拉格朗日函數(shù):

r是變異函數(shù)，上述方程組可用矩陣表示如下:解上述方程組，求得λi代入普通Kriging插值公式可計算出待估未知點的值。

1.2 通用計算圖形處理器

GPGPU將原本只用于圖形處理的硬件用于通用計算領域，OpenCL和統(tǒng)一計算設備架構(Compute Unified Device Architecture，CUDA)是實現(xiàn)這個理念的兩種主要軟件開發(fā)平臺，OpenCL由于能夠在不同的硬件平臺上運行，更具應用前景。目前，GPGPU已經(jīng)被廣泛地應用于各領域的科學研究，其相關知識已經(jīng)被相關文獻大量描述，本文僅介紹在其他文獻中被較少描述且被應用于本研究的相關內(nèi)容。AMD公司和NVIDIA公司對同一硬件概念有不同的稱謂，由于本研究基于AMD公司的GPU實現(xiàn)，因此在以后的描述中以AMD公司的概念為準。在AMD公司的GPU中，一個計算單元(Compute Unit，CU)可以同時運行多個workgroup，具體硬件運行的workgroup數(shù)及workgroup內(nèi)的線程數(shù)由GPU的硬件環(huán)境決定，主要是內(nèi)存數(shù)量，更少的內(nèi)存消耗，可以讓每個CU容納更多的workgroup和每個workgroup容納更多線程。workgroup內(nèi)線程被分成若干個wavefront(等同于NVIDIA硬件中的warp)，它是硬件調(diào)度線程執(zhí)行的最小單位，其內(nèi)部的線程以嚴格鎖同步方式執(zhí)行。如果wavefront內(nèi)線程存在不同分支，則每個線程會因為鎖同步方式而執(zhí)行全部分支(無計算任務線程執(zhí)行空操作)，導致執(zhí)行時間增加。數(shù)量足夠的wavefront可以保證在一個wavefront等待內(nèi)存訪問結果時，其他的wavefront可以被調(diào)入執(zhí)行單元執(zhí)行(即交叉訪問)，這種方式大大提高了GPU的利用效率。因此wavefront對程序的執(zhí)行速度有著巨大的影響。OpenCL并未對開發(fā)者提供wavefront抽象，但是通過對不同線程處理數(shù)據(jù)的合理安排可以提高GPU的運行效率，即提高wavefront數(shù)量以及減少單個wavefront內(nèi)的線程分支。

2 基于計算性能與數(shù)據(jù)分布的負載均衡方法

在Kriging插值過程中，相鄰的未知點通常有相同的鄰近點。由于變異函數(shù)矩陣是基于待估未知點的鄰近點構造的，因此相鄰的未知點通常具有相同的變異函數(shù)矩陣。根據(jù)這一特點，在Kriging插值算法中增加一個步驟用于檢查相鄰未知點間是否具有相同的鄰近點:若存在，則只存儲其中一個未知點的鄰近點，并且僅對其構造變異函數(shù)矩陣和求逆矩陣，即將相同鄰近矩陣合并存儲及計算，其他的未知點在計算過程中使用這個計算的結果。

采樣點一般是不規(guī)則(不均勻)分布的，局部Kriging插值在對每個未知點進行插值時，需要求出離它最近的若干個采樣點，在采樣點密度較小的情況下，相鄰未知點具有相同鄰近點的概率更大，如圖1左上角所示，兩個未知點(*)具有相同的臨近采樣點(·)；而采樣點較密集時則相反(如圖1右下角所示):兩個未知點(*)具有不同的采樣點(·)。對于GPU，較多的具有相同鄰近點的線程使wavefront內(nèi)產(chǎn)生更少的分支；由于合并存儲和計算，采樣點稀疏區(qū)域的未知點預測會消耗更少的內(nèi)存和計算；內(nèi)存消耗的減少讓GPU能夠容納更多的線程并發(fā)執(zhí)行以便GPU線程的交叉執(zhí)行減少內(nèi)存延遲，即不同采樣點密度會影響GPU的性能。因而在負載均衡過程中不僅需考慮CPU與GPU自身的計算能力差異，還應考慮CPU與GPU對不同數(shù)據(jù)空間分布特征的適應性差異。在算法中，將插值區(qū)域劃分為若干個網(wǎng)格，統(tǒng)計每一個網(wǎng)格里面的采樣點個數(shù)，包含較多采樣點的網(wǎng)格中的未知點預測交給CPU處理，反之，交給GPU處理(如圖2所示，灰色表示密度較大，白色較小)。

圖1 鄰近點查找與采樣點分布的關系Fig.1 Relation between searching nearest neighbor and the distribution of samples

CPU與GPU之間任務的分配比例通常依據(jù)設備浮點計算能力，然而，考慮到Kriging插值的復雜性及數(shù)據(jù)傳輸?shù)挠绊?，浮點計算能力不能體現(xiàn)出CPU與GPU在計算Kriging插值時的性能差異，因此，依據(jù)CPU與GPU單獨計算Kriging插值時的時間初步估計GPU任務分配比例:

圖2 CPU與GPU負載分配Fig.2 Load allocation between CPU and GPU

然后根據(jù)結果向GPU適當傾斜任務量。上述任務分配估計方法體現(xiàn)了CPU和GPU對本算法的性能差異，因此僅需CPU或GPU單獨運行一次就可以估算出結果，在之后的計算中即便針對不同的數(shù)據(jù)集也不需要再重新估算任務分配比例。由于CPU端的串行執(zhí)行并不能完全利用CPU的計算能力，因此式(5)中TimeCPU指的是OpenMP并行方式下的執(zhí)行時間。

在LBCPDD算法的實現(xiàn)中，CPU端以OpenMP方式并行。為了對提出的GPU算法的性能進行比較，實現(xiàn)了同樣采用局部插值的Shi和Ye描述的Kriging算法［8］。

3 試驗與分析

采用的GPU為華碩R9 280，顯示核心為AMD Tahiti XTL次世代圖形核心(Graphics Core Next，GCN)架構，擁有3G GDDR5顯存，流處理器數(shù)量為1792個；CPU為AMD FX-8150主頻為3.6GHz，8核心，8線程；操作系統(tǒng)為Win7 64位。實驗數(shù)據(jù)為從美國國家海洋和大氣局(National Oceanic and Atmospheric Administration，NOAA)下載的全球降雨數(shù)據(jù)，為了保證魯棒性，根據(jù)實驗所需采樣點數(shù)據(jù)集的大小，隨機選擇相應數(shù)量的點作為采樣點。所有的實驗中插值的輸出網(wǎng)格大小均為1280×800個像素。

試驗1：CPU性能。為了測試算法在不同數(shù)據(jù)量下的表現(xiàn)，選取了5組數(shù)據(jù)集進行測試，數(shù)據(jù)量每次增加1000個采樣點；OpenMP分別選擇2核、4核和8核進行測試(在表1中表示為OMP2，OMP4和OMP8)。算法執(zhí)行時間見表1:隨著核數(shù)的增加，運行時間持續(xù)減少；隨著數(shù)據(jù)量增加，不同算法的執(zhí)行時間都有所增加。圖3給出了OpenMP在不同線程(核心)下的加速比，從圖3可以看出隨著線程數(shù)的增加，OpenMP的加速比也隨之增加，最高達到了3.9倍左右；不同線程的OpenMP方法加速比對數(shù)據(jù)集大小不敏感，只是數(shù)據(jù)量比較大時有些輕微波動。

表1 串行與OpenMP運行時間對比Tab.1 Time comparison of sequential code and OpenMP ms

圖3 OpenMP加速比比較Fig.3 Comparison of speed-up of OpenMP

試驗2：GPU性能。表2給出了GPU并行的運行時間及加速比。GPU1并行程序依據(jù)Shi和Ye描述的方法實現(xiàn)［8］，GPU2是改進后的GPU算法，增加了檢查相鄰未知點間是否具有相同的鄰近點的步驟。從表2可以看出:GPU算法在不同的數(shù)據(jù)集下都表現(xiàn)出很好的加速效果，GPU2最高達到21.53倍；GPU2由于對Kriging插值過程中各未知點的鄰近點進行了比較，對擁有相同鄰近點的未知點，只計算其中的一個未知點的變異函數(shù)矩陣和逆矩陣，因此GPU2的加速比明顯高于GPU1的加速比；同OpenMP算法類似，GPU算法的運行時間隨數(shù)據(jù)集的大小增加而增加，但其加速比同數(shù)據(jù)集的大小無關，波動不大。

試驗3：CPU-GPU性能。表3給出了通過LBCPDD進行負載均衡后的CPU-GPU協(xié)同執(zhí)行時間。表中，OpenMP表示CPU端并行方式，采用8核并行，任務量占1/8；GPU2L表示GPU端使用表2中GPU2算法且通過LBCPDD方法負載均衡分配數(shù)據(jù)，任務量占7/8。由于GPU并行與OpenMP并行同時執(zhí)行，因此算法的執(zhí)行時間取決于較長的時間(負載分配時間小于3ms，未列入統(tǒng)計)。表3中性能提升指總時間相對GPU單獨執(zhí)行時間(表2中GPU2)的性能提升((GPU2-總時間)/GPU2)。從表3中可以看出，CPU-GPU的運算總時間比GPU單獨運算時在不同的數(shù)據(jù)集下都有提高，最高達到了29.95%。CPU-GPU算法減少的時間一部分來自于CPU分擔了部分計算，另外一部分來自于GPU由于處理數(shù)據(jù)的優(yōu)化而節(jié)約的時間。

表2 GPU程序運行時間及加速比Tab.2 Time and speed-up of GPU program

表3 引入LBCPDD后CPU-GPU協(xié)同運行時間Tab.3 Run time of CPU-GPU cooperation with LBCPDD

試驗4：LBCPDD方法對性能的影響。表4給出了OpenMP及GPU并行在應用LBCPDD前后完成各自分配任務量的計算時間及性能變化。標準方法指在分配任務時根據(jù)插值點在二維平面上的坐標順序?qū)⑶?/8分配給GPU，剩余部分分配給CPU。表中GPU2S和GPU2L分別表示經(jīng)過標準方法和LBCPDD方法進行任務分配后的GPU端執(zhí)行，任務量為7/8；OpenMP1和OpenMP2分別表示經(jīng)過標準方法和LBCPDD方法進行任務分配后的CPU端執(zhí)行，任務量為1/8；GPU2L性能提升指GPU2L相對于GPU2S的性能提升((GPU2S-GPU2L)/GPU2S)；OpenMP2性能提升指OpenMP2相對OpenMP1的性能提升((OpenMP1-OpenMP2)/OpenMP1)。從表4中可以看出，LBCPDD方法有效減少了GPU的計算時間，性能最多提高50.43%，而對OpenMP方法的性能影響相對較小。這主要是由于將適合GPU計算的數(shù)據(jù)分給了GPU，而將其余數(shù)據(jù)分給了對數(shù)據(jù)分布并不敏感的CPU。由于OpenMP運行時，各線程間以無序方式異步運行，這樣做有利于OpenMP線程間動態(tài)負載均衡，但如果各線程間需要同步，則會極大地影響其速度，故OpenMP程序無法從減少冗余矩陣和數(shù)據(jù)重組中獲得好處。

表4 LBCPDD對CPU和GPU的影響Tab.4 Impact of LBCPDD on GPU and CPU

在CPU與GPU協(xié)同計算中，對采樣點的密度進行了分析，將采樣點分布密度較低的區(qū)域分給了GPU。這樣做有很多好處:更多的冗余矩陣分布在一個workgroup中，將有利于單個wavefront內(nèi)的線程具有更少甚至沒有分支。大量的冗余矩陣分布在一個workgroup內(nèi)會使后續(xù)計算讀取同一塊內(nèi)存區(qū)域，從而可以基于GPU內(nèi)存存取機制(廣播)減少延遲。CPU與GPU之間解決內(nèi)存延遲使用兩種不同策略，CPU用了大量的芯片空間構造片內(nèi)緩存，主要通過多級緩存來減少內(nèi)存訪問的延遲；而GPU芯片大量空間被用來構造計算單元，通過大量線程的交叉訪問來隱藏內(nèi)存延遲，更多的冗余矩陣導致更少的內(nèi)存消耗，從而可以讓更多的線程調(diào)度到CU中執(zhí)行，為線程的交叉訪問提供條件。在冗余矩陣合并存儲時，通過一個全局變量記錄每個線程對應的矩陣編號，由于OpenCL不支持全局同步，因此，通過對全局變量的原子操作來完成編號的記錄，原子操作具有排他性，會嚴重影響性能，冗余矩陣減少后，原子操作減少。一個workgroup內(nèi)冗余矩陣增多，意味著線程計算時的內(nèi)存消耗減少(大量線程空操作)，這樣會減少變量溢出到全局存儲器中，全局存儲器的存儲周期約500個時鐘周期，由于等待時間還可能翻倍，而寄存器存儲周期只有1個時鐘周期，馬安國等［15］在通過預取對GPU程序進行優(yōu)化時，性能影響最高達38%?；谝陨显?，數(shù)據(jù)重分配后GPU計算的時間明顯減少。從表4可以看到，隨著采樣點數(shù)據(jù)的增多，采樣點的密度增加，導致冗余矩陣減少，這將導致性能改善降低，表4中數(shù)據(jù)重分配與按序分配時GPU性能改善隨采樣點數(shù)據(jù)增加而減少就驗證了這一點。

4 結論

Kriging插值算法是一種廣泛應用的算法。利用CPU與GPU協(xié)同計算對Kriging插值進行加速，基于CPU與GPU的計算性能差異和插值數(shù)據(jù)采樣點的不均勻分布，對CPU與GPU進行負載均衡。測試結果表明，采用了LBCPDD后的CPU-GPU協(xié)同加速方法優(yōu)于單純依賴CPU的OpenMP加速方法、完全GPU的算法和未采用LBCPDD的CPU-GPU協(xié)同加速方法。LBCPDD法對其他數(shù)據(jù)分布敏感的CPU-GPU協(xié)同加速算法在進行負載均衡時也同樣具有參考價值。

References)

［1］Krige D G.A statistical approach to some basic mine valuation problems on the witwatersrand［J］.Journal of the Chemical Metallurgical＆Mining Society of South Africa，1951，94(3):95-111.

［2］De Baar JH S，Dwight R P，Bijl H.Speeding up Kriging through fast estimation of the hyperparameters in the frequency-domain［J］.Computers＆Geosciences，2013，54:99-106.

［3］吳立新，楊宜舟，秦成志，等.面向新型硬件構架的新一代GIS基礎并行算法研究［J］.地理與地理信息科學，2013，29(4):1-8.WU Lixin，YANG Yizhou，QIN Chengzhi，et al.On basic geographic parallel algorithms of new generation GIS for new hardware architectures［J］.Geography and Geo-Information Science，2013，29(4):1-8.(in Chinese)

［4］程果，陳犖，吳秋云，等.一種面向復雜地理空間柵格數(shù)據(jù)處理算法并行化的任務調(diào)度方法［J］.國防科技大學學報，2012，34(6):61-65.CHENG Guo，CHEN Luo，WU Qiuyun，et al.A task schedulingmethod for parallelization of complicated geospatial raster data processing algorithms［J］.Journal of National University of Defense Technology，2012，34(6):61-65.(in Chinese)

［5］Hu H D，Shu H.An improved coarse-grained parallel algorithm for computational acceleration of ordinary Kriging interpolation［J］.Computers＆Geosciences，2015，78: 44-52.

［6］Cheng T P，Li D D，Wang Q.On parallelizing universal Kriging interpolation based on OpenMP［C］//Proceedings of the9th International Symposium on Distributed Computing and Applications to Business，Engineering and Science，Hong Kong:IEEE，2010:36-39.

［7］Cheng T P.Accelerating universal Kriging interpolation algorithm using CUDA-enabled GPU［J］.Computers＆Geosciences，2013，54:178-183.［8］Shi X，Ye F.Kriging interpolation over heterogeneous computer architectures and systems［J］.GIScience＆Remote Sensing，2013，50(2):196-211.

［9］盧風順，宋君強，銀?？担?CPU/GPU協(xié)同并行計算研

究綜述［J］.計算機科學，2011，38(3):5-9.

LU Fengshun，SONG Junqiang，YIN Fukang，et al.Survey of CPU/GPU synergetic parallel computing［J］.Computer Science，2011，38(3):5-9.(in Chinese)

［10］方留楊，王密，李德仁，等.負載分配的CPU/GPU高分

辨率衛(wèi)星影像調(diào)制傳遞補償方法［J］.測繪學報，2014，43(6):598-606.FANG Liuyang，WANG Mi，LIDeren，et al.A workloaddistribution based CPU/GPU MTF compensation approach for high resolution satellite images［J］.Acta Geodaetica et Cartographica Sinica，2014，43(6):598-606.(in Chinese)［11］Wang S，Armstrong M P.A theoretical approach to the use of cyberinfrastructure in geographical analysis［J］.International Journal of Geographical Information Science，2009，23(2): 169-193.

［12］Dong B，Li X，Wu QM，et al.A dynamic and adaptive load balancing strategy for parallel file system with large-scale I/O servers［J］.Journal of Parallel and Distributed Computing，2012，72(10):1254-1268.

［13］夏飛，朱強華，金國慶.基于CPU-GPU混合計算平臺的RNA二級結構預測算法并行化研究［J］.國防科技大學學報，2013，35(6):138-146.XIA Fei，ZHU Qianghua，JIN Guoqing.Accelerating RNA secondary structure prediction applications based on CPUGPU hybrid platforms［J］.Journal of National University of Defense Technology，2013，35(6):138-146.(in Chinese)

［14］趙斯思，周成虎.GPU加速的多邊形疊加分析［J］.地理科學進展，2013，32(1):114-120.ZHAO Sisi，ZHOU Chenghu.Accelerating polygon overlay analysis by GPU［J］.Progress in Geography，2013，32(1): 114-120.(in Chinese)

［15］馬安國，成玉，唐遇星，等.GPU異構系統(tǒng)中的存儲層次和負載均衡策略研究［J］.國防科技大學學報，2009，31(5):38-43.MA Anguo，CHENG Yu，TANG Yuxing，et al.Research on memory hierarchy and load balance strategy in heterogeneous system based on GPU［J］.Journal of National University of Defense Technology，2009，31(5):38-43.(in Chinese)

A load balancing method in accelerating K riging algorithm on CPU-GPU heterogeneous p latforms

JIANG Chunlei1，2，ZHANG Shuqing1
(1.Northeast Institute of Geography and Agroecology，Chinese Academy of Sciences，Changchun 130102，China；2.University of Chinese Academy of Sciences，Beijing100049，China)

Kriging interpolation algorithm is of great practical significance and is widely applied to various fields of geoscience.However，Kriging interpolation would inevitably encounter the performance bottleneck when the output grid or input samples increase.Implemented with OpenCL and OpenMP，the ordinary Kriging interpolation was accelerated on heterogeneous platforms:GPU and CPU.By considering the performance difference of CPU and GPU on the densities of samples，a new load balancing method of LBCPDD(Load Balancing based on Computation Performance and Data Distribution)was proposed，in which not only hardware performance but also data distribution characteristics were taken into account.Experiment results show that LBCPDDmethod can effectively enhance the speed of ordinary Kriging，savememory space and improve the efficiency ofmemory access.

general purpose graphics processor units；open computing language；Kriging interpolation；load balancing

F209

A

1001-2486(2015)05-035-05

10.11887/j.cn.201505006

http://journal.nudt.edu.cn

2015-06-24

國家自然科學基金資助項目(41271196)；中國科學院重點部署資助項目(KZZD-EW-07-02)

姜春雷(1981—)，男，吉林德惠人，博士研究生，E-mail:jiangchunlei@iga.ac.cn；張樹清(通信作者)，男，研究員，博士，博士生導師，E-mail:zhangshuqing@iga.ac.cn