基于OpenCL的Lammps短程力算法優(yōu)化研究*

趙成龍 ，施慧彬 ，俞忻峰

(1.南京航空航天大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,江蘇 南京 210016;2.南京理工大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,江蘇 南京 210094)

基于OpenCL的Lammps短程力算法優(yōu)化研究*

趙成龍1，施慧彬1，俞忻峰2

(1.南京航空航天大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,江蘇 南京 210016;2.南京理工大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,江蘇 南京 210094)

Lammps是用于分子動(dòng)力學(xué)模擬及其相關(guān)問題的一款開源軟件，可利用其了解固體、液體性質(zhì)，應(yīng)用廣泛。支持使用CUDA及OpenCL進(jìn)行GPU加速。因OpenCL具有跨平臺(tái)特性，將其作為研究重點(diǎn)?？偨Y(jié)了OpenCL內(nèi)核編程中需要注意的設(shè)計(jì)原則并闡述了一種改進(jìn)的阿姆達(dá)爾定律用于衡量異構(gòu)平臺(tái)理論加速性能。測(cè)試了Lammps短程力計(jì)算在Y485P平臺(tái)下的性能參數(shù)。通過對(duì)短程力計(jì)算中的關(guān)鍵部分如鄰接表的建立及短程力計(jì)算部分的內(nèi)核代碼進(jìn)行優(yōu)化，使其取得了更好的加速效果。

GPU;Lammps;OpenCL;分子動(dòng)力學(xué);短程力;內(nèi)核優(yōu)化

1 引言

分子動(dòng)力學(xué)模擬的基本原理為建立一個(gè)合適的粒子系統(tǒng)，利用牛頓運(yùn)動(dòng)方程確定系統(tǒng)中各個(gè)粒子的運(yùn)動(dòng)情況，通過求取粒子動(dòng)力學(xué)方程組的數(shù)值解，決定系統(tǒng)中的各個(gè)粒子在相空間中的運(yùn)動(dòng)規(guī)律，最終通過統(tǒng)計(jì)物理及熱力學(xué)原理得出系統(tǒng)相應(yīng)的宏觀物理特性。其初始研究可以追溯到1957年Alder B J和Wainwright T E的硬球模型[1]。隨著其發(fā)展，應(yīng)用領(lǐng)域逐漸涵蓋化學(xué)、物理、能源、生物、材料、醫(yī)學(xué)等多個(gè)領(lǐng)域，因其運(yùn)算量龐大、耗時(shí)長的特點(diǎn)，一直制約了其發(fā)展。

大規(guī)模原子分子并行模擬器Lammps(Large-scale Atomic/Molecular Massively Parallel Simulator)主要用在分子動(dòng)力學(xué)相關(guān)的一些計(jì)算和模擬上，是美國Sandia國家實(shí)驗(yàn)室開發(fā)并以GPL licence發(fā)布的開放源代碼的分子動(dòng)力學(xué)軟件。Lammps可以支持氣態(tài)、液態(tài)或者固態(tài)的各種系統(tǒng)下，千萬級(jí)以及上億級(jí)規(guī)模的原子分子體系，并支持多種勢(shì)函數(shù)，具有良好的并行擴(kuò)展性[2]。

除Lammps之外,還有流行的分子動(dòng)力學(xué)模擬軟件如NAMD[3]、Amber[4]、Gromacs[5]也采用了GPU進(jìn)行加速。OpenCL因其可跨平臺(tái)使用的特點(diǎn)越來越多地被分子動(dòng)力學(xué)軟件采用。 Lammps中基于OpenCL的GPU加速包由Brown W M等在文獻(xiàn)[6～8]中提出。同時(shí)，基于OpenCL的分子動(dòng)力學(xué)模擬加速還有文獻(xiàn)[9]提出的一種性能自適應(yīng)的分子動(dòng)力學(xué)模擬內(nèi)核。對(duì)于OpenCL內(nèi)核優(yōu)化工作，目前大多數(shù)的研究主要針對(duì)內(nèi)存的讀取、資源利用率及算法本身等的優(yōu)化，而對(duì)內(nèi)核代碼的優(yōu)化，目前沒有出現(xiàn)過詳細(xì)闡述，如文獻(xiàn)[10]闡述了訪存效率、資源利用率、硬件資源限制對(duì)性能的影響及優(yōu)化方式。文獻(xiàn)[11,12]僅對(duì)內(nèi)核編程中需要注意的極少部分進(jìn)行了闡述。

本文詳細(xì)總結(jié)了OpenCL內(nèi)核編程過程中的原則。闡述了一種改進(jìn)的阿姆達(dá)爾定律，可用于異構(gòu)平臺(tái)下的性能評(píng)估。通過研究Lammps中的短程力LJ(Lennard-Jhones)力實(shí)現(xiàn)部分，對(duì)鄰接表建立、力的計(jì)算部分分別進(jìn)行優(yōu)化，對(duì)比優(yōu)化前后的性能,驗(yàn)證了內(nèi)核程序的優(yōu)化有助于算法性能的提升。

