基于CUDA的熱傳導GPU并行算法研究

孟小華，黃叢珊，朱麗莎

(暨南大學 a. 計算機科學系；b. 天體測量、動力學與空間科學中法聯合實驗室，廣州 5 10632)

基于CUDA的熱傳導GPU并行算法研究

孟小華a,b，黃叢珊a，朱麗莎a,b

(暨南大學 a. 計算機科學系；b. 天體測量、動力學與空間科學中法聯合實驗室，廣州 5 10632)

在熱傳導算法中，使用傳統的CPU串行算法或MPI并行算法處理大批量粒子時，存在執行效率低、處理時間長的問題。而圖形處理單元(GPU)具有大數據量并行運算的優勢，為此，在統一計算設備架構(CUDA)并行編程環境下，采用CPU和GPU協同合作的模式，提出并實現一個基于CUDA的熱傳導GPU并行算法。根據GPU硬件配置設定Block和Grid的大小，將粒子劃分為若干個block，粒子輸入到GPU顯卡中并行計算，每一個線程執行一個粒子計算，并將結果傳回CPU主存，由CPU計算出每個粒子的平均熱流。實驗結果表明，與CPU串行算法在時間效率方面進行對比，該算法在粒子數到達16 000時，加速比提高近900倍，并且加速比隨著粒子數的增加而加速提高。

熱傳導算法；圖形處理單元；統一計算設備架構；并行計算；時間效率；加速比

1 概述

熱傳導算法可以運用于很多領域，如建筑工程對各種材料的研究和選擇，以及地質勘探研究等方面，幾乎所有的材料成形模擬過程都涉及熱傳導，但是國內的材料成形模擬大多仍停留在單機串行求解問題的水平上。隨著模擬工程的規模不斷增長，工程應用對材料成形模擬軟件的性能要求不斷提高。為了滿足實際需求，需要實現并行化，以提高性能和求解大型問題的能力[1]，因此如何提高算法的執行效率已經成為當前的一個重要的研究課題。在最近十年中，傳統的單CPU串行算法或者MPI并行算法被廣泛用來計算熱傳導中粒子的平均熱流和平均溫度，當粒子的數量太大時，這2種算法的執行時間太長。雖然圖形處理單元(Graphic Processing Unit, GPU)原本只用在圖形處理方面，但是現在它們已經越來越多地用于通用的科學和工程計算[2]，盡管如此，關于熱傳導的GPU并行處理的研究仍然不多。文獻[3]利用GPU的強大浮點運算能力來求解有關線形熱傳導和各向異性擴散的有限元方程。文獻[1]研究了二維熱傳導GPU并行化處理方案，總結出一個最優策略：將整個粒子計算域劃分成恰當數量的網格塊，每個網格塊粒子的計算由GPU并行執行，粒子處理效率將會提高。

本文通過參考已有的研究，即利用GPU強大的并行處理和浮點運算能力，以有效地提高熱傳導算法的執行效率，并基于CUDA編程環境提出一種熱傳導GPU并行算法。

2 熱傳導CPU串行計算研究

2.1 熱傳導算法的理論背景

熱傳導指的是，由大量物質的分子熱運動互相撞擊，溫度較高的分子(粒子)與相鄰的溫度較低的分子碰撞，而使能量從物體的高溫部分傳至低溫部分，或由高溫物體傳給低溫物體，直到溫度(熱流)平衡的過程[4]。

由熱的基礎理論知道，熱和物體的內能總是有密切的關系，在熱力學中，物體的內能是屬于物質內分子和原子的取向以及與運動有關的能量。因此，物體的熱傳導本質上和組成物體的微觀粒子的作用力、運動位移、運動速度等密切相關。本文熱傳導系統中的研究對象是分子結構(粒子)。

2.2 C PU串行算法的實現

CPU串行算法是過去解決熱傳導問題常用的算法，根據熱傳導的理論的計算特點，本文用一維分子結構系統來模擬一維熱傳導的過程，算法的流程如圖1所示。

