彭 龍,陳俊仕,安 虹
(中國(guó)科學(xué)技術(shù)大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,合肥 230031)
我國(guó)自主研發(fā)的神威太湖之光超級(jí)計(jì)算機(jī)[1]共配置40 960個(gè)SW26010異構(gòu)眾核處理器[2],內(nèi)部包括1 040萬(wàn)個(gè)計(jì)算核心,整機(jī)性能位居世界前列。SW26010異構(gòu)眾核處理器由我國(guó)自主設(shè)計(jì),采用主從硬件架構(gòu),共有4個(gè)核組,每個(gè)核組包含1個(gè)控制核心和64個(gè)從核組成的加速陣列,峰值性能高達(dá)每秒12.5億億次。由于神威太湖之光超級(jí)計(jì)算機(jī)特殊的硬件架構(gòu)和編譯環(huán)境,現(xiàn)有商用平臺(tái)上的軟件需要經(jīng)過(guò)重新實(shí)現(xiàn)與優(yōu)化才能充分利用其計(jì)算資源。AMBER[3-4]是主流分子動(dòng)力學(xué)(Molecular Dynamics,MD)模擬軟件之一,能夠有效描述蛋白質(zhì)、核酸以及藥物小分子之間的相互作用,軟件自帶勢(shì)能模型,具有一系列與分子動(dòng)力學(xué)模擬軟件相關(guān)的工具鏈。由于每個(gè)分子的相互作用都由分子力場(chǎng)描述,而AMBER包含60多個(gè)程序,根據(jù)模擬需要協(xié)調(diào)使用,因此其從AMBER11開(kāi)始支持GPU平臺(tái),并利用GPU加速卡加速軟件非鍵力的計(jì)算,同時(shí)支持Intel的MIC平臺(tái)。目前,AMBER在神威太湖之光超級(jí)計(jì)算機(jī)上的研究工作還處于空白狀態(tài),并且國(guó)內(nèi)對(duì)AMBER模擬軟件的并行化研究較少。中國(guó)科學(xué)院上海藥物研究所在開(kāi)展老藥新制的過(guò)程中,使用AMBER軟件進(jìn)行藥物分子的動(dòng)力學(xué)模擬,為能更快地完成藥物分子的模擬計(jì)算,縮短新藥研制周期,對(duì)AMBER軟件性能的優(yōu)化成為亟待解決的問(wèn)題。
本文基于神威太湖之光國(guó)產(chǎn)超級(jí)計(jì)算機(jī),開(kāi)展AMBER分子動(dòng)力學(xué)模擬軟件的移植與優(yōu)化研究。通過(guò)調(diào)整AMBER軟件中熱點(diǎn)函數(shù)的數(shù)據(jù)結(jié)構(gòu),建立主從加速模型,并利用SW26010處理器從核陣列加速熱點(diǎn)函數(shù)計(jì)算。同時(shí),結(jié)合SW26010眾核處理器的硬件特性,提出主從異步流水方案,使主核與從核協(xié)同計(jì)算,解決從核陣列訪存帶寬受限和從核訪存速度過(guò)慢的問(wèn)題,以提升AMBER軟件性能。
分子動(dòng)力學(xué)模擬[5-6]通過(guò)計(jì)算分子體系結(jié)構(gòu)中原子的受力情況,利用經(jīng)典牛頓力學(xué)或量子力學(xué)更新原子位置,從而得到原子運(yùn)行軌跡的模擬過(guò)程。通過(guò)指定模擬步數(shù)不斷進(jìn)行迭代計(jì)算,而單位迭代的時(shí)間步長(zhǎng)通常在飛秒尺度上。一般科學(xué)研究或者實(shí)驗(yàn)需要在微秒甚至更長(zhǎng)的時(shí)間跨度上進(jìn)行分析,因此需要上千萬(wàn)甚至上億次的迭代計(jì)算,而且在每一次迭代過(guò)程中需要計(jì)算系統(tǒng)中每個(gè)原子的受力情況,計(jì)算量非常龐大[7-9]。在整個(gè)流程中需要先進(jìn)行系統(tǒng)原子數(shù)據(jù)初始化,然后根據(jù)各原子位置等信息計(jì)算系統(tǒng)勢(shì)能,從而求出每個(gè)原子的受力情況,再由牛頓力學(xué)求出加速度進(jìn)行速度和坐標(biāo)更新,最后進(jìn)行下一次迭代計(jì)算。AMBER力場(chǎng)勢(shì)能計(jì)算如式(1)所示:
(1)
在式(1)中:等式右邊的前3項(xiàng)表示通過(guò)鍵力項(xiàng)求鍵力勢(shì)能,分別為鍵的伸縮項(xiàng)、鍵角的彎曲項(xiàng)和二面角項(xiàng),其中,Kr、req、Kθ、θeq、Vn和γ為力場(chǎng)相關(guān)參數(shù),r、θ和φ為鍵長(zhǎng)、鍵角和二面角;等式右邊的后2項(xiàng)表示范德華作用項(xiàng)和庫(kù)倫作用項(xiàng),代表非鍵力項(xiàng),其中,Rij為原子對(duì)之間的距離,qi和qj為原子i、j攜帶的電荷量,ε為介電常數(shù)。此處范德華作用項(xiàng)采用Lennard-Jones勢(shì)(簡(jiǎn)稱(chēng)為L(zhǎng)-J勢(shì)),又稱(chēng)12-6勢(shì)能[10],當(dāng)Rij非常大時(shí),該勢(shì)能趨近于0,因此在求得各原子i在該位置上的勢(shì)能Utotal后,通過(guò)對(duì)該勢(shì)能求梯度i得到原子i的受力情況,如式(2)所示:
Fi=-iUtotal
(2)
由牛頓第二定律得到原子i的加速度,如式(3)所示:
(3)
根據(jù)加速度al及原子初始位置和速度,在經(jīng)過(guò)單位時(shí)間步長(zhǎng)t下求出原子速度和位置(如式(4)和式(5)所示),更新速度和坐標(biāo)后可以進(jìn)行下一次的迭代計(jì)算。
(4)
(5)
經(jīng)過(guò)不斷循環(huán)迭代計(jì)算,直到滿足迭代終止條件時(shí)為止。在不斷迭代計(jì)算的過(guò)程中,能得到每個(gè)單位時(shí)間步長(zhǎng)下所有原子運(yùn)動(dòng)的勢(shì)能、位置和速度等信息。
神威太湖之光計(jì)算機(jī)系統(tǒng)是一臺(tái)10億億次量級(jí)的超大規(guī)模并行處理計(jì)算機(jī)系統(tǒng)。神威太湖之光的整體系統(tǒng)結(jié)構(gòu)如圖1所示,其中整機(jī)包含40×4個(gè)超級(jí)節(jié)點(diǎn),每個(gè)超級(jí)節(jié)點(diǎn)包含64個(gè)插件板,每個(gè)插件板包含4個(gè)計(jì)算節(jié)點(diǎn),超級(jí)節(jié)點(diǎn)內(nèi)部采用全連接的模式,可實(shí)現(xiàn)高效的消息廣播。每個(gè)計(jì)算節(jié)點(diǎn)內(nèi)集成1個(gè)SW26010處理器,每個(gè)SW26010處理器包含256個(gè)核,因此共有40 960個(gè)SW26010處理器。SW26010處理器由我國(guó)自主研發(fā),采用片上計(jì)算陣列從核和分布式共享存儲(chǔ)相結(jié)合的異構(gòu)眾核系統(tǒng)結(jié)構(gòu)。
圖1 神威太湖之光系統(tǒng)結(jié)構(gòu)
如圖2所示,1個(gè)處理器由4個(gè)核組組成,每個(gè)核組包含1個(gè)主核和64個(gè)從核,核組間支持緩存(Cache)一致性。主核和從核的工作頻率為1.5 GHz,其中主核擁有32 KB的一級(jí)指令和數(shù)據(jù)Cache以及512 KB的二級(jí)Cache,從核擁有16 KB的指令Cache和64 KB的可重構(gòu)局部數(shù)據(jù)存儲(chǔ)。SW26010集成了4路128位DDR3存儲(chǔ)控制器、8通道PCIe3.0和千兆以太網(wǎng)接口。
圖2 SW26010眾核處理器架構(gòu)
AMBER是由加利福尼亞大學(xué)的KOLLMAN教授專(zhuān)門(mén)針對(duì)生物分子系統(tǒng)模擬而開(kāi)發(fā)的分子動(dòng)力學(xué)應(yīng)用軟件。從1975年到現(xiàn)在,最新版本為AMBER2019。AMBER軟件由60多個(gè)程序組成,這些程序可以分類(lèi)成模擬前的處理程序、模擬程序和模擬后的處理分析程序。此外,AMBER還包含多種分子的勢(shì)能力場(chǎng),用來(lái)描述大分子或者小分子之間的相互作用。
AMBER軟件從AMBER2011開(kāi)始支持GPU加速[11-13],由LEGRAND等人[14-16]共同開(kāi)發(fā)和維護(hù),從AMBER2014開(kāi)始支持Intel Xeno Phi架構(gòu)[17-19],由BHUIYAN等人基于MPI和OpenMP雙重分解策略進(jìn)行研發(fā),將計(jì)算熱點(diǎn)部分放置于Intel協(xié)處理器中進(jìn)行加速,從而提升AMBER軟件性能[20-21]。AMBER2018的pmemd版本增強(qiáng)了整體性能,相比基于Pascal架構(gòu)的AMBER2016性能提升25%~42%。截至2018年10月,AMBER2018已提供對(duì)Turing架構(gòu)(RTX2060、RTX2070、RTX2080和RTX2080TI)的支持。
在移植與優(yōu)化AMBER軟件前需先對(duì)其進(jìn)行測(cè)試與分析,找出sander模擬程序的熱點(diǎn)計(jì)算部分[22]。本文通過(guò)在Intel商用平臺(tái)上模擬二氫葉酸還原酶(DHFR)得到熱點(diǎn)函數(shù)執(zhí)行時(shí)間占比如表1所示。
表1 AMBER_sander模擬程序中熱點(diǎn)函數(shù)計(jì)算時(shí)間占比Table 1 The proportion of hotspot function calculation time in AMBER_sander simulation program
可以看出,sander模擬程序中模擬計(jì)算時(shí)間runmd占整體計(jì)算時(shí)間的92.05%,計(jì)算系統(tǒng)原子之間非鍵力的nonbond熱點(diǎn)函數(shù)占模擬計(jì)算時(shí)間runmd的93.11%,計(jì)算短程非鍵力的get_nb_energy熱點(diǎn)函數(shù)占非鍵力計(jì)算時(shí)間的77.06%,該函數(shù)調(diào)用short_ene函數(shù)來(lái)計(jì)算原子與其鄰居原子的短程非鍵力。因此,本文先將AMBER整個(gè)程序集移植到SW26010主核上,然后利用SW26010處理器的主從硬件特性來(lái)加速AMBER模擬軟件中短程非鍵力的計(jì)算部分,從而提升軟件整體性能。
本文移植工作是將AMBER軟件及其在編譯過(guò)程中依賴(lài)的第三方庫(kù)均移植到神威太湖之光超級(jí)計(jì)算機(jī)的主核上。SW26010主核能夠獨(dú)立執(zhí)行完整的進(jìn)程,而從核陣列作為加速核心,只能執(zhí)行由主核進(jìn)程創(chuàng)建的加速線程。因此,在將熱點(diǎn)函數(shù)的計(jì)算任務(wù)加載到從核前,需要先將AMBER移植到主核上,具體的編譯步驟為:1)修改configure,重新設(shè)置編譯器及相關(guān)參數(shù),配置相關(guān)庫(kù);2)設(shè)置并運(yùn)行configure參數(shù),生成環(huán)境配置腳本和configure頭文件;3)執(zhí)行Makefile完成編譯安裝過(guò)程。由于SW26010處理器的架構(gòu)和編譯環(huán)境與商用CPU有較大的區(qū)別,因此在configure時(shí)需要對(duì)代碼及若干參數(shù)進(jìn)行調(diào)整。此外,AMBER在編譯過(guò)程中需要NETCDF、FFTW等第三方庫(kù)的支持,因此可提前將這些第三方庫(kù)進(jìn)行交叉編譯移植到神威太湖之光超級(jí)計(jì)算機(jī)的主核上,以供AMBER編譯使用。
本文利用從核對(duì)get_nb_energy熱點(diǎn)函數(shù)的計(jì)算過(guò)程進(jìn)行并行實(shí)現(xiàn),首先應(yīng)用OpenACC將熱點(diǎn)函數(shù)從核化,然后采用Athread線程庫(kù)設(shè)計(jì)主從加速和主從異步流水方案,最后通過(guò)從核訪存優(yōu)化進(jìn)一步加速AMBER模擬軟件熱點(diǎn)函數(shù)的計(jì)算。
通過(guò)對(duì)AMBER模擬程序sander熱點(diǎn)計(jì)算部分的分析可知,計(jì)算短程非鍵力的get_nb_energy熱點(diǎn)函數(shù)需要并行加速熱點(diǎn)部分,其中g(shù)et_nb_energy熱點(diǎn)函數(shù)會(huì)在循環(huán)中調(diào)用short_ene函數(shù),然后計(jì)算每個(gè)系統(tǒng)原子i受到的短程非鍵力,具體算法如下:
算法1短程非鍵力計(jì)算算法
輸入系統(tǒng)原子及其鄰居原子的參數(shù)信息
輸出系統(tǒng)原子短程非鍵力
1.do cell index = myindexlo to myindexhi
2.計(jì)算原子索引ncell_lo and ncell_hi
3.do i = ncell_lo to ncell_hi
4.計(jì)算原子i的鄰居原子個(gè)數(shù)ntot;
5.對(duì)ntot為原子i的所有鄰居原子分配存儲(chǔ)中間結(jié)果的內(nèi)存空間;
6.調(diào)用shor_ene函數(shù)計(jì)算原子i對(duì)其所有鄰居原子的非鍵力作用;
7.釋放為原子i的鄰居原子分配的內(nèi)存空間;
8.更新鄰居列表ipairs的偏移指針numpack
9.end do
10.end do
如算法1所示,整個(gè)計(jì)算過(guò)程由兩重循環(huán)組成。在外層循環(huán)中,所有網(wǎng)格單元index通過(guò)內(nèi)層循環(huán)遍歷每個(gè)網(wǎng)格單元中的所有原子i,并對(duì)遍歷到的每個(gè)原子i調(diào)用short_ene函數(shù)計(jì)算其與所有鄰居原子的非鍵力,同時(shí)釋放已為原子i分配的內(nèi)存空間且更新鄰居列表ipairs的偏移指針numpack。在內(nèi)層循環(huán)中,根據(jù)原子i的鄰居列表大小為該原子申請(qǐng)臨時(shí)內(nèi)存空間存放中間數(shù)據(jù),為減少內(nèi)存消耗,AMBER軟件自身維護(hù)一個(gè)stack模塊,用于在計(jì)算原子i與其鄰居原子的非鍵力前動(dòng)態(tài)分配內(nèi)存,在數(shù)據(jù)更新后釋放已分配的內(nèi)存空間。在申威編譯器支持下的Fortran語(yǔ)言中,使用SW26010處理器從核并行時(shí)需要將核心計(jì)算相關(guān)的數(shù)據(jù)通過(guò)Common區(qū)進(jìn)行共享,即主核和從核均可以直接訪問(wèn)和修改共享存儲(chǔ)空間的數(shù)據(jù)。然而,通過(guò)該方式實(shí)現(xiàn)數(shù)據(jù)共享需要確保Common區(qū)的數(shù)組大小固定,而AMBER采用動(dòng)態(tài)內(nèi)存分配方式,根據(jù)原子參數(shù)的不同分配不同大小的內(nèi)存空間,顯然該內(nèi)存空間分配方式并不適用于從核并行化。
另外,在計(jì)算原子非鍵力時(shí),由于某些原子的鄰居原子存在重疊,如果使用從核并行計(jì)算就會(huì)產(chǎn)生寫(xiě)寫(xiě)數(shù)據(jù)依賴(lài),導(dǎo)致錯(cuò)誤的計(jì)算結(jié)果,因此在并行計(jì)算前需要解決動(dòng)態(tài)內(nèi)存和寫(xiě)寫(xiě)依賴(lài)的問(wèn)題。本文通過(guò)調(diào)整數(shù)據(jù)結(jié)構(gòu)的方式來(lái)解決這兩個(gè)問(wèn)題,首先將stack動(dòng)態(tài)內(nèi)存分配方式調(diào)整為固定內(nèi)存分配方式,將分配的內(nèi)存設(shè)置為Common區(qū)以供從核并行計(jì)算訪問(wèn),然后分配一個(gè)全局force_sw數(shù)組用來(lái)存放并行計(jì)算時(shí)所有鄰居原子的非鍵力信息,在從核計(jì)算結(jié)束后,通過(guò)主核串行的方式更新到全局force數(shù)組中,從而解決相同鄰居原子的寫(xiě)寫(xiě)依賴(lài)問(wèn)題。采用固定內(nèi)存分配方式的主存Common區(qū)和Global區(qū)數(shù)據(jù)分配如圖3所示。
圖3 主存數(shù)據(jù)分配
由圖3可以看出,在計(jì)算非鍵力前,首先將AMBER自身維護(hù)的動(dòng)態(tài)內(nèi)存分配方式改為固定內(nèi)存分配方式,分為存放常量參數(shù)的Constant區(qū)(包含用于輸入的3個(gè)整型(integer)數(shù)據(jù)和3個(gè)實(shí)型(real)數(shù)據(jù))和存放系統(tǒng)原子參數(shù)的System區(qū)以及存放鄰居原子參數(shù)的Neighbor區(qū)(包含系統(tǒng)總原子個(gè)數(shù)(num_atoms)和ipairs鄰居列表長(zhǎng)度(lastatom)的若干integer和real數(shù)據(jù),用于相關(guān)輸入和輸出),這些內(nèi)存空間均被標(biāo)識(shí)為Common,以供從核計(jì)算訪問(wèn)。在從核并行計(jì)算前,由主核訪問(wèn)主存中的Global區(qū)從而對(duì)input標(biāo)記的Common區(qū)的部分內(nèi)存進(jìn)行數(shù)據(jù)初始化,Global區(qū)中包含存儲(chǔ)最后結(jié)果的force數(shù)組、鄰居列表ipairs及其他相關(guān)參數(shù)信息。在從核并行計(jì)算完成后,從核將計(jì)算結(jié)果更新到output標(biāo)記的Common區(qū)。在調(diào)整算法中的數(shù)據(jù)結(jié)構(gòu)后,將進(jìn)行熱點(diǎn)計(jì)算部分的從核并行化。
在調(diào)整熱點(diǎn)函數(shù)的數(shù)據(jù)結(jié)構(gòu)后,使用OpenACC實(shí)現(xiàn)主核與從核協(xié)作計(jì)算的執(zhí)行模型,從而充分利用SW26010的主從硬件資源加速熱點(diǎn)函數(shù)?;贠penACC編寫(xiě)程序?qū)?shù)據(jù)管理和傳輸設(shè)置為隱式,由編譯器根據(jù)編譯指示信息,自動(dòng)生成空間管理、數(shù)據(jù)傳輸?shù)瓤刂拼a,并且通過(guò)調(diào)整數(shù)據(jù)結(jié)構(gòu)解決熱點(diǎn)函數(shù)并行時(shí)產(chǎn)生的寫(xiě)寫(xiě)依賴(lài)問(wèn)題后并行計(jì)算原子之間的短程非鍵力,OpenACC實(shí)現(xiàn)算法具體如下:
算法2OpenACC并行化短程非鍵力計(jì)算算法
輸入Common區(qū)數(shù)據(jù)參數(shù)
輸出系統(tǒng)原子短程非鍵力
1.主核進(jìn)行Common區(qū)數(shù)據(jù)初始化
2.!$ACC parallel loop/local/copyin/copyout/reduction/annotate(參數(shù)變量)//從核化
3.do atom i = gstart to gend
4.!$ACC data copyin/copy//加速數(shù)據(jù)區(qū)指示
5.從Common區(qū)/LDM中獲取原子數(shù)據(jù)
6.從核并行計(jì)算ki與其鄰居原子的非鍵力
7.Update force_sw data
8.!$ACC end data
9.end do
10.!$ACC END parallel loop//從核計(jì)算結(jié)束
11.Update nonbond force data
將整個(gè)系統(tǒng)原子的短程非鍵力計(jì)算循環(huán)外加OpenACC指示,當(dāng)主核運(yùn)行到該函數(shù)時(shí),按照相關(guān)指示將整個(gè)循環(huán)交給從核并行完成,從核完成計(jì)算后,按照相關(guān)指示將計(jì)算結(jié)果拷貝或者歸約到主存Common區(qū)。在從核局部數(shù)據(jù)緩存(Local Data Memory,LDM)中為循環(huán)中的一些變量分配存儲(chǔ)空間,主要是從核計(jì)算過(guò)程中需要用到的存儲(chǔ)中間數(shù)據(jù)的變量,這些變量要足夠小且在從核計(jì)算過(guò)程中被頻繁訪問(wèn),如原子電荷、相對(duì)距離等,通過(guò)local指示在LDM中分配存儲(chǔ)空間。
local中的部分變量需要作為結(jié)果在計(jì)算完成后拷貝到主存中,這些變量通過(guò)copyout指示。copy和copyin指示變量在計(jì)算前由主存讀入LDM中,在默認(rèn)情況下根據(jù)循環(huán)劃分方式將對(duì)應(yīng)數(shù)據(jù)進(jìn)行分割拷貝到LDM中,主要為原子索引數(shù)組、存儲(chǔ)結(jié)果的非鍵力數(shù)組以及原子坐標(biāo)數(shù)組等全局?jǐn)?shù)組。copy指示的變量在計(jì)算結(jié)束后需要將LDM中的值拷貝到主存中,而copyin指示的變量在計(jì)算結(jié)束后并不將值拷貝到主存中。在循環(huán)中通過(guò)加速數(shù)據(jù)區(qū)指示data定義加速計(jì)算的數(shù)據(jù)區(qū),用來(lái)描述變量在LDM內(nèi)的屬性,指示在計(jì)算完成后是否需要將值從LDM中更新回主存,其中使用reduction規(guī)約子句完成部分變量在各從核計(jì)算結(jié)束后的加法規(guī)約操作,具體代碼如下:
reduction(+:vxx_si,eelt,evdw,vzz_si,vyz_si,vyy_si,vxz_si,vxy_si,ee_vir_iso_si)
編譯器會(huì)在每個(gè)線程上為reduction指示的變量創(chuàng)建一個(gè)私有副本,并根據(jù)加法操作以及變量數(shù)據(jù)類(lèi)型對(duì)該私有副本進(jìn)行初始化。在循環(huán)結(jié)束后,使用加法操作對(duì)同一變量的所有私有副本進(jìn)行一次并行規(guī)約的加法操作,并將規(guī)約結(jié)果寫(xiě)回原始變量中。至此,根據(jù)循環(huán)中的變量及其計(jì)算特點(diǎn)通過(guò)添加OpenACC編譯指示使得整個(gè)循環(huán)計(jì)算在64個(gè)從核中并行執(zhí)行,完成熱點(diǎn)函數(shù)在從核上的并行計(jì)算。
本文通過(guò)調(diào)用神威定制的Athread加速線程庫(kù)函數(shù)接口編寫(xiě)從核代碼,利用主從加速模型加速AMBER軟件的熱點(diǎn)函數(shù)計(jì)算。根據(jù)SW26010處理器特殊的主從硬件架構(gòu),并將從核作為加速核心,用于加速AMBER軟件熱點(diǎn)函數(shù)的計(jì)算,通過(guò)64個(gè)從核的并行計(jì)算提升軟件整體性能,主從加速模型的計(jì)算流程如圖4所示。
圖4 熱點(diǎn)函數(shù)的主從加速模型
主核在完成Common區(qū)內(nèi)存分配后,循環(huán)遍歷系統(tǒng)中所有原子,由原子索引號(hào)訪問(wèn)Global區(qū)獲得原子及其鄰居原子的參數(shù)信息,進(jìn)而對(duì)Common區(qū)的部分內(nèi)存進(jìn)行數(shù)據(jù)初始化。Common區(qū)常量參數(shù)部分包括系統(tǒng)原子的索引號(hào)、第一次分配的標(biāo)志號(hào)以及短程力的截?cái)喟霃降?。系統(tǒng)原子參數(shù)部分包括該原子的鄰居原子個(gè)數(shù)、電荷量以及在全局鄰居列表中的偏移量等。鄰居原子參數(shù)部分包括原子間的相對(duì)距離、evdw參數(shù)和用于存儲(chǔ)中間計(jì)算結(jié)果的force數(shù)組等。主核在完成Common區(qū)數(shù)據(jù)初始化工作后,會(huì)通過(guò)Athread接口函數(shù)來(lái)調(diào)用從核代碼,64個(gè)從核通過(guò)訪問(wèn)Common區(qū)分配的數(shù)據(jù),并行計(jì)算原子的短程非鍵力。
本文將系統(tǒng)原子均勻地分配給64個(gè)從核進(jìn)行計(jì)算,保證每個(gè)從核的計(jì)算量相同,這樣在若干計(jì)算周期后所有從核能同時(shí)完成數(shù)據(jù)計(jì)算。在從核計(jì)算過(guò)程中主核處于阻塞狀態(tài),等待全部從核完成計(jì)算后,主核開(kāi)始進(jìn)行最終數(shù)據(jù)的更新,此時(shí)從核處于空閑狀態(tài),通過(guò)訪問(wèn)Common區(qū)的output部分更新Global區(qū)的force數(shù)組。至此,利用Athread完成熱點(diǎn)函數(shù)計(jì)算的從核化。在主從加速并行方案下,從核計(jì)算過(guò)程中的主核處于阻塞狀態(tài),需要等待從核計(jì)算完成后才能進(jìn)行數(shù)據(jù)更新。另外,主核在從核并行計(jì)算前需要循環(huán)遍歷系統(tǒng)原子進(jìn)行Common區(qū)數(shù)據(jù)初始化,增加了額外的時(shí)間開(kāi)銷(xiāo)。因此,本文將利用主從異步并行方案,最大限度地減少額外的時(shí)間開(kāi)銷(xiāo),實(shí)現(xiàn)熱點(diǎn)函數(shù)的性能加速。
本文采用主從異步并行方案實(shí)現(xiàn)主核和從核并行計(jì)算,通過(guò)流水化設(shè)計(jì)隱藏主核Common區(qū)數(shù)據(jù)初始化帶來(lái)的額外時(shí)間開(kāi)銷(xiāo)。該流水化方案是在從核進(jìn)行并行計(jì)算的同時(shí),使主核不等待,而是進(jìn)行數(shù)據(jù)初始化。流水化前后的主核算法具體如下:
算法3流水化前的主核算法
功能主從加速方案
1.內(nèi)存分配allocate(參數(shù))
2.do i = gstart,gend
3.查找Global區(qū)
4.Common區(qū)數(shù)據(jù)初始化
5.end do
6.call athread_spawn(sw_fun,1)//調(diào)用從核代碼
7.call athread_join()//從核計(jì)算結(jié)束
8.全局?jǐn)?shù)據(jù)更新force
算法4流水化后的主核算法
功能主從異步方案
1.內(nèi)存分配allocate(參數(shù))
2.Call athread_spawn(sw_fun,1)//調(diào)用從核代碼
3.do i=gstart,gend//主核循環(huán)遍歷系統(tǒng)原子
4.查找Global區(qū)
5.Common區(qū)數(shù)據(jù)初始化
6.標(biāo)志Flage()更新//設(shè)置全局標(biāo)記
7.end do
8.call athread_join()//從核計(jì)算結(jié)束
9.全局?jǐn)?shù)據(jù)更新force//主核進(jìn)行全局?jǐn)?shù)據(jù)更新
在流水化前的主核算法中,由主核分配Common區(qū)并對(duì)其進(jìn)行數(shù)據(jù)初始化,然后調(diào)用從核代碼執(zhí)行從核計(jì)算。在流水化后的主核算法中,在主核進(jìn)行數(shù)據(jù)初始化前,使用athread_spawn函數(shù)接口創(chuàng)建從核線程組開(kāi)始執(zhí)行從核代碼。與此同時(shí),主核進(jìn)行主存Common區(qū)的數(shù)據(jù)初始化,接著調(diào)用athread_join等待所有從核完成短程非鍵力的計(jì)算,最后由主核進(jìn)行數(shù)據(jù)更新。為實(shí)現(xiàn)主核數(shù)據(jù)初始化與從核計(jì)算的并行化,本文設(shè)置一個(gè)Flage全局?jǐn)?shù)組作為標(biāo)志放在Common區(qū),用來(lái)判斷為原子i分配的Common區(qū)數(shù)據(jù)是否初始化完成。主核循環(huán)遍歷系統(tǒng)所有原子,通過(guò)原子索引訪問(wèn)主存Global區(qū)為原子i的Common區(qū)初始化,同時(shí)從核檢測(cè)到原子i-1初始化完成后開(kāi)始計(jì)算短程非鍵力,具體流程如圖5所示。
圖5 主從流水化并行執(zhí)行
當(dāng)主核為原子i的Commoni區(qū)初始化完成后,該區(qū)標(biāo)志為真,從核檢測(cè)到標(biāo)志為真即可開(kāi)始計(jì)算Computei,而當(dāng)從核計(jì)算Computei時(shí),主核進(jìn)行Common(i+1)區(qū)的數(shù)據(jù)初始化。因此,主核Common(i+1)區(qū)數(shù)據(jù)初始化可與從核計(jì)算Computei并行執(zhí)行,實(shí)現(xiàn)主核數(shù)據(jù)初始化與從核計(jì)算兩個(gè)模塊的流水化操作,從而消除部分等待時(shí)間。
在進(jìn)行主從異步流水化設(shè)計(jì)后,使主核和從核并行工作,從而消除主從加速方案導(dǎo)致的額外時(shí)間開(kāi)銷(xiāo),進(jìn)一步提升熱點(diǎn)函數(shù)的性能加速。下文將針對(duì)從核訪存進(jìn)行測(cè)試與分析,由于AMBER_sander模擬程序在模擬分子時(shí),根據(jù)模擬步數(shù)不斷迭代計(jì)算,因此需要不斷地從內(nèi)存中讀取和更新數(shù)據(jù)。在從核并行計(jì)算時(shí),64個(gè)從核并行訪存主存,這樣會(huì)導(dǎo)致主存帶寬壓力過(guò)大,而且在上述實(shí)現(xiàn)中,從核通過(guò)Global方式直接訪問(wèn)主存,訪存速度很慢,因此本文針對(duì)訪存帶寬和訪存速度進(jìn)行以下優(yōu)化:
1)針對(duì)從核訪存速度過(guò)慢的問(wèn)題,本文選擇將從核計(jì)算中需要用到的部分?jǐn)?shù)據(jù)暫存在從核LDM中。首先這些數(shù)據(jù)要在從核計(jì)算過(guò)程中被頻繁訪問(wèn),從核訪問(wèn)LDM要比直接訪問(wèn)主存快很多,如果被頻繁訪問(wèn)的數(shù)據(jù)放到從核LDM中,則在從核計(jì)算過(guò)程中可以直接從LDM中讀取數(shù)據(jù),無(wú)需多次訪問(wèn)主存。其次這些數(shù)據(jù)要足夠小,從核LDM只有64 KB,不足以存儲(chǔ)從核計(jì)算過(guò)程中的所有數(shù)據(jù),因此選擇數(shù)據(jù)量較小且在計(jì)算過(guò)程中被頻繁訪問(wèn)的數(shù)據(jù)存放到從核LDM中,從而加速訪存時(shí)間,提升程序性能。這些數(shù)據(jù)主要包括原子電荷(cgi和cgj)和系統(tǒng)原子坐標(biāo)參數(shù)(Xk、Yk和Zk),其中電荷以real類(lèi)型、坐標(biāo)以integer類(lèi)型存儲(chǔ)在Common區(qū),在計(jì)算短程非鍵力的過(guò)程中被訪問(wèn)的次數(shù)與系統(tǒng)原子總個(gè)數(shù)(atoms)成正比,如表2所示。
表2 存放在從核LDM中的數(shù)據(jù)Table 2 Data stored in the LDM of slave core
2)針對(duì)從核并行訪存的帶寬競(jìng)爭(zhēng)問(wèn)題,采用直接內(nèi)存存取(Direct Memory Access,DMA)通道的方式實(shí)現(xiàn)主存與從核LDM間的數(shù)據(jù)交換。一方面,由于DMA可以控制主從與從核LDM之間的數(shù)據(jù)交換,在數(shù)據(jù)交換的同時(shí)從核可以分配一些與所傳輸數(shù)據(jù)無(wú)關(guān)的計(jì)算操作,實(shí)現(xiàn)從核通信與計(jì)算的重疊。另一方面,由于AMBER_sander在求解短程非鍵力時(shí)需依次遍歷系統(tǒng)中所有的原子,數(shù)據(jù)空間局部性大,而在DMA數(shù)據(jù)傳輸時(shí)是將一塊連續(xù)數(shù)據(jù)段從主存讀取到LDM中,在數(shù)據(jù)更新時(shí)也是更新一塊連續(xù)的數(shù)據(jù)段到主存中,因此通過(guò)DMA的異步性和連續(xù)性可以加快從核訪問(wèn)主存的速度。
神威太湖之光超級(jí)計(jì)算機(jī)支持256 bit向量指令,1次可以加載4個(gè)double類(lèi)型。為進(jìn)一步提升熱點(diǎn)函數(shù)的加速性能,本文使用向量指令集對(duì)從核計(jì)算過(guò)程中的部分無(wú)依賴(lài)且規(guī)整的代碼進(jìn)行手動(dòng)向量化操作。在計(jì)算從核的原子非鍵力作用后,得到原子在x、y和z方向上的作用分量dumx、dumy和dumz,利用這3個(gè)部分的作用分量更新臨時(shí)變量force_sw。這部分計(jì)算是一個(gè)3×3的雙重循環(huán),為對(duì)其進(jìn)行向量化操作加速,需要為dum和force_sw數(shù)組添加一個(gè)維度,該維度的值沒(méi)有意義,其不進(jìn)行主核的全局force數(shù)組更新,但參與從核計(jì)算。通過(guò)添加一個(gè)維度將該迭代計(jì)算變?yōu)?×4的兩層for循環(huán),然后使用4個(gè)向量取代兩層for循環(huán)?;谙蛄炕瘜?shí)現(xiàn),從核的執(zhí)行效率得到進(jìn)一步提升,其中參與向量化的部分從核代碼的運(yùn)行速度在理論上能達(dá)到4倍的線性加速,對(duì)于整個(gè)程序的計(jì)算性能也有所提升。
本文完成AMBER軟件在神威太湖之光平臺(tái)上的并行實(shí)現(xiàn)與性能優(yōu)化工作。本節(jié)將針對(duì)神威太湖之光平臺(tái)上移植與優(yōu)化后的AMBER軟件與Intel商用平臺(tái)上的AMBER軟件做性能測(cè)試對(duì)比。測(cè)試平臺(tái)的配置信息如表3和表4所示。
表3 神威太湖之光測(cè)試平臺(tái)配置信息Table 3 Configuration information of Sunway TaihuLighttest platform
表4 Intel商用測(cè)試平臺(tái)的配置信息Table 4 Configuration information of Intel commercial test platform
測(cè)試算例選擇AMBER軟件的DHFR測(cè)試基準(zhǔn),結(jié)構(gòu)如圖6所示,包含23 558個(gè)原子,單位迭代時(shí)間步長(zhǎng)為2 fs,短程非鍵力截?cái)喟霃綖? angstrom。
圖6 DHFR結(jié)構(gòu)
本文測(cè)試了AMBER軟件的熱點(diǎn)函數(shù)在神威太湖之光平臺(tái)的SW26010處理器上采取不同優(yōu)化方案后的執(zhí)行時(shí)間,并對(duì)比其在Intel Xeon Platinum 8163上的執(zhí)行時(shí)間,具體測(cè)試數(shù)據(jù)如表5所示。
表5 AMBER熱點(diǎn)函數(shù)在神威太湖之光與Intel商用平臺(tái)上的測(cè)試數(shù)據(jù)Table 5 Test data of AMBER hotspot function on theSunway TaihuLight and Intel commercial platform s
在神威太湖之光平臺(tái)的SW26010處理器上采取不同優(yōu)化方案后的熱點(diǎn)函數(shù)加速情況,如圖7所示。可以看出,優(yōu)化前AMBER_sander熱點(diǎn)函數(shù)在主核上的執(zhí)行時(shí)間為15.57 s。在實(shí)現(xiàn)熱點(diǎn)函數(shù)主從加速并行方案后,由于未進(jìn)行訪存優(yōu)化且在從核計(jì)算過(guò)程中主核處于阻塞狀態(tài),因此函數(shù)執(zhí)行時(shí)間增加至32.23 s。在實(shí)現(xiàn)主從異步流水并行方案后,由于主核和從核之間存在并行協(xié)作計(jì)算,使得熱點(diǎn)函數(shù)的執(zhí)行時(shí)間下降至19.04 s,然后進(jìn)行訪存優(yōu)化后,熱點(diǎn)函數(shù)的執(zhí)行時(shí)間下降至1.25 s。通過(guò)上述所有優(yōu)化方案后,本文使用SW26010單核組的64個(gè)從核將AMBER_sander的熱點(diǎn)函數(shù)計(jì)算性能約提升了15倍。
圖7 熱點(diǎn)函數(shù)性能加速
由于不同優(yōu)化方案對(duì)AMBER_sander有不同的優(yōu)化效果,因此本文統(tǒng)計(jì)了不同優(yōu)化方案對(duì)熱點(diǎn)函數(shù)性能加速的執(zhí)行時(shí)間占比,測(cè)試數(shù)據(jù)如表6所示。在實(shí)現(xiàn)主從加速并行化后,采用主從異步流水、訪存優(yōu)化(LDM和DMA通道)以及SIMD向量化對(duì)并行化版本的AMBER做性能優(yōu)化,使得AMBRE軟件的熱點(diǎn)函數(shù)執(zhí)行時(shí)間降低至31.20 s(主從+訪存+SIMD)。
表6 優(yōu)化方案對(duì)AMBER_sander的性能加速情況Table 6 Performance acceleration of AMBER_sander by optimization schemes
可以看出,不同優(yōu)化方案對(duì)性能的加速效果不同,如圖8所示,針對(duì)從核訪存(LDM+DMA)的優(yōu)化占了總優(yōu)化時(shí)間(31.20 s)的50%以上,通過(guò)測(cè)試分析可以為后續(xù)在神威太湖之光平臺(tái)上移植與優(yōu)化其他軟件提供借鑒與指導(dǎo)。
圖8 優(yōu)化方案對(duì)AMBER_sander的性能加速占比
本文將AMBER_sander在神威太湖之光平臺(tái)上單核組的單進(jìn)程性能與Intel商用平臺(tái)上單進(jìn)程性能進(jìn)行對(duì)比,如圖9所示。
圖9 AMBER_sander在神威太湖之光與Intel商用平臺(tái)上的性能對(duì)比
可以看出,AMBER_sander熱點(diǎn)函數(shù)在SW26010單核組上優(yōu)化后的性能相比Intel Xeon Platinum 8163上約提升了4.6倍。值得注意的是,優(yōu)化后AMBER_sander熱點(diǎn)函數(shù)在單核組上的整體性能略高于Intel商用平臺(tái)上的性能。這是因?yàn)殡S著優(yōu)化過(guò)程的不斷深入,熱點(diǎn)函數(shù)的執(zhí)行時(shí)間越來(lái)越短,非熱點(diǎn)函數(shù)的執(zhí)行時(shí)間在程序總執(zhí)行時(shí)間中所占的比重越來(lái)越大,根據(jù)Amdahl定律熱點(diǎn)函數(shù)執(zhí)行效率提升所帶來(lái)的整體性能收益也變得越來(lái)越小。
本文基于神威太湖之光處理器的主從硬件結(jié)構(gòu),建立主從加速模型,并實(shí)現(xiàn)AMBER軟件的從核并行化設(shè)計(jì)。在并行化版本的基礎(chǔ)上,提出主從異步流水化方案,針對(duì)從核的并行訪存進(jìn)行一系列性能優(yōu)化,極大提升AMBER軟件在神威太湖之光平臺(tái)上的計(jì)算性能。應(yīng)用所有優(yōu)化方案后,AMBER熱點(diǎn)函數(shù)性能較優(yōu)化前約提升15倍,其單核組性能相比Intel Xeon Platinum 8163約提升4.6倍。隨著計(jì)算節(jié)點(diǎn)的增多,節(jié)點(diǎn)間的數(shù)據(jù)傳輸將成為性能瓶頸,因此后續(xù)將開(kāi)展AMBER軟件在神威太湖之光平臺(tái)上的強(qiáng)擴(kuò)展性研究,解決節(jié)點(diǎn)間的消息傳輸帶寬限制問(wèn)題。