2 OpenCL內(nèi)核優(yōu)化編程

2.1 OpenCL簡(jiǎn)介及平臺(tái)信息

開放計(jì)算語言O(shè)penCL[13](Open Computing Language)是第一個(gè)面向異構(gòu)系統(tǒng)的通用并行編程標(biāo)準(zhǔn)，其特點(diǎn)是免費(fèi)、跨平臺(tái)以及良好的互操作特性。它為軟件開發(fā)人員提供了統(tǒng)一的編程環(huán)境，廣泛適用于多核處理器(CPU)、圖形處理器(GPU)、Cell類型架構(gòu)等并行處理器，擁有廣闊的發(fā)展前景。其標(biāo)準(zhǔn)可利用平臺(tái)模型、執(zhí)行模型、內(nèi)存模型、編程模型進(jìn)行描述。

本文所使用的測(cè)試平臺(tái)為聯(lián)想Y485P筆記本，CPU型號(hào)為A10-5750M，主頻2.5 GHz。平臺(tái)顯卡信息如表1所示，其中HD8650G為集顯。測(cè)試系統(tǒng)為Ubuntu 13.10。

2.2 并行加速理論模型

在并行編程中，衡量多核環(huán)境下并行應(yīng)用在固定負(fù)載下并行處理效果的加速比S，通常使用阿姆達(dá)爾定律(Amdahl’s Law)。該定律描述了系統(tǒng)中

Table1 Test platform information表1 測(cè)試平臺(tái)信息

并行以及串行所占系統(tǒng)計(jì)算比重對(duì)于加速比的影響,公式如下：

(1)

其中，s為系統(tǒng)中串行部分所占比例，p為整個(gè)應(yīng)用中并行所占比例(為1-s)，N為處理器個(gè)數(shù)。

阿姆達(dá)爾定律對(duì)多核系統(tǒng)中負(fù)載恒定的理想情況下是成立的，而當(dāng)N→∞時(shí)，取得最大值：

S=1/s

(2)

式(2)說明，理想情況下，即使串行部分的工作很小，系統(tǒng)的加速比也只能夠限制在1/s。

阿姆達(dá)爾定律的首要條件是系統(tǒng)中的信息傳輸開銷可以忽略。但是，在使用GPU作為加速器件的情況這種開銷將變得不可忽視。從平臺(tái)整體架構(gòu)來說，使用GPU進(jìn)行加速的應(yīng)用可以作為非對(duì)稱處理器的一種。文獻(xiàn)[14]就這種情況下對(duì)阿姆達(dá)爾定律進(jìn)行了改進(jìn)，引入了兩個(gè)新的參數(shù)：

(1)已并行加速比例p′:對(duì)于傳統(tǒng)對(duì)稱多核處理器，并行比例p最優(yōu)化后的比例為p/N。然而，非對(duì)稱多核處理器的加速內(nèi)核相較串行內(nèi)核在計(jì)算能力上有很大的不同，故使用p和N在公式中將變得沒有任何意義。而真正在GPU上所實(shí)現(xiàn)的效果也會(huì)因?yàn)樗惴捌溆成涞郊铀賰?nèi)核上的方式的不同而有很大的區(qū)別。所以，應(yīng)當(dāng)將p轉(zhuǎn)換成p′。

(2)并行損耗o:在傳統(tǒng)對(duì)稱多核處理器情況下，這種并行損耗通常是可以忽略的，因?yàn)橥ǔ?nèi)核間的通信代價(jià)很小。而對(duì)于加速器件平臺(tái)(如GPU)來說，這部分的損耗就變得不可忽視，這主要是因?yàn)镃PU內(nèi)存與GPU內(nèi)存之間的數(shù)據(jù)傳輸時(shí)間占程序執(zhí)行時(shí)間的比重較大。

最終的加速比S′可以表示為：

(3)

從式(3)中可看出，若希望系統(tǒng)的加速比盡可能高，則必須找出系統(tǒng)中足夠多的并行化部分代碼，盡量減小s。同時(shí)，需要注意到的是，由于GPU的數(shù)據(jù)傳輸造成o比較大，將弱化并行化p′所獲得的加速效果,這在后續(xù)的實(shí)驗(yàn)數(shù)據(jù)中也可以看出。這也為OpenCL的性能優(yōu)化提供了一定的理論基礎(chǔ)。

利用AMD公司提供的SDK中的BufferBandWidth例子，測(cè)得系統(tǒng)中兩個(gè)GPU與CPU在不同數(shù)據(jù)大小(從2 KB～64 MB)情況下的傳輸速率。最終結(jié)果如圖1所示。

Figure 1 Bandwidth test