圖1 C PU串行算法流程

CPU串行算法主要步驟如下：

(1)先初始化速度、作用力、位移、熱流、溫度等數組。

(2)讀入所用到的初始變量值，時間間隔dt=0.01，質量m=1.0，系數k=1.0，kb=1.0，第一個粒子的初始溫度T1=2，最后一個粒子的初始溫度T2=1。

(3)分步計算每個粒子的位移、速度、相互作用力，以及溫度和熱流；由于龍格-庫塔(Runge-Kutta)方法可以提高算法的精確度[5]，因此采用龍格-庫塔法求粒子的位移和速度。

(4)計算每個粒子的平均熱流和平均溫度。

3 GPU并行計算

3.1 GP U與CPU架構

隨著計算機技術的不斷發展，人們嘗試在每個計算節點裝備多路CPU，而在每個CPU芯片中也開始采用多核技術來提升CPU芯片的整體計算效率。CPU和GPU都是具有運算能力的芯片，但兩者功能的側重點不同，CPU是一個計算機的核心部件，它既執行指令運算又可以進行數值運算。GPU是目前使用很廣泛的顯卡的圖形處理器，已經成為強大的通用并行處理單元，它是顯卡的核心，專注于浮點運算、并行數值計算等方面。圖2說明了CPU和GPU在架構上的區別[6]。

圖2 C PU和GPU的架構區別

從圖2中看出，CPU進行串行運算，負責處理一條一條的數據，控制單元、邏輯單元和存儲單元三大部分相互協調，進行分析、判斷、運算并控制計算機各部分協調工作。GPU使用更多的單元進行并行數值計算，從收到數據到完成所有處理的全過程都可以完全獨立的并行處理。因此，對于大規模運算平行處理，GPU相對CPU具有優勢[7]。

CUDA是NVIDIA推出的通用并行計算架構，廣泛應用在Windows、Linux等不同的應用環境。CUDA采用C語言作為編程語言，且提供大量高性能計算指令的開發能力，可在GPU強大計算能力的基礎上建立一種效率更高的密集數據計算解決方案[8]。CUDA使得GPU對科學應用來說完全可編程并對高級語言(如C、C++、Fortran等)加以支持。CUDA將GPU作為數據并行計算設備，操作系統的多任務機制負責管理多個并發運行的CUDA和圖像應用程序對GPU的訪問，應用程序中精深復雜的計算和圖形部分運行在GPU上，而簡單的、連續的部分運行在CPU上。

3.2 GP U-CPU并行計算

GPU-CPU協作計算模式已成為普通PC機所擁有的強大、高效的計算模式。根據NVIDIA公司最近的測試結果顯示：利用GPU-CPU模式進行傅里葉變換、排序等科學計算比單獨用CPU計算的速度提高了19倍[9]。但是GPU-CPU協同模式也存在很多未解決的問題。雖然GPU在通用計算中有了很大的發展，但是和真正通用的CPU在結構、指令流處理等方面差別還較大。其次，用GPU-CPU模式進行并行計算缺少統一的編程模型。

4 熱傳導并行算法

4.1 M PI——一種典型的熱傳導并行算法

在GPU通用計算未普及之前，MPI并行算法是最常用的熱傳導算法。MPI并行算法是基于消息傳遞模式的并行計算方法。在MPI這種模式下，使用一個或多個進程間的通信調用庫函數發送、接收消息從而完成并行計算。消息傳遞模式使用通信協同的方式，源進程調用庫函數發送數據，目的進程調用相應的庫函數接收數據[10]。

MPI的優點在于可在任何并行機上運行，并且用戶可以顯式地操作并行程序的存儲；但是熱傳導MPI并行算法也很明顯：(1)消息傳遞接口的不同有時對算法產生很大的影響；(2)實際消息傳遞中的帶寬、延遲以及計算于通信的重疊會極大地影響并行算法性能；(3)MPI算法對硬件的要求比較高，而一般熱傳導算法在應用上可能用不到那么高要求的硬件。

