孟小華，黃叢珊，朱麗莎

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

基于CUDA的熱傳導(dǎo)GPU并行算法研究

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

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

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

熱傳導(dǎo)算法；圖形處理單元；統(tǒng)一計算設(shè)備架構(gòu)；并行計算；時間效率；加速比

1 概述

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

本文通過參考已有的研究，即利用GPU強(qiáng)大的并行處理和浮點(diǎn)運(yùn)算能力，以有效地提高熱傳導(dǎo)算法的執(zhí)行效率，并基于CUDA編程環(huán)境提出一種熱傳導(dǎo)GPU并行算法。

2 熱傳導(dǎo)CPU串行計算研究

2.1 熱傳導(dǎo)算法的理論背景

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

由熱的基礎(chǔ)理論知道，熱和物體的內(nèi)能總是有密切的關(guān)系，在熱力學(xué)中，物體的內(nèi)能是屬于物質(zhì)內(nèi)分子和原子的取向以及與運(yùn)動有關(guān)的能量。因此，物體的熱傳導(dǎo)本質(zhì)上和組成物體的微觀粒子的作用力、運(yùn)動位移、運(yùn)動速度等密切相關(guān)。本文熱傳導(dǎo)系統(tǒng)中的研究對象是分子結(jié)構(gòu)(粒子)。

2.2 C PU串行算法的實現(xiàn)

CPU串行算法是過去解決熱傳導(dǎo)問題常用的算法，根據(jù)熱傳導(dǎo)的理論的計算特點(diǎn)，本文用一維分子結(jié)構(gòu)系統(tǒng)來模擬一維熱傳導(dǎo)的過程，算法的流程如圖1所示。

圖1 C PU串行算法流程

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

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

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

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

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

3 GPU并行計算

3.1 GP U與CPU架構(gòu)

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

圖2 C PU和GPU的架構(gòu)區(qū)別

從圖2中看出，CPU進(jìn)行串行運(yùn)算，負(fù)責(zé)處理一條一條的數(shù)據(jù)，控制單元、邏輯單元和存儲單元三大部分相互協(xié)調(diào)，進(jìn)行分析、判斷、運(yùn)算并控制計算機(jī)各部分協(xié)調(diào)工作。GPU使用更多的單元進(jìn)行并行數(shù)值計算，從收到數(shù)據(jù)到完成所有處理的全過程都可以完全獨(dú)立的并行處理。因此，對于大規(guī)模運(yùn)算平行處理，GPU相對CPU具有優(yōu)勢[7]。

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

3.2 GP U-CPU并行計算

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

4 熱傳導(dǎo)并行算法

4.1 M PI——一種典型的熱傳導(dǎo)并行算法

在GPU通用計算未普及之前，MPI并行算法是最常用的熱傳導(dǎo)算法。MPI并行算法是基于消息傳遞模式的并行計算方法。在MPI這種模式下，使用一個或多個進(jìn)程間的通信調(diào)用庫函數(shù)發(fā)送、接收消息從而完成并行計算。消息傳遞模式使用通信協(xié)同的方式，源進(jìn)程調(diào)用庫函數(shù)發(fā)送數(shù)據(jù)，目的進(jìn)程調(diào)用相應(yīng)的庫函數(shù)接收數(shù)據(jù)[10]。

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

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

4.2 熱傳導(dǎo)GPU并行算法

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

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

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

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

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

4.2.1 算法原理

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

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

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

上述設(shè)計得到GPU并行算法的內(nèi)核函數(shù)執(zhí)行時的線程塊、線程于計算粒子間的對應(yīng)關(guān)系如圖3所示。

圖3 內(nèi)核函數(shù)執(zhí)行時的線程塊、線程與粒子間的對應(yīng)關(guān)系

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

圖4 熱傳導(dǎo)GPU并行算法流程

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

4.2.2 線程分配的設(shè)計

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

5 實驗結(jié)果與分析

5.1 實驗結(jié)果

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

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

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

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

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

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

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

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

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

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

(1)加速比

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

根據(jù)5.2節(jié)的實驗數(shù)據(jù)，得到加速比如表3所示。

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

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

(2)并行算法的伸縮性

當(dāng)處理器數(shù)目一定時，算法規(guī)模越大，加速比越高，則稱此并行算法具有良好的伸縮性。從實驗結(jié)果上很容易看出，本文的熱傳導(dǎo)GPU并行算法具有伸縮性。

6 結(jié)束語

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

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

[1] 王 梁. 二維穩(wěn)態(tài)熱傳導(dǎo)問題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] 百度百科. 熱傳導(dǎo)[EB/OL]. [2013-05-13]. http://baike.baidu. com/view/348360.htm.

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

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

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

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

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

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

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

[12] 朱宗柏. 多重網(wǎng)格方法的并行化及其在傳熱數(shù)值分析中的應(yīng)用[J]. 武漢交通科技大學(xué)學(xué)報, 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

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

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

2013-02-26

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

1000-3428(2014)05-0041-04

A

TP399