圖1a描述了從主機(jī)(CPU)到設(shè)備(GPU)的帶寬?？梢钥吹奖M管HD8650G作為AMD推出的加速處理器APU(Accelerated Processing Unit)中的集顯，排除了PCIE訪問，但是仍然可以看出其數(shù)據(jù)傳輸相較HD8790M的獨(dú)顯來說，只有在超出一個(gè)閾值(Threshold)時(shí)這種優(yōu)勢(shì)才變得比較明顯。圖1b描述了相反方向的測(cè)試結(jié)果，即從設(shè)備(GPU)到主機(jī)(CPU)的測(cè)試結(jié)果。設(shè)備到主機(jī)的傳輸，該閾值比主機(jī)到設(shè)備要小。同時(shí)，設(shè)備到主機(jī)的最大帶寬也略高于主機(jī)到設(shè)備的帶寬。

這也從側(cè)面印證了對(duì)于APU設(shè)備來說，在小數(shù)據(jù)量的情況下，數(shù)據(jù)的傳輸并沒有克服PCIE帶寬所帶來的限制，但是從大數(shù)據(jù)量上來說，情況要好很多。這可能是由于數(shù)據(jù)在傳輸之前需要進(jìn)行DMA傳輸?shù)慕⒁约皟?nèi)存的固定(Pinning)操作有關(guān)。同時(shí)，對(duì)于本測(cè)試平臺(tái)，HD8650G可以支持AMD的零拷貝特性(Zero Copy)但是直接的帶寬量難以測(cè)得，使用Map(clEnqueueMapBuffer)內(nèi)存的方式對(duì)寫入和讀取操作進(jìn)行了進(jìn)一步的測(cè)試，得到使用memset寫零拷貝內(nèi)存的帶寬可達(dá)到4.179 GB/s，而讀取操作則只有1.214 GB/s。

為了能夠更好地發(fā)揮APU零拷貝特性，在使用APU進(jìn)行異構(gòu)加速運(yùn)算時(shí)，對(duì)于內(nèi)存的創(chuàng)建和讀寫建議通過以下方式進(jìn)行：

(1)在創(chuàng)建內(nèi)存對(duì)象時(shí)，除基本可訪問特性外，通過增加CL_MEM_ALLOC_HOST_PTR、CL_MEM_USE_HOST_PTR或CL_MEM_USE_PERSISTENT_MEM_AMD進(jìn)行內(nèi)存創(chuàng)建。CL_MEM_ALLOC_HOST_PTR和CL_MEM_USE_HOST_PTR為OpenCL 1.2標(biāo)準(zhǔn)中所定義的內(nèi)存分配方式，CL_MEM_USE_PERSISTENT_MEM_AMD為AMD平臺(tái)擴(kuò)展實(shí)現(xiàn)的內(nèi)存分配方式。通過這些方式創(chuàng)建的內(nèi)存具備零拷貝特性。

(2)在讀取和寫入時(shí)，不采用clEnqueueReadBuffer以及clEnqueueWriteBuffer的方式進(jìn)行內(nèi)存訪問，使用clEnqueueMapBuffer以及clEnqueueUnmapMemObject的方式進(jìn)行數(shù)據(jù)讀寫操作。注意，對(duì)于圖像內(nèi)存對(duì)象來說，只有使用CL_MEM_USE_PERSISTENT_MEM_AMD特性創(chuàng)建的內(nèi)存對(duì)象才具備此特性。

通過上述方式，可以直接減少PCIE傳輸帶來的內(nèi)存?zhèn)鬏旈_銷，從而帶來傳輸性能的提升。同時(shí)，查閱OpenCL 2.0標(biāo)準(zhǔn)還可發(fā)現(xiàn)，新的內(nèi)存模型中集成的GPU和CPU將采用相同的內(nèi)存地址范圍，也就意味著真正零拷貝的實(shí)現(xiàn)。

2.3 內(nèi)核優(yōu)化原則

通過整理與測(cè)試，發(fā)現(xiàn)在OpenCL內(nèi)核編程中，也有一些需要在編碼過程中注意的地方，將其整理為如下內(nèi)核編程原則：

原則1展開小規(guī)模循環(huán)：如果已知一個(gè)循環(huán)的邊界值，并且該循環(huán)的規(guī)模是比較小的，比如指令為16～32個(gè)，可以展開該循環(huán)，這樣避免了循環(huán)中的比較操作。展開的操作可以使用手工的方式或者使用編譯命令#pragma unrollN展開循環(huán)。使用編譯器命令時(shí)必須確定在編譯時(shí)編譯器可以獲得循環(huán)的次數(shù)，否則不起作用。而使用手工方式則可以進(jìn)行更加精細(xì)的控制，比如加入向量化運(yùn)算特性、自主決定循環(huán)中同步操作的次數(shù)及位置。

原則2使用編譯器命令傳遞常量：如果在內(nèi)核中需要使用到一個(gè)常量，應(yīng)當(dāng)使用編譯條件進(jìn)行定義，而不是浪費(fèi)一個(gè)內(nèi)核參數(shù)進(jìn)行傳遞。這樣可以減少對(duì)本地內(nèi)存或者私有內(nèi)存的讀取操作。例如，你想傳入一個(gè)LEN的宏，可以通過如下方式：