根據以上分析知道，熱傳導MPI并行算法的不足是由MPI并行計算的特點決定，因此在MPI機制下無法改進。為了提高算法運行的可靠性及降低對硬件的要求，本文提出一個新的并行計算的方法——GPU并行算法。

4.2 熱傳導GPU并行算法

在本文提出的GPU并行算法[11]中，GPU并行是指GPU硬件上的并行，CPU負責控制程序的主流程，由GPU負責并行計算。該算法能彌補MPI并行算法存在的缺陷，解決以往熱傳導算法存在的問題，同時在較小的硬件開銷下大大提高算法的運行效率。

要進行GPU并行計算，必須具備以下的3個條件：

(1)一臺具備GPU并行計算平臺的計算機，該計算機具有支持CUDA編程的顯卡。

(2)并行的算法必須具備可并行性，可以找出可并行執行的子任務。

(3)在上述2個條件滿足的前提下，在GPU并行計算機上編寫GPU并行程序，并且運行該程序。

4.2.1 算法原理

本文針對GPU的計算特點及熱傳導粒子間相互作用關系的特點，設計一種GPU并行程序采用元胞單元法來組織數據。元胞單元法是將模擬計算區域劃分為一系列線程網格Grid，每個線程網格中包含多個粒子；由于粒子的狀態只和相鄰的2個粒子狀態有關，它具有近程性，因此與一個粒子有關的粒子只需要在同一個線程網格或是相鄰的線程網格中尋找即可，而不需要在整個計算區域里尋找[12]。

考慮到如果粒子在不同的線程網格中，它的狀態計算會使對應的內存開銷以及粒子狀態傳遞的操作數增多，因此，本文算法將所有的粒子劃分在一個線程網格Grid中，一個粒子的計算對應一個線程。多個線程組成一個Block(塊)，Block(塊)數目和每個Block中的Thread數目可根據粒子的具體數目進行調整，以達到最好的算法效率。算法設計在GPU中將粒子按照數目規模劃分為若干個網格，每執行一步共享內存中數據同步一次。200 000步后計算結果數據傳回CPU，計算出每個粒子的平均熱流量。

考慮到本文所設計的GPU并行算法要求在同一個線程網格內進行操作，可將對粒子信息的存儲用預先固定大小的數組存儲線程網格內粒子的狀態信息，這種方式可以很好地避免上一種方案的缺陷，還可以連續、并行地訪問內存。但是這些數組的大小必須夠大，可以存儲所有粒子的所有狀態信息。這必然導致這些數組占用的內存信息較大。

上述設計得到GPU并行算法的內核函數執行時的線程塊、線程于計算粒子間的對應關系如圖3所示。

圖3 內核函數執行時的線程塊、線程與粒子間的對應關系

基于上述理論，可以把整個熱傳導GPU并行算法的執行過程描述如下：(1)初始化及粒子狀態信息讀入；(2)建立粒子到線程網格的映射；(3)在具體的每個步中計算粒子的狀態。設計的熱傳導GPU并行算法的算法流程如圖4所示。

圖4 熱傳導GPU并行算法流程

在此算法的3個主要部分中，計算粒子的狀態是整個算法中耗時最長的部分，約占據整個算法的80%，這部分由GPU完成，而其他部分由CPU完成，這正好體現了GPU程序設計的初衷：將并行部分放在GPU上執行，串行主線部分放在CPU上執行。由于GPU不能直接訪問主儲存器，而CPU不能直接訪問顯存，因此計算完粒子狀態后，需要將顯存中的數據信息傳到CPU中，以便CPU完成接下來的操作。為了避免每個時間步都進行CPU到GPU及GPU 到CPU的2次數據交換操作，CPU在開始計算前讀入所有參與計算的粒子狀態信息，再將這些信息映射到相應的線程塊，傳遞給GPU，由GPU處理完算法中的所有模塊。

