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

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

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

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

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

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

Kriging插值算法是一種考慮區域變量空間變異結構(自相關結構)特點的線性無偏最優估計插值方法［1］，它在地學的很多領域都有著重要的應用。但是，當Kriging插值的輸出網格增大及輸入采樣點增多時，其計算時間急劇增加［2］。當前，集成Kriging插值算法的主流地理資訊系統(Geographic Information System，GIS)仍以串行計算為基礎框架，不能充分利用多核處理器的優勢［3］。隨著各種并行技術的不斷成熟，地理空間柵格數據算法的并行化研究成為新的熱點［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共同承擔計算任務的方式是未來協同并行計算的發展方向［9］。

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

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

1 背景知識

1.1 普通K riging

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

a.無偏性條件。

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

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

依據拉格朗日乘數法原理，建立拉格朗日函數:

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

1.2 通用計算圖形處理器

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

2 基于計算性能與數據分布的負載均衡方法

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

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

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

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

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

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

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

3 試驗與分析

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

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

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

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

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

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

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

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

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

4 結論

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

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］程果，陳犖，吳秋云，等.一種面向復雜地理空間柵格數據處理算法并行化的任務調度方法［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協同并行計算研

究綜述［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高分

辨率衛星影像調制傳遞補償方法［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異構系統中的存儲層次和負載均衡策略研究［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