clBuildProgram(program,0,NULL,”-DLEN=256”,NULL,NULL);

在編譯這個(gè)內(nèi)核時(shí)，程序?qū)?huì)把LEN的宏自動(dòng)定義為256。

原則3在無需進(jìn)行共享數(shù)據(jù)的情況下，使用私有內(nèi)存(Private Memory)來代替本地內(nèi)存(Local Memory)。工作項(xiàng)訪問私有內(nèi)存的速率比訪問本地內(nèi)存的速率要快，但是如果這些數(shù)據(jù)需要在工作組之間共享，則只能使用本地內(nèi)存。需要注意的是，如果使用OpenCL操作的設(shè)備是CPU，則該原則不一定成立。例如，AMD中CPU設(shè)備的OpenCL實(shí)現(xiàn)，使用Cache進(jìn)行模擬，則不一定會(huì)有優(yōu)化效果。

原則4盡量避免使用%操作符：求余操作符在顯卡上需要大量的處理時(shí)間，應(yīng)當(dāng)避免或者尋求其他的解決方式。

原則5內(nèi)聯(lián)化非內(nèi)核函數(shù)：使用inline聲明告訴編譯器，在調(diào)用函數(shù)處對(duì)函數(shù)進(jìn)行展開，可以免除上下文切換以及調(diào)用堆棧所造成的損失。

原則6使用內(nèi)核函數(shù)：這里舉一些例子進(jìn)行說明：

(1)if-then-else語句會(huì)造成分支，使用select函數(shù)避免分支的產(chǎn)生，如：

if(x== 1)r=0.5;

if(x== 2)r=1.0;

可以改成：

r=select(r,0.5,x==1);

r=select(r,1.0,x==2);

(2)浮點(diǎn)運(yùn)算的乘加運(yùn)算，如果對(duì)速度要求大于精度要求，可使用mad以及fma函數(shù)進(jìn)行乘加運(yùn)算。

(3)帶native_前綴的函數(shù)通常比不帶native_前綴或帶half_前綴的相應(yīng)函數(shù)要性能更優(yōu)，其精度由具體實(shí)現(xiàn)決定。

原則7使用分支預(yù)測(cè)而不是流程控制，如下面情況的代碼：

if(A>B){

C+=D;

}else{

C-=D;

}

可改成：

intfactor=(A>B)?:1:-1;C+=factor*D;

這是因?yàn)榍耙粋€(gè)代碼段會(huì)生成IF/ELSE/ENDIF的控制語句，每一個(gè)都將會(huì)花費(fèi)～8的周期，整個(gè)代碼可能花費(fèi)到～36個(gè)周期。而下面部分的代碼段只需要12個(gè)周期，可以提高至少1.3倍。

原則8通過順序訪問本地內(nèi)存避免訪問沖突：AMD的本地內(nèi)存一般由32個(gè)庫(Bank)組成，如果多個(gè)工作組訪問同一個(gè)庫，則訪問操作將會(huì)串行化。而當(dāng)同組中的相鄰工作項(xiàng)訪問相鄰的庫，可盡可能地減少串行化訪問帶來的損失。

原則9盡量采用float4、int4這樣的向量數(shù)據(jù)類型，在內(nèi)核參數(shù)中使用向量數(shù)據(jù)類型可以有效提高數(shù)據(jù)帶寬。而在內(nèi)核代碼中使用該數(shù)據(jù)類型也更符合GPU架構(gòu)并行化特點(diǎn)。同時(shí)，需要注意的是，對(duì)于AMD中的GCN(Graphics Core Next)架構(gòu)的GPU來說，由于設(shè)備架構(gòu)的改變，這一特性已經(jīng)弱化，但是在數(shù)據(jù)的傳輸時(shí)，采用向量類型仍然能夠獲取一定的性能提升。

原則10如一個(gè)表達(dá)式在內(nèi)核中重復(fù)出現(xiàn)，應(yīng)考慮用臨時(shí)變量存儲(chǔ)，以避免多次計(jì)算造成的性能損失。

原則11避免不必要的barrier操作，在內(nèi)核中barrier操作對(duì)性能影響很大，應(yīng)考慮調(diào)優(yōu)算法來減少barrier的使用。

上述是比較通用的11條編程原則。另外，還有一些需要注意的事項(xiàng)需在具體編碼中根據(jù)測(cè)試和并行化編程特點(diǎn)具體進(jìn)行分析才能決定。

同時(shí)，在編碼過程中，合理地選擇數(shù)據(jù)訪問和傳輸方式也能夠在一定程度上提升程序性能。不良的編碼方式可能導(dǎo)致數(shù)據(jù)的重復(fù)復(fù)制，浪費(fèi)帶寬。而從圖1中可以看出，數(shù)據(jù)的傳輸速率并不是很快。對(duì)于AMD平臺(tái)來說，AMD的OpenCL編程手冊(cè)[15]中根據(jù)不同的情況給出了建議，可以做為數(shù)據(jù)傳輸?shù)膮⒖肌?/p>