4.2.2 線程分配的設計

在CUDA中所用到的線程只是輕量級的線程，因此可以有上千個線程同時在執行。在GPU并行算法中，對于相同的粒子數目，每個線程網格內的粒子數越多，相應的線程塊數越少，則計算的時間越少。在GPU中，每個線程網格Grid中的線程塊Block數最多有65 535×65 535個，每個線程塊中的線程個數最大限制為512個。根據上述GPU特性，本文算法設計線程網格數為1，線程塊數也盡量少。每個粒子各占一個線程，這樣使得GPU中的資源可以實現最大化的并行計算。

5 實驗結果與分析

5.1 實驗結果

實驗環境為聯想S20的圖形工作站；CPU：Intel(R) Xeon(R) CPU W3550@3.07 GHz 3.06 GHz；GPU：NVIDIA Quadro FX 1800。在同一個實驗平臺下分別對本文所涉及的2種算法——熱傳導CPU串行程序及熱傳導GPU并行程序進行測試及對比。

對于CPU串行算法，主要測試它在不同粒子數時的時間性能。粒子數取值分別是16、160、1 600、16 000、160 000，測試結果如表1所示。

表1 C PU串行算法的時間性能測試結果

對于GPU并行算法，主要測試它在不同粒子數時的總時間性能、通信協同時間及GPU上的并行計算時間。粒子數取值同CPU串行算法一致，所得測試結果如表2所示。

表2 G PU并行算法的時間性能測試結果

5.2 串行和并行算法在時間效率上的比較

根據5.1節的實驗數據比較2種算法的時間效率，如圖5所示，粒子的個數分別為16×10、16×100、16×1 000、 16×10 000、16×100 000。從圖中容易看出，當系統中粒子數目不大時，并行算法沒有較明顯的優勢，但隨著粒子數目成指數倍地增長時，GPU并行算法處理粒子的效率遠高于CPU串行算法。

圖5 2種算法的時間效率比較

5.3 GP U并行算法的性能分析

GPU并行算法的性能優劣主要由加速比及并行算法的伸縮性決定，因此，本文實驗重點測試并比較這2個因素。

(1)加速比

加速比的公式為：C(n)=S(n)/P(n)，其中，S(n)是CPU串行算法的運行時間；P(n)是GPU并行算法的運行時間。

根據5.2節的實驗數據，得到加速比如表3所示。

表3 G PU并行算法在不同粒子數下的加速比

由表3可知，隨著粒子數的增加，GPU并行算法的加速比大幅增加，因此可以得出，在有較大規模粒子數時，GPU并行算法能更好地發揮其優勢。

(2)并行算法的伸縮性

當處理器數目一定時，算法規模越大，加速比越高，則稱此并行算法具有良好的伸縮性。從實驗結果上很容易看出，本文的熱傳導GPU并行算法具有伸縮性。

6 結束語

本文提出一個GPU并行算法，并基于CUDA的編程環境使用C語言實現了一維體系熱傳導算法。通過對一維分子結構體系模型進行研究，引入龍格庫塔法，運用它求出每一個粒子的速度、力、位移，最后計算出每個粒子的平均熱流量。通過實驗對比分析，結果表明，在大批量處理粒子數時，相對于傳統CPU串行算法，GPU并行算法在加速比方面可以提高近900倍，并且加速比隨著粒子數的增加而加速提高。這充分證明了GPU的并行處理能力的優勢。

GPU的發展勢頭異常迅猛，隨著時間的推移，以及人們對計算機并行處理能力要求的進一步提高，通過基于CUDA的GPU并行處理熱問題會越來越多，因此，如何提高GPU的并行計算能力將是未來學者們的一個重要研究課題。相信隨著硬件的不斷發展，GPU的處理能力也將進一步提高。雖然本文提出的一維熱傳導GPU并行算法較之前的CPU串行算法在時間效率上提高了幾百倍，但這還遠沒有發揮出GPU并行計算的優勢。若將此算法擴展為GPU集群上的并行算法，那么計算能力會比單GPU的算法能力強得多，下一步將研究工作環境更高的加速比。