3 優(yōu)化及結(jié)果

本文采用Lammps自帶的LJ力計(jì)算測(cè)試腳本進(jìn)行測(cè)試，測(cè)試的粒子個(gè)數(shù)為262 144個(gè)，測(cè)試的步長分別為5 000和10 000步。

分別使用4CPU、4CPU加一個(gè)集成顯卡(HD8650G，4CPU+1GPU)以及加兩個(gè)顯卡(HD8650G+HD8790M，4CPU+2GPU)三種情況進(jìn)行測(cè)試，記錄其總計(jì)算時(shí)間Time、LJ力計(jì)算用時(shí)占的百分比PT(Pair Time)、建立鄰接表用時(shí)所占的百分比NT(Neigh Time)以及系統(tǒng)通信所占的百分比CT(Comm Time)。為了便于描述，將百分比部分的數(shù)據(jù)精度設(shè)定為一位，對(duì)于較小的數(shù)字采用直接設(shè)其為零。測(cè)試的初始結(jié)果如表2所示。

Table 2 Initial test results of LJ force表2 LJ力初始測(cè)試結(jié)果

可以看出，在只用CPU進(jìn)行運(yùn)算時(shí)，其主要的時(shí)間開銷為LJ力計(jì)算部分以及鄰接表的建立。而當(dāng)使用GPU進(jìn)行計(jì)算后，雖然在總時(shí)間上略少于只使用CPU進(jìn)行運(yùn)算，但是通信時(shí)間所占的比重增大了很多，這也導(dǎo)致鄰接表建立的時(shí)間占運(yùn)算總時(shí)間的比重微乎其微，如公式(3)所描述的那樣弱化了p′帶來的性能提升。同時(shí)，通過測(cè)試得到在LJ力計(jì)算部分加速比僅比使用4CPU的情況加速了1.2～1.24倍。本文后續(xù)將闡述對(duì)鄰接表建立的研究及優(yōu)化、LJ力計(jì)算研究及優(yōu)化并給出優(yōu)化后的結(jié)果進(jìn)行對(duì)比。因原Lammps代碼中，OpenCL加速部分是作為靜態(tài)庫進(jìn)行調(diào)用的，這樣每次更新測(cè)試時(shí)均需要對(duì)整個(gè)代碼進(jìn)行編譯。為了方便，我們首先改動(dòng)了OpenCL庫的編譯文件Makefile，使其生成動(dòng)態(tài)庫供Lammps主程序調(diào)用。

3.1 鄰接表建立研究及優(yōu)化

對(duì)于分子動(dòng)力學(xué)模擬來說，仿真的粒子數(shù)目通?？梢赃_(dá)到萬、百萬甚至更高。對(duì)每一個(gè)計(jì)算步來說，直接計(jì)算其他粒子對(duì)某一粒子的作用力，計(jì)算復(fù)雜度將達(dá)到O(n2)。從短程力本身具有截?cái)嗵攸c(diǎn)出發(fā)，人們研究出了兩種降低計(jì)算復(fù)雜度的方式，即元胞列表法和網(wǎng)格搜索法。

(1)元胞列表法:將系統(tǒng)中的所有粒子，劃分成N*N*N的網(wǎng)格，這樣只需計(jì)算每個(gè)網(wǎng)格中粒子與相鄰網(wǎng)格的粒子之間的作用力便可，比如二維情況下只需進(jìn)行8個(gè)鄰近元胞的搜索，而三維情況下為27個(gè)。網(wǎng)格大小通常選擇為截?cái)喟霃絩cut。

(2)網(wǎng)格搜索法：與元胞列表法類似，不同的是，不把鄰近網(wǎng)格中所有的粒子都加入當(dāng)前粒子的鄰居列表，而是將相同或者相鄰網(wǎng)格里的小于截?cái)喟霃降牧Ｗ蛹尤氲洁従恿斜碇?，進(jìn)一步減小鄰居列表的長度。

Lammps使用網(wǎng)格搜索法進(jìn)行粒子鄰接表的創(chuàng)建，分兩步進(jìn)行，首先建立元胞列表，然后使用網(wǎng)格搜索法建立鄰居列表?？梢杂枚S的模型進(jìn)行描述，如圖2所示。

Figure 2 Neighbor-list building

將所有的粒子進(jìn)行編號(hào)，并劃分元胞區(qū)域。通過粒子當(dāng)前坐標(biāo)，更新粒子所在的元胞位置。在建立鄰居列表時(shí)，首先更新元胞情況，以圖2為例，搜索圓圈區(qū)域內(nèi)部的粒子并加入到圓心粒子的鄰居列表中。為了建立元胞列表，Lammps使用了基數(shù)排序算法。具體過程為將粒子所在的元胞編號(hào)作為關(guān)鍵字，粒子編號(hào)作為記錄，按照關(guān)鍵字排序后，便可得出每個(gè)元胞中所具有粒子的有序表。我們首先使用OpenCL實(shí)現(xiàn)了一種并行基數(shù)排序算法。算法的核心是求前綴和(Prefix Sum)。需要說明的是，因鄰接表的建立并不是程序主要耗時(shí)部分，對(duì)減少計(jì)算總時(shí)間并沒有顯著優(yōu)勢(shì)。

Lammps鄰居列表建立部分OpenCL內(nèi)核代碼存儲(chǔ)在neighbor_gpu_cl.h頭文件中。優(yōu)化內(nèi)核實(shí)現(xiàn)代碼工作如下：

(1)利用中間變量去除代碼中求取元胞坐標(biāo)的取余(%)操作。源代碼中的計(jì)算過程為：

intix=BLOCK_ID_X+cells_in_cutoff;

intiy=BLOCK_ID_Y% (ncelly-cells_in_cutoff*2)+cells_in_cutoff;

intiz=BLOCK_ID_Y/ (ncelly-cells_in_cutoff*2)+cells_in_cutoff;

改進(jìn)后過程為：

intix=BLOCK_ID_X+cells_in_cutoff;

inttmp=ncelly-cells_in_cutoff*2;

intdevide=BLOCK_ID_Y/tmp;

intiz=devide+cells_in_cutoff;

intiy=BLOCK_ID_Y-devide*tmp+cells_in_cutoff;

這里可以看出違反了原則10和原則4，通過定義一個(gè)臨時(shí)變量存儲(chǔ)中間結(jié)果以及調(diào)整計(jì)算順序避免了取余操作。

(2)在查詢粒子i對(duì)應(yīng)的鄰接元胞jcell的序號(hào)時(shí)，需要對(duì)三個(gè)坐標(biāo)軸方向進(jìn)行掃描，在原代碼中采用了三重for循環(huán)進(jìn)行計(jì)算，并且在最內(nèi)層循環(huán)計(jì)算jcell，因計(jì)算jcell操作會(huì)帶來很多連乘及累加操作，有一定的性能損失。元胞序號(hào)有一定的規(guī)律，第一個(gè)優(yōu)化方式為使用序列值的規(guī)律進(jìn)行簡(jiǎn)化，可以將三重for循環(huán)簡(jiǎn)化成一重for循環(huán),jcell的計(jì)算也可以簡(jiǎn)化到加法運(yùn)算。這在串行程序中，可以獲得比較好的加速效果。但是，因循環(huán)過程有一些判斷語句，在異構(gòu)并行環(huán)境下反而造成了性能的損失，我們的測(cè)試也證明了這一點(diǎn)。通過分析，使用初始jcell，將jcell輪詢中的變化分割到每層循環(huán)中，獲得了較好的效果。源代碼中計(jì)算過程為：

for (intnborz=nborz0;nborz≤nborz1;nborz++)

for (intnbory=nbory0;nbory≤nbory1;nbory++)

for (intnborx=nborx0;nborx≤nborx1;nborx++)

intjcell=nborx+nbory*ncellx+nborz*ncellx*ncelly;

end for

end for

end for

改進(jìn)后過程為：

intnborz=nborz0,nbory=nbory0,nborx=nborx0;

inttmp2=ncellx*ncelly;

intjcell0=nborx0 +nbory0*ncellx+nborz0*tmp2;

intjcelly0=jcell0,jcellz0=jcell0;

intjcell=jcell0;

while(nborz≤nborz1)

while(nbory≤nbory1)

while(nborx≤nborx1)

//do scan

jcell++;

nborx++;

end while

nbory++;jcelly0+=nx;jcell=jcelly0;

end while

nborz++;

jcellz0+=tmp2;jcelly0=jcellz0;jcell=jcellz0;

end while

(3)根據(jù)原則9將源代碼中部分可使用向量方式優(yōu)化的代碼進(jìn)行優(yōu)化，如坐標(biāo)的運(yùn)算，源代碼采用的是標(biāo)量方式進(jìn)行計(jì)算。需要三條語句完成，而通過向量化操作后則只需一條語句即可實(shí)現(xiàn)。

3.2 LJ力計(jì)算研究及優(yōu)化

OpenCL加速庫中LJ力計(jì)算部分在lj_cl.h文件中。主要過程為根據(jù)查詢需計(jì)算的i粒子鄰接粒子j的范圍，輪詢?cè)摲秶械牧Ｗ有畔ⅰＭㄟ^使用公式(4)中的Lennard-Jones勢(shì)函數(shù)計(jì)算公式計(jì)算粒子間的作用力。

V(r)=4ε[(σ/r)12-(σ/r)6]

(4)

其中，ε為勢(shì)能阱的深度，σ是相互作用的勢(shì)能剛好為0時(shí)的兩體距離，即截?cái)喟霃絩cut。這兩個(gè)參數(shù)通過查找表得到。r為兩粒子間的距離。對(duì)該部分內(nèi)核的優(yōu)化如下：