[1] 王 梁. 二維穩態熱傳導問題CPU/GPU并行求解[EB/OL]. [2013-05-13]. http://tech.it168.com/a2010/0722/1081/0000010 81200.shtml.

[2] Frezzotti A, Ghiroldi G P. Solving the Boltzmann Equation on GPUs[EB/OL]. (2010-05-28). http://arxiv.org/abs/1005.5405.

[3] Rumpf M, Strzodka R. Usin g Graphics Cards for Quantized FEM Computations[C]//Proceedings of VIIP’01. Marbella, Spain: [s. n.], 2001: 98-107.

[4] 百度百科. 熱傳導[EB/OL]. [2013-05-13]. http://baike.baidu. com/view/348360.htm.

[5] 李夏云, 陳傳淼. 用龍格-庫塔法求解非線性方程組[J].數學理論與應用, 2008, 28(2): 62-65.

[6] 林 茂, 董玉敏, 鄒 杰. GPGPU編程技術初探[J]. 軟件開發與設計, 2010, (2): 15-17, 23.

[7] 孫敏杰. C PU架構和技術的演變看GPU未來發展[EB/OL]. [2013-05-13]. http://www.pcpop.com/doc/0/521/521832_all.sh tml.

[8] 孟小華, 劉堅強，區業祥. 基于CUDA的拉普拉斯邊緣檢測算法[J]. 計算機工程, 2012, 38(18): 190-193.

[9] 李 超. 論GPU-CPU協作計算模式的應用研究[J]. 電子商務, 2010, (11): 54.

[10] 郁志輝. 高性能計算并行編程技術-MPI并行程序設計[M].北京: 清華大學出版社, 2001.

[11] 朱麗莎. 基于GPU的一維體系熱傳導算法研究[D]. 廣州:暨南大學, 2011.

[12] 朱宗柏. 多重網格方法的并行化及其在傳熱數值分析中的應用[J]. 武漢交通科技大學學報, 2000, 24(4): 351-354.

編輯 任吉慧

Research on GPU Parallel Algorithm of Heat Conduction Based on CUDA

MENG Xiao-huaa,b, HUANG Cong-shana, ZHU Li-shaa,b

(a. Department of Computer Science; b. Sino-France Joint Laboratory for Astrometry, Dynamics and Space Science, Jinan University, Guangzhou 510632, China)

For real applications processing lar ge volume of particles in one-dimensional heat conduction problem, the response time of CPU serial algorithm and MPI parallel algorithm is too long. Considering Grap hic Processing Unit(GPU) of fers powerful paral lel processing capabilities, it implements a GPU parallel heat conduction algorithm on Compute Unifie d Device Architecture(CUDA) parallel programming environment using CPU and GPU collaborative mode. The algorithm sets the block and grid size based on GPU hardware configuration. Particles are divided into a plurality of blocks, the particle is into the GPU graphics for parallel computing, and one thread performs a calculation of a particle. It retrieves the processed data to CPU main me mory and calculates the average heat flow o f each particle. Experimental results sho w that, compared with CPU serial algorithm, GPU parallel algorithm has a great advantage in t ime efficiency, the speedup is close to 900, and speedup can improve as the particle number size increases.

heat con duction algorithm; Graph ic Processing Unit(GPU); C ompute U nified D evice A rchitecture(CUDA); parallel computing; time efficiency; speedup ratio

10.3969/j.issn.1000-3428.2014.05.009

國家自然科學基金資助項目(61073064)。

孟小華(1965－)，男，副教授、碩士，主研方向：并行分布式系統；黃叢珊，碩士研究生；朱麗莎，碩士。

2013-02-26

2013-05-21E-mail：xhmeng@163.com

1000-3428(2014)05-0041-04

A

TP399