(1)根據(jù)原則1展開小規(guī)模循環(huán)操作。在內(nèi)核代碼中的store_answers、store_answers_q、k_lj、k_lj_fast部分可以發(fā)現(xiàn)多處已知邊界的小循環(huán)。通過測(cè)定，相較原代碼，僅展開操作所帶來的LJ力計(jì)算用時(shí)性能提升可達(dá)到1.25～1.303倍。

(2)根據(jù)原則6將代碼中部分可以使用內(nèi)核函數(shù)優(yōu)化的代碼進(jìn)行優(yōu)化，如在代碼中存在一些使用條件判斷進(jìn)行賦值的單行語句，可以使用原則6所舉例子中的select內(nèi)核函數(shù)對(duì)其進(jìn)行優(yōu)化。

(3)根據(jù)原則10將代碼中重復(fù)計(jì)算部分通過提取到循環(huán)外部、使用臨時(shí)變量等方式進(jìn)行優(yōu)化。

(4)根據(jù)原則9將代碼中部分標(biāo)量操作轉(zhuǎn)換成向量運(yùn)算形式。

3.3 優(yōu)化結(jié)果

通過對(duì)短程力實(shí)現(xiàn)部分優(yōu)化后再次測(cè)試得到的數(shù)據(jù)如表3所示。

Table 3 LJ test results after optimization表3 LJ力優(yōu)化后測(cè)試結(jié)果

對(duì)比發(fā)現(xiàn)，相較原OpenCL實(shí)現(xiàn)代碼，改進(jìn)后的總時(shí)間在使用單GPU測(cè)試時(shí)5 000步及10 000步分別達(dá)到了原來的1.3倍及1.72倍。而采用雙GPU同時(shí)測(cè)試時(shí)，計(jì)算5 000步及10 000步分別達(dá)到了原來的1.39倍及1.74倍。對(duì)于LJ力計(jì)算部分，單GPU 5 000步提升了1.39倍，10 000步提升了1.48倍。雙GPU 5 000步提升1.75倍，10 000步提升了1.79。鄰接表建立部分代碼單GPU最高提升1.35倍，雙GPU最高提升1.64倍。

從上述結(jié)果可以看出，雖然GCN顯卡架構(gòu)從結(jié)構(gòu)上使得GPU更加適用于異構(gòu)編程，但通過合理的優(yōu)化代碼仍然能夠在性能上有所提升。我們發(fā)現(xiàn)影響內(nèi)核程序性能的主要原因在于原則1、原則8、原則9、原則11;而當(dāng)特定語句如原則4、原則6、原則7、原則10中相關(guān)情況出現(xiàn)較多時(shí)，也會(huì)造成性能的損失；其他原則對(duì)性能的影響并不明顯。其主要思想可以歸結(jié)為減少?zèng)_突、增加訪問帶寬以及合理利用資源。另外，針對(duì)如鄰接表建立中(2)中所述情況的特殊情況應(yīng)合理調(diào)整計(jì)算方式。為最大化性能提升，可繼續(xù)研究算法中可并行化部分，如本文中增加了鄰接表部分基數(shù)排序算法的OpenCL實(shí)現(xiàn)。由于OpenCL內(nèi)存模型中全局內(nèi)存和本地內(nèi)存讀取速度差異較大，而本地內(nèi)存所能夠存儲(chǔ)的信息容量一般為32 KB～64 KB，這在大規(guī)模分子并行模擬時(shí)，使得模擬粒子信息存儲(chǔ)在全局內(nèi)存中，加重了傳輸損耗，造成性能瓶頸。比較好的方式是，結(jié)合消息傳遞接口MPI(Message Passing Interface)進(jìn)一步分冶，減小單節(jié)點(diǎn)負(fù)擔(dān)，獲得更好的性能提升，這在文獻(xiàn)[2]中也有所體現(xiàn)。

4 結(jié)束語

本文首先歸納總結(jié)了在OpenCL內(nèi)核編程中可以提升性能的一些原則，闡述了一種衡量異構(gòu)平臺(tái)情況下加速效果的理論公式。首先對(duì)Lammps短程力在Y485P平臺(tái)下的加速效果進(jìn)行了測(cè)試，然后通過研究Lammps分子動(dòng)力學(xué)中的OpenCL短程力實(shí)現(xiàn)代碼，根據(jù)總結(jié)出的內(nèi)核代碼編寫原則進(jìn)行優(yōu)化，使得最終的總時(shí)間較原代碼最高提高了1.74倍，LJ力計(jì)算部分較原代碼最高提高了1.79倍，鄰接表計(jì)算較原代碼最高提高了1.64倍。證明通過合理優(yōu)化內(nèi)核代碼能夠提高其性能。

[1] Alder B J,Wainwright T E.Phase transition for a hard sphere system[J]. The Journal of Chemical Physics,1957,27(5):1208-1209.

[2] Li Bo-yang,Nie Feng-guang,Li Xiao-xia,et al.Performance test of lammps molecular dynamics simulation on GPU parallel computing cluster[J].Computers and Applied Chemistry,2011,28(10):1229-1233.(in Chinese)

[3] http://www.ks.uiuc.edu/Research/namd/.2014.

[4] http://ambermd.org/.2014.

[5] http://www.gromacs.org/.2014.

[6] Brown W M,Wang P,Plimpton S J,et al. Implementing molecular dynamics on hybrid high performance computers-short range forces[J]. Computer Physics Communications,2011,182(4):898-911.

[7] Brown W M,Kohlmeyer A,Plimpton S J,et al. Implementing molecular dynamics on hybrid high performance computers-particle-particle particle-mesh[J]. Computer Physics Communications,2012,183(3):449-459.

[8] Brown W M,Yamada M. Implementing molecular dynamics on hybrid high performance computers—Three-body potentials[J]. Computer Physics Communications,2013,184(12):2785-2793.

[9] Davis S,Loyola C,González F,et al. Las palmeras molecular dynamics:A flexible and modular molecular dynamics code[J]. Computer Physics Communications,2010,181(12):2126-2139.

[10] Jia Hai-peng,Zhang Yun-quan, Long Guo-ping, et al. Research on Laplace image enhancement algorithm optimization based on OpenCL[J]. Computer Science, 2012,39(5):271-277.(in Chinese)

[11] Jiang Li-yuan,Zhang Yun-quan,Long Guo-ping, et al. Parallelism and research on functions with continuously independent data and intensive memory access using OpenCL[J]. Computer Science, 2013,40(3):111-115.(in Chinese)

[12] Pang Xu,Zhang Yun-quan,Long Guo-ping, et al. Research on mean shift algorithm using OpenCL on multiple many-core platforms[J]. Computer Science, 2013,40(3):79-85.(in Chinese)

[13] Khronos OpenCL Working Group. The OpenCL specification,version 1.2,revision 16[S]. 2011.

[14] Daga M,Aji A M,Wu Chun-feng. On the efficacy of a fused CPU+GPU processor (or APU) for parallel computing[C]∥Proc of 2011 Symposium on Application Accelerators in High-Performance Computing (SAAHPC),2011:141-149.

[15] AMD.AMD accelerated parallel processing OpenCL programming guide[Z].2013.

附中文參考文獻(xiàn)：

[2] 李伯楊,聶峰光,李曉霞,等. GPU 并行計(jì)算集群上的 LAMMPS 分子動(dòng)力學(xué)模擬性能測(cè)試[J]. 計(jì)算機(jī)與應(yīng)用化學(xué),2011,28(10):1229-1233.

[10] 賈海鵬,張?jiān)迫?龍國平,等. 基于 OpenCL 的拉普拉斯圖像增強(qiáng)算法優(yōu)化研究[J]. 計(jì)算機(jī)科學(xué),2012,39(5):271-277.

[11] 蔣麗媛,張?jiān)迫?龍國平,等. 基于 OpenCL 的連續(xù)數(shù)據(jù)無關(guān)訪存密集型函數(shù)并行與優(yōu)化研究[J]. 計(jì)算機(jī)科學(xué),2013,40(3):111-115.

[12] 龐旭,張?jiān)迫?龍國平,等. 基于 OpenCL 的均值平移算法在多個(gè)眾核平臺(tái)的性能優(yōu)化研究[J]. 計(jì)算機(jī)科學(xué),2013,40(3):79-85.

趙成龍(1989-),男，江蘇徐州人，碩士生，研究方向?yàn)橛?jì)算機(jī)系統(tǒng)結(jié)構(gòu)。E-mail:Zhaochenglong08@126.com

ZHAO Cheng-long,born in 1989,MS candidate,his research interest includes computer architecture.

Short-range force algorithm optimization in Lammps based on OpenCL

ZHAO Cheng-long1，SHI Hui-bin1，YU Xin-feng2

(1.College of Computer Science and Technology,Nanjing University of Aeronautics and Astronautics,Nanjing 210016;2.College of Computer Science and Technology,Nanjing University of Science and Technology,Nanjing 210094,China)

Lammps is a widely used open source software for molecular dynamics simulations and related issues, through which we can understand solid and liquid nature. It enables the CUDA and the OpenCL to accelerate calculation. Since the OpenCL has cross-platform features, it becomes a research focus. We summarize the OpenCL kernel programming principles and propose an improved Amdahl's law which can measure the acceleration performance on heterogeneous platforms. The short-range force algorithm in Lammps on Y485P platform is tested. We achieve better acceleration effect by optimizing the kernel source codes of the key parts, such as neighbor-list building and short-range force calculation.

GPU;Lammps;OpenCL;molecular dynamics;short-range force;kernel optimization

1007-130X(2015)09-1614-07

2014-08-25;

2014-10-31

TP302

A

10.3969/j.issn.1007-130X.2015.09.002

通信地址：210016 江蘇省南京市秦淮區(qū)御道街29號(hào)南京航空航天大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院

Address:College of Computer Science and Technology,Nanjing University of Aeronautics and Astronautics,29 Yudao St,Qinhuai District,Nanjing 210016,Jiangsu,P.R.China