面向異構(gòu)計(jì)算的高性能計(jì)算算法與軟件?

徐 順 ,王 武 ,張 鑒 ,姜金榮 ,金 鐘,3 ,遲學(xué)斌

1(中國(guó)科學(xué)院 計(jì)算機(jī)網(wǎng)絡(luò)信息中心,北京 100190)

2(中國(guó)科學(xué)院 計(jì)算科學(xué)應(yīng)用研究中心,北京 100190)

3(廣東省生物醫(yī)藥計(jì)算重點(diǎn)實(shí)驗(yàn)室(中國(guó)科學(xué)院 廣州生物醫(yī)藥與健康研究院),廣東 廣州 510530)

高性能計(jì)算應(yīng)用對(duì)國(guó)家工業(yè)、國(guó)家安全、生命健康和科學(xué)發(fā)現(xiàn)均具有重要的作用.作為信息科技政策制定咨詢權(quán)威組織,美國(guó)信息技術(shù)與創(chuàng)新基金會(huì)(Information Technology and Innovation Foundation,簡(jiǎn)稱ITIF)早在2016 年就作了題為《高性能計(jì)算對(duì)美國(guó)競(jìng)爭(zhēng)力至關(guān)重要》[1]的專題分析報(bào)告.近年來,我國(guó)也探索出了“天河”、“神威”和“曙光”三大國(guó)產(chǎn)高性能計(jì)算機(jī)的發(fā)展之路,均具備自主研制高性能處理器和加速器的能力;同時(shí)也計(jì)劃在2020～2021 年前后設(shè)計(jì)實(shí)現(xiàn)E 級(jí)(每秒百億億次浮點(diǎn)運(yùn)算)計(jì)算機(jī),在全球高性能計(jì)算性能Top500 排行榜上[2],可看到國(guó)產(chǎn)超級(jí)計(jì)算機(jī)性能處于全球前列水平.但不難看出,我國(guó)超算系統(tǒng)軟件和應(yīng)用軟件的研發(fā)步調(diào)相對(duì)落后,國(guó)產(chǎn)超算上的軟件生態(tài)沒有同步形成[3].一個(gè)客觀的困難是:國(guó)際上主流基于Intel x86 或Nvidia GPU 的科學(xué)計(jì)算軟件無法直接在部分國(guó)產(chǎn)超算上編譯運(yùn)行,特別是面向異構(gòu)計(jì)算加速這方面,需要有大量的移植適配工作.對(duì)于神威超算這種特殊的主從核架構(gòu),大部分軟件需要修改甚至重寫以實(shí)現(xiàn)從核加速計(jì)算.

盡管高性能計(jì)算可應(yīng)用于工業(yè)、國(guó)防、生命健康等諸多領(lǐng)域,但這些應(yīng)用主要以科學(xué)計(jì)算為基礎(chǔ).隨著計(jì)算機(jī)軟硬件技術(shù)的高速發(fā)展,特別是對(duì)稱處理器核數(shù)、內(nèi)存容量、網(wǎng)絡(luò)帶寬的提升引發(fā)(并行)計(jì)算算法需要考慮新的平衡因素,促使在科學(xué)計(jì)算與數(shù)值模擬應(yīng)用方面提出了大量新的軟件架構(gòu)和理論算法.比如在高能物理領(lǐng)域,形成了美國(guó)USQCD[4]和MILC[5]、英國(guó)Grid[6]和日本Bridge++軟件等諸多代表性的格點(diǎn)量子色動(dòng)力學(xué)(quantum chromodynamics,簡(jiǎn)稱QCD)模擬軟件,這些軟件形成了自身特色的軟件架構(gòu),或者是面向特定體系結(jié)構(gòu)而優(yōu)化設(shè)計(jì)的.特別是USQCD 的Chroma 軟件構(gòu)建了一套多層次、松耦合的軟件架構(gòu)體系,取得了有效整合軟件模塊和增強(qiáng)硬件適應(yīng)性的顯著效果;他們也開發(fā)了基于GPU 加速計(jì)算的軟件包QUDA[7],充分利用了GPU異構(gòu)計(jì)算的優(yōu)勢(shì).在分子模擬應(yīng)用領(lǐng)域,LAMMPS[8]和NAMD[9]軟件為代表的應(yīng)用軟件在高性能計(jì)算方面具有顯著的性能優(yōu)勢(shì),但它們的軟件架構(gòu)完成不同,而兩者都具有代表性.LAMMPS 是典型的MPI 分布式并行程序,具有非常強(qiáng)的軟件功能和計(jì)算性能的可擴(kuò)展性;NAMD 更是采用Charm++并行計(jì)算框架支持新的并行編程模式,數(shù)據(jù)通信采用遠(yuǎn)程對(duì)象的函數(shù)調(diào)用方式實(shí)現(xiàn).此外,兼?zhèn)淞孔佑?jì)算化學(xué)研究的NWChem[10]軟件引入了Global Arrays 并行計(jì)算框架來支撐計(jì)算化學(xué)類的特殊問題,更多的是考慮高可伸縮的計(jì)算化學(xué)算法使用傳統(tǒng)的消息傳遞通信模型很繁瑣,而使用可移植共享內(nèi)存編程模型更簡(jiǎn)便[11].以上這些高性能計(jì)算典型應(yīng)用軟件不僅著眼于科學(xué)問題的研究,而且通過其成功的應(yīng)用大大推動(dòng)了高性能計(jì)算應(yīng)用技術(shù)的發(fā)展,在高性能異構(gòu)計(jì)算方面取得了諸多理論和應(yīng)用成果.比如格點(diǎn)QCD 應(yīng)用和NAMD 軟件均獲得過高性能計(jì)算應(yīng)用領(lǐng)域重要的“戈登?貝爾”獎(jiǎng)(Gordon Bell Prize,簡(jiǎn)稱GB 獎(jiǎng)).美國(guó)的E 級(jí)計(jì)算計(jì)劃(Exascale Computing Project,簡(jiǎn)稱ECP)[12]中均毫無懸念地將這些軟件作為重點(diǎn)發(fā)展的對(duì)象.

相形之下,我國(guó)在相關(guān)的科學(xué)計(jì)算應(yīng)用領(lǐng)域的自研程序代碼普遍依托于相應(yīng)領(lǐng)域的課題組開發(fā)和繼承、大多是以科學(xué)問題研究為主,缺乏專業(yè)的開發(fā)人員來維護(hù)軟件.更是缺乏前瞻性的軟件頂層設(shè)計(jì),也難以形成規(guī)?；目赏茝V使用的軟件成果,軟件的并行計(jì)算性能更是難以保證.近年來,GPU 異構(gòu)計(jì)算和國(guó)產(chǎn)自主處理器的興起使得軟件框架和軟件模塊不斷推陳出新,導(dǎo)致軟件程序無法跟上硬件架構(gòu)和軟件應(yīng)用環(huán)境快速發(fā)展的腳步,并行計(jì)算設(shè)計(jì)優(yōu)化加大了困難.慶幸的是,莫?jiǎng)t堯團(tuán)隊(duì)針對(duì)以上問題探索出了一條創(chuàng)新的技術(shù)路徑.其基本方法是凝練數(shù)值模擬應(yīng)用領(lǐng)域的共性需求來建模,并設(shè)計(jì)模型驅(qū)動(dòng)的并行算法和采用了并行編程技術(shù)[13].他們提出了面向結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格的標(biāo)準(zhǔn)化網(wǎng)格數(shù)據(jù)模型,分別開發(fā)了結(jié)構(gòu)網(wǎng)格的編程框架JASMIN[14]和非結(jié)構(gòu)網(wǎng)格的編程框架JAUMIN[15],均一定程度支持并行算法和并行編程的軟件復(fù)用.雖然我國(guó)科學(xué)計(jì)算軟件的開發(fā)積累了不少經(jīng)驗(yàn),但從串行程序到大規(guī)模并行程序、從算法程序到軟件產(chǎn)品還有相當(dāng)長(zhǎng)的道路要走.從《國(guó)家高性能計(jì)算環(huán)境發(fā)展報(bào)告》[16]可得,在我國(guó)各大超級(jí)計(jì)算中心的超級(jí)計(jì)算系統(tǒng)上所運(yùn)行應(yīng)用軟件大部分來自國(guó)外,一旦出現(xiàn)“斷供”,我國(guó)計(jì)算科學(xué)研究即顯現(xiàn)出被“卡脖子”的被動(dòng)局面.

因此,為了發(fā)展國(guó)產(chǎn)超級(jí)計(jì)算系統(tǒng)的應(yīng)用生態(tài)并適應(yīng)未來E 級(jí)計(jì)算的需求,融合國(guó)際主流軟件的優(yōu)點(diǎn)并實(shí)現(xiàn)支撐科學(xué)計(jì)算應(yīng)用研究的最終目的,非常有必要梳理適配國(guó)產(chǎn)高性能計(jì)算系統(tǒng)的科學(xué)計(jì)算應(yīng)用軟件的發(fā)展情況.本文將從高性能計(jì)算應(yīng)用算法和軟件的角度,闡述一些典型應(yīng)用軟件核心算法的并行計(jì)算特征,并討論面對(duì)國(guó)產(chǎn)異構(gòu)系統(tǒng)的算法對(duì)策及其發(fā)展思路.

1 典型高性能計(jì)算應(yīng)用的核心算法與軟件特征

高性能計(jì)算應(yīng)用軟件種類眾多,不可能一一羅列.現(xiàn)以美國(guó)國(guó)家能源研究科學(xué)計(jì)算中心(National Energy Research Scientific Computing Center,簡(jiǎn)稱NERSC)的計(jì)算資源使用情況作為代表來分類.NERSC 是美國(guó)能源部科學(xué)局的旗艦科學(xué)計(jì)算單位,在通過計(jì)算加速科學(xué)發(fā)現(xiàn)的領(lǐng)域處于世界領(lǐng)先地位.圖1 是2018 年NERSC 超級(jí)計(jì)算機(jī)時(shí)使用量按學(xué)科分布圖(數(shù)據(jù)來源:https://www.nersc.gov/assets/Uploads/2018NERSCAnnualReport.pdf),可見其超算應(yīng)用主要集中在地球與環(huán)境、化學(xué)、材料科學(xué)、高能物理、天文/宇宙學(xué)等學(xué)科領(lǐng)域.

Fig.1 Usage amount in subjects on NERSC supercomputer in 2018圖1 2018 年NERSC 超級(jí)計(jì)算機(jī)按學(xué)科分的使用量

本文也大致選擇以上主要計(jì)算應(yīng)用領(lǐng)域來討論其中的應(yīng)用軟件核心算法和軟件特征.這些領(lǐng)域的研究對(duì)象具有不同的計(jì)算問題研究尺度,按尺度量級(jí)從大到小排序見表1.“計(jì)算模擬”較“實(shí)驗(yàn)”和“理論”兩大科學(xué)研究手段有很多優(yōu)勢(shì):可以對(duì)不同尺度問題建模,可以對(duì)不同數(shù)據(jù)維度建模,也可以對(duì)確定性問題(如偏微分方程)或非確定性問題(如概率統(tǒng)計(jì)問題)建模.以下列舉的幾大類領(lǐng)域應(yīng)用將說明相關(guān)的建模方式.

Table 1 Scale size of different scientific subjects表1 不同科學(xué)研究對(duì)象的尺度大小

1.1 宇宙N體模擬

N 體問題又稱“多體問題”,兼具數(shù)據(jù)密集和計(jì)算密集的特點(diǎn),是物理科學(xué)高端模擬中7 個(gè)主要問題之一[17].此類問題在天體物理學(xué)、分子動(dòng)力學(xué)等領(lǐng)域處于核心的地位.天文學(xué)的一些關(guān)鍵現(xiàn)象(如星系成團(tuán)及演化、暗物質(zhì)暈等)需采用大尺度、高解析度的大規(guī)模N 體模擬結(jié)果來解釋和指導(dǎo)天文觀測(cè);宇宙模型和理論也通過計(jì)算結(jié)果得以驗(yàn)證和完善.隨著先進(jìn)快速算法和高性能計(jì)算技術(shù)發(fā)展,萬億量級(jí)規(guī)模N 體問題的高效大規(guī)模并行計(jì)算模擬已成為可能.

基于樹遍歷的算法如快速多極子方法(fast multipole method,簡(jiǎn)稱FMM),不僅具有近似線性的復(fù)雜度和較高的精度,而且適合并行.其中短程力的求解具有計(jì)算密集的特征,可有效利用高性能異構(gòu)平臺(tái)加速器件(如GPU、MIC 或眾核處理器等)的超強(qiáng)計(jì)算能力.基于快速傅里葉變換(fast Fourier transform,簡(jiǎn)稱FFT)的粒子網(wǎng)格方法(particle-mesh method,簡(jiǎn)稱PM)雖速度更快,但由于受限于精度,只適用于計(jì)算長(zhǎng)程力,而且將帶來全局通信.為兼顧計(jì)算速度和精度,重疊短程力的計(jì)算與長(zhǎng)程力的通信,混合算法成為高性能并行N 體模擬的主流[18,19].如圖2 所示,虛線內(nèi)為FMM 計(jì)算區(qū)域,其中粒子間的短程力相互作用(P2P)占FMM 計(jì)算時(shí)間的90%以上,適合在加速卡上實(shí)現(xiàn);虛線外的區(qū)域采用PM 計(jì)算長(zhǎng)程力,可在CPU端調(diào)用FFT 實(shí)現(xiàn).

Fig.2 A schematic diagram of FMM-PM algorithm圖2 FMM-PM算法示意圖

中國(guó)科學(xué)院國(guó)家天文臺(tái)和計(jì)算機(jī)網(wǎng)絡(luò)信息中心團(tuán)隊(duì)合作,在國(guó)產(chǎn)高性能異構(gòu)計(jì)算平臺(tái)上采用FMM 與PM耦合算法實(shí)現(xiàn)了大規(guī)模的N 體問題高效并行模擬.計(jì)算的并行效率達(dá)80%,浮點(diǎn)峰值效率達(dá)30%以上,為實(shí)現(xiàn)更高分辨率的超大規(guī)模天文應(yīng)用模擬奠定了基礎(chǔ).

1.2 大氣與地球系統(tǒng)模式

大氣與地球系統(tǒng)模式中,Stencil 計(jì)算是核心算法,也是物理科學(xué)高端模擬中7 個(gè)主要問題之一[17]并廣泛存在與科學(xué)與工程計(jì)算中.Stencil 計(jì)算(如圖3 所示)絕大部分是訪存密集型計(jì)算,計(jì)算訪存比低,主存帶寬受限嚴(yán)重[20].已經(jīng)有很多 Stencil 計(jì)算優(yōu)化方面的工作:如 Strumpen 等人[21]利用時(shí)間局部性優(yōu)化三級(jí)存儲(chǔ)系統(tǒng);Micikevicus 等提出了GPU 平臺(tái)上的高階Stencil 優(yōu)化方法.此外,也有一些基于Stencil 計(jì)算的應(yīng)用在國(guó)產(chǎn)機(jī)器上的移植優(yōu)化工作,如Fu 等人[22]將大氣模式CAM 和大型非線性地震模擬系統(tǒng)移植到“神威·太湖之光”超級(jí)計(jì)算機(jī)上.

與大氣模式CAM 類似,中國(guó)科學(xué)院地球系統(tǒng)模式中的有限差分Stencil 計(jì)算的運(yùn)行時(shí)間占比超過90%,但是其浮點(diǎn)計(jì)算效率僅有不到5%.而且這些Stencil 計(jì)算包含更多更復(fù)雜的數(shù)據(jù)結(jié)構(gòu),Stencil 的種類和差分格式、算法求解過程和特點(diǎn)也各不相同.其主要特點(diǎn)有:通信和同步多、復(fù)雜計(jì)算多、變量多而計(jì)算網(wǎng)格點(diǎn)少、計(jì)算訪存比和主存帶寬受限問題更突出.這對(duì)當(dāng)前的國(guó)產(chǎn)并行計(jì)算平臺(tái)的計(jì)算和訪存能力提出了新需求.而在高性能計(jì)算機(jī)普遍采用異構(gòu)架構(gòu)的當(dāng)下,性能調(diào)優(yōu)顯得尤其重要.代碼的移植優(yōu)化工作是一項(xiàng)具有挑戰(zhàn)性的工作.主要的優(yōu)化策略有:通過通信優(yōu)化手段減少通信占比,提高程序的可擴(kuò)展性;通過計(jì)算子區(qū)域規(guī)模劃分,調(diào)節(jié)計(jì)算資源和高速寄存器及緩存使用率,達(dá)到較好的均衡;通過改進(jìn)算法或算法實(shí)現(xiàn),減少數(shù)據(jù)依賴和分支計(jì)算,提高數(shù)據(jù)訪問連續(xù)性和cache 命中率;提高程序的向量化率等.

Fig.3 Discrete structure in Stencil computation圖3 Stencil 計(jì)算的離散結(jié)構(gòu)

1.3 相場(chǎng)模擬微結(jié)構(gòu)演化

合金材料的微觀組織對(duì)材料的諸多性能起到?jīng)Q定性的影響,因此深入研究其微觀組織的演化機(jī)理具有重要的意義.相場(chǎng)方法是模擬合金材料微觀組織演化的主要手段,主要通過求解關(guān)于相場(chǎng)、成分場(chǎng)等變量的一組偏微分方程來進(jìn)行微組織演化過程的模擬.方程具有強(qiáng)剛性和強(qiáng)非線性的特點(diǎn),大體系長(zhǎng)時(shí)間的模擬是非常具有挑戰(zhàn)性的課題.傳統(tǒng)的有限差分顯式算法雖可在當(dāng)代高性能計(jì)算平臺(tái)上有較高效的實(shí)現(xiàn)[23],但由于時(shí)間步長(zhǎng)受到算法穩(wěn)定性條件的限制而無法滿足長(zhǎng)時(shí)間模擬的需求.而隱格式的設(shè)計(jì)則面臨著相場(chǎng)模型方程強(qiáng)非線性以及相應(yīng)的大規(guī)模線性和非線性代數(shù)方程組并行求解等問題.中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心和金屬研究所團(tuán)隊(duì)合作,針對(duì)國(guó)產(chǎn)異構(gòu)高性能計(jì)算平臺(tái)設(shè)計(jì)了可擴(kuò)展指數(shù)時(shí)間差分算法,可高效地進(jìn)行相場(chǎng)模型的大體系長(zhǎng)時(shí)間模擬.2016 年在“神威·太湖之光”超級(jí)計(jì)算機(jī)上使用全機(jī)約1 000 萬計(jì)算核心實(shí)現(xiàn)了持續(xù)性能為50 PFLOPS(每秒5 億億次浮點(diǎn)運(yùn)算)的超大規(guī)模相場(chǎng)模擬[24].基于“神威”超級(jí)計(jì)算機(jī)的“申威”處理器的主從核架構(gòu),數(shù)據(jù)通信使用主核DMA 方式,而七點(diǎn)Stencil 計(jì)算則在從核上計(jì)算.最終應(yīng)用測(cè)試(如圖4 所示)的并行效率達(dá)到70%,浮點(diǎn)性能達(dá)到峰值的41%,遠(yuǎn)超同類的計(jì)算軟件.

Fig.4 Configuration coarsening evolution of two-phase system governed by the Cahn-Hilliard equation圖4 兩相體系由Cahn-Hilliard 方程控制在粗化過程中的構(gòu)型演化

1.4 分子動(dòng)力學(xué)模擬

分子動(dòng)力學(xué)模擬(molecular dynamics,簡(jiǎn)稱MD)是在原子水平上利用經(jīng)典牛頓力學(xué)方程研究分子狀態(tài)的計(jì)算機(jī)模擬方法[25].隨著計(jì)算機(jī)計(jì)算能力的提高,可實(shí)現(xiàn)全原子模型的動(dòng)力學(xué)模擬,在原子結(jié)構(gòu)層次揭示某些現(xiàn)象隨時(shí)間的發(fā)展過程,已廣泛應(yīng)用于計(jì)算化學(xué)、生物化學(xué)、納米材料和藥物分子設(shè)計(jì)等領(lǐng)域.

分子動(dòng)力學(xué)模擬的計(jì)算核心在于力場(chǎng)的計(jì)算,具體可細(xì)分為長(zhǎng)程作用力計(jì)算、短程力計(jì)算與近鄰列表的計(jì)算、成鍵作用力計(jì)算等多種形式;除了核心的力場(chǎng)計(jì)算,還存在著一些物理系統(tǒng)約束需要實(shí)現(xiàn),比如分子結(jié)構(gòu)約束方法SHAKE 和LINCS、周期性邊界條件和統(tǒng)計(jì)系綜理論約束等.加速分子動(dòng)力學(xué)模擬,在方法上可采用增強(qiáng)抽樣算法,在實(shí)現(xiàn)上則可借助大規(guī)模并行計(jì)算.在并行計(jì)算時(shí)任務(wù)劃分常采用體系空間分解法,可采用MPI+X的多級(jí)并行計(jì)算設(shè)計(jì),其中X 可為多核CPU 計(jì)算或GPU 異構(gòu)計(jì)算等形式.比如美國(guó)能源部主導(dǎo)開發(fā)的MD 軟件LAMMPS(large-scale atomic/molecular massively parallel simulator),可利用多核CPU、GPU 和Intel Xeon Phi等多種方式加速其計(jì)算性能.該軟件提供了GPU package 和KOKKOS package 可利用GPU 部件加速計(jì)算,使用MPI 計(jì)算任務(wù)下的GPU 任務(wù)offload 工作方式[26].這種方式將計(jì)算密集性任務(wù)交給GPU 來計(jì)算完成,并獲得返回結(jié)果,充分利用了異構(gòu)計(jì)算環(huán)境來提高并行效率.結(jié)合區(qū)域空間分解和MPI 并行技術(shù),可進(jìn)行數(shù)萬核以上的超大規(guī)模并行計(jì)算.GPU 作為異構(gòu)計(jì)算設(shè)備用來加速分子動(dòng)力學(xué)模擬的一種實(shí)現(xiàn)流程如圖5 所示,不同軟件在實(shí)現(xiàn)上可能存在差異,但基本用于加速分子作用力計(jì)算以及相關(guān)的近鄰列表構(gòu)建等計(jì)算密集型任務(wù).

Fig.5 A flowchart for an implementation of GPU accelerated molecular dynamics simulation[27]圖5 一種實(shí)現(xiàn)GPU 加速分子動(dòng)力學(xué)模擬的流程圖[27]

1.5 量子化學(xué)

量子化學(xué)源自理論化學(xué),應(yīng)用量子力學(xué)的基本原理和方法來研究化學(xué)問題.經(jīng)過多年的發(fā)展,已經(jīng)有多種量子化學(xué)計(jì)算軟件被開發(fā)出來.科學(xué)家們用這些工具計(jì)算和求解分子結(jié)構(gòu)問題,不斷地?cái)U(kuò)大其可解問題的深度和廣度.

自洽場(chǎng)(self-consistent field,簡(jiǎn)稱SCF)迭代計(jì)算是量子化學(xué)計(jì)算軟件的核心方法.它主要包括計(jì)算雙電子積分、構(gòu)造Fock 矩陣和求解廣義特征值方程,兼具計(jì)算密集型和訪存密集型的特點(diǎn).計(jì)算雙電子積分一般采用解析的方法,是計(jì)算密集型的,計(jì)算量與體系大小的4 次方成正比.在引入了prescreening 方法后,可以有效降低積分計(jì)算的scaling,同時(shí)利用分子劃片的方法,可以得到近似線性的scaling 結(jié)果.構(gòu)造Fock 矩陣是訪存密集型的,有大量的隨機(jī)不連續(xù)訪存.如何協(xié)調(diào)處理計(jì)算積分和構(gòu)造矩陣是異構(gòu)優(yōu)化自洽場(chǎng)計(jì)算的難點(diǎn).2008年,Yasuda[28]首先提出利用GPU 計(jì)算雙電子積分.同年利用積分公式自動(dòng)展開的方法將積分的計(jì)算轉(zhuǎn)化為適合GPU 運(yùn)行的方式,取得了較好的結(jié)果.后來,Ufimtsev 等人[29,30]提出了預(yù)先排序的方法,實(shí)現(xiàn)了在GPU 上構(gòu)造Fock 矩陣(如圖6 所示).盡管算法相對(duì)于傳統(tǒng)的CPU算法在計(jì)算量上有一定的損失,但仍然取得了近百倍的加速效果(相對(duì)于1 個(gè)CPU 核).

Fig.6 A flowchart of GPU accelerated SCF calculation圖6 GPU 加速SCF 計(jì)算流程圖

后自洽場(chǎng)(post-SCF)計(jì)算在提高量子化學(xué)計(jì)算結(jié)果精度和研究前沿科學(xué)問題上有著非常重要的作用.作為多參考方法起點(diǎn)的多組態(tài)(multi-configurational)自洽場(chǎng)的計(jì)算熱點(diǎn)非常類似Hartree-Fock 的計(jì)算,區(qū)別之處是在多組態(tài)自洽場(chǎng)使用的積分常為分子軌道下的積分,自洽迭代過程中積分的多個(gè)指標(biāo)需要進(jìn)行持續(xù)的變換操作.此部分計(jì)算同樣兼具訪存密集型和計(jì)算密集型的特點(diǎn).基于多行列式的微擾計(jì)算則是成功地實(shí)現(xiàn)了GPU 加速,極具代表性的有NWchem 軟件中的CCSD(T)方法.該法中的CCSD 計(jì)算復(fù)雜,由CPU 處理;后續(xù)基于CCSD波函數(shù)的微擾計(jì)算(即T 部分)則交由計(jì)算能力較強(qiáng)的GPU 來實(shí)現(xiàn),需針對(duì)偶合簇波函數(shù)中的每一個(gè)行列式表示的系數(shù)進(jìn)行微擾校正,計(jì)算熱點(diǎn)為密度矩陣同電子積分的乘積運(yùn)算.

1.6 格點(diǎn)量子色動(dòng)力學(xué)

量子色動(dòng)力學(xué)(即QCD)是一種強(qiáng)相互作用的規(guī)范理論,它和電弱統(tǒng)一理論(EW)通稱為標(biāo)準(zhǔn)模型,是高能物理研究的基本理論體系.QCD 描述了組成強(qiáng)作用粒子(強(qiáng)子)的夸克和與色量子數(shù)相聯(lián)系的規(guī)范場(chǎng)的相互作用,可統(tǒng)一地描述強(qiáng)子的結(jié)構(gòu)和它們之間的強(qiáng)相互作用,被認(rèn)為是有希望的強(qiáng)作用基本理論.格點(diǎn)量子色動(dòng)力學(xué)(格點(diǎn)QCD)是從第一性原理出發(fā)的研究QCD 的非微擾方法[31],對(duì)于高能物理中的重大物理問題,如夸克禁閉、宇宙早期演化、標(biāo)準(zhǔn)模型精確檢驗(yàn)及超出標(biāo)準(zhǔn)模型新物理需尋找等研究有不可替代的作用.格點(diǎn)QCD 方法的計(jì)算模擬過程對(duì)高性能計(jì)算資源需求巨大,其計(jì)算過程體現(xiàn)出計(jì)算密集型特征,而其數(shù)據(jù)分析則體現(xiàn)出數(shù)據(jù)密集型特征.計(jì)算精度和數(shù)據(jù)分析效率直接反映了格點(diǎn)QCD 計(jì)算方法的應(yīng)用水平.大規(guī)模計(jì)算應(yīng)用中通常采用空間區(qū)域分解的并行計(jì)算方式(如圖7 所示),格點(diǎn)子區(qū)域上的計(jì)算任務(wù)分配到不同的處理器上計(jì)算,各處理器核負(fù)責(zé)格點(diǎn)邊界上的數(shù)據(jù)通信,總體上使得數(shù)值模擬計(jì)算具有高可并行性和高可擴(kuò)展性等特點(diǎn),是高性能科學(xué)計(jì)算的典型應(yīng)用領(lǐng)域[32],曾獲得1988、1998、2006 年的Gordon Bell 獎(jiǎng),2018 年也入圍該獎(jiǎng)項(xiàng)的決賽[33].

Fig.7 A schematic diagram of domain decomposition in parallel computing lattice QCD圖7 格點(diǎn)QCD 區(qū)域劃分并行計(jì)算示意圖

在格點(diǎn)QCD 軟件設(shè)計(jì)上,抽象分層軟件架構(gòu)的設(shè)計(jì)方法得到了重要應(yīng)用,比如美國(guó)的USQCD 軟件框架分為4 層:最上層應(yīng)用層包括Chroma,之下3 層為支撐庫(kù)層,分別進(jìn)行消息通信、數(shù)據(jù)結(jié)構(gòu)定義和基本算法實(shí)現(xiàn).再如英國(guó)的Grid 軟件[6]也是抽象了數(shù)據(jù)并行計(jì)算模型,針對(duì)格點(diǎn)QCD 的熱點(diǎn)計(jì)算抽象了統(tǒng)一的數(shù)據(jù)并行模型,支持多種平臺(tái)的SIMD 矢量化計(jì)算,如AVX2 和AVX512 等形式.這種數(shù)據(jù)并行模型可以有效支撐上層應(yīng)用的開發(fā),保證較高的高性能計(jì)算效率.軟件框架的這種分層抽象的設(shè)計(jì)方法極大地提高了代碼的開發(fā)效率和運(yùn)行效率,大大地促進(jìn)了整個(gè)格點(diǎn)QCD 研究領(lǐng)域的高速發(fā)展.針對(duì)異構(gòu)計(jì)算環(huán)境,美國(guó)為主的團(tuán)隊(duì)開發(fā)了針對(duì)GPU加速的格點(diǎn)QCD 軟件包QUDA,軟件包主要采用C++模版類支持不同數(shù)據(jù)類型的加速計(jì)算,已在美國(guó)Summit超級(jí)計(jì)算機(jī)上實(shí)現(xiàn)大規(guī)模GPU 加速的格點(diǎn)QCD 計(jì)算,取得了顯著的加速效果.基于CUDA 設(shè)計(jì)的QUDA 軟件如要移植到曙光國(guó)產(chǎn)超算平臺(tái),需要經(jīng)過轉(zhuǎn)碼為ROCm 支持的HIP 代碼形式.

2 共性問題的提煉及解決對(duì)策

2.1 共性算法

2004 年著名高性能計(jì)算專家Colella[17]從計(jì)算物理問題中提出了科學(xué)計(jì)算的7 種算法模型,即結(jié)構(gòu)性網(wǎng)格、非結(jié)構(gòu)性網(wǎng)格、快速傅里葉變換(FFT)、稠密線性代數(shù)(dense linear algebra)、稀疏線性代數(shù)(sparse linear algebra)、粒子動(dòng)力學(xué)(particles)、蒙特卡洛(Monte Carlo).這7 種算法模型還可以進(jìn)一步歸類為求解數(shù)值微分法方程類和求解數(shù)值線性代數(shù)類,實(shí)現(xiàn)形式上有顯性計(jì)算方式和隱性計(jì)算方式.這些核心的共性計(jì)算模型對(duì)高性能計(jì)算機(jī)提出了不同方面和層次的需求,也基本覆蓋了常見高性能計(jì)算應(yīng)用.如N 粒子相互作用中等設(shè)計(jì)FFT計(jì)算長(zhǎng)程相互作用,粒子動(dòng)力學(xué)是一個(gè)微分方程的求解,格點(diǎn)QCD 計(jì)算也會(huì)涉及大量的線性系統(tǒng)求解,通常為千萬到億級(jí)矩陣維度的巨大計(jì)算量,同時(shí)存在海量3x3 低維向量矩陣乘法.提煉以上共性算法特征,有利于集中研發(fā)力量開發(fā)可共享使用的軟件算法包;也有利于突破計(jì)算關(guān)鍵技術(shù)問題,提高高性能計(jì)算軟件研發(fā)的效率.

國(guó)產(chǎn)軟件編程框架JASMIN[14]和JAUMIN[15]核心正是從共性算法提煉出發(fā),對(duì)大量物理模型計(jì)算問題進(jìn)行抽象和提煉共性特征,形成編程框架的核心.編程框架底層以網(wǎng)格數(shù)據(jù)模型為基礎(chǔ),提供數(shù)值算法構(gòu)件,中間層實(shí)現(xiàn)并行計(jì)算模式和負(fù)載平衡建模,頂層為數(shù)值代數(shù)解法器庫(kù)或者獨(dú)立于編程框架的第三方解法器庫(kù)接口.

2.2 浮點(diǎn)計(jì)算精度

不論何種科學(xué)計(jì)算應(yīng)用,浮點(diǎn)計(jì)算精度的問題均扮演著重要的角色.從應(yīng)用問題的角度來看,數(shù)值結(jié)果精度往往越高越好,但高精度浮點(diǎn)計(jì)算較低精度浮點(diǎn)計(jì)算需要更多的時(shí)間.一般64 位的雙精度計(jì)算所需時(shí)間是32位的單精度浮點(diǎn)計(jì)算時(shí)間的兩倍,因此常常遇到在高精度和高計(jì)算性能之間的取舍問題.高低精度混合計(jì)算是一般的處理方式,但大型科學(xué)計(jì)算模擬中浮點(diǎn)運(yùn)算過程非常復(fù)雜,很難從理論上推導(dǎo)出合理的精度位,不得不憑經(jīng)驗(yàn)判斷,但這不可靠更不科學(xué).一個(gè)相對(duì)科學(xué)的方式是先評(píng)估科學(xué)計(jì)算應(yīng)用算法中最少需要多少位數(shù)值有效位數(shù),而這個(gè)有效位數(shù)能滿足對(duì)結(jié)果需求精度.這可以借助軟件工具Verificarlo[34]來推算出計(jì)算算法中保持結(jié)果所需的最低有效位,然后選擇相匹配的浮點(diǎn)數(shù)表示精度.另外,在某些體系結(jié)構(gòu)支持多個(gè)低精度浮點(diǎn)拼湊為高精度浮點(diǎn)的運(yùn)算,但這增加了算法實(shí)現(xiàn)的復(fù)雜度.在通用處理器中一般只提供浮點(diǎn)精度計(jì)算,對(duì)應(yīng)于C 語(yǔ)言中的double 和float 精度;而專用FPGA 加速模塊卻可以自定義浮點(diǎn)計(jì)算精度;最近的Nvidia GPU 也支持半精度的運(yùn)算.科學(xué)計(jì)算中浮點(diǎn)運(yùn)算追求的目標(biāo)之一是所謂的“高準(zhǔn)確率低精度(high-accuracy low-precision,簡(jiǎn)稱HALP)”[35].

2.3 海量數(shù)據(jù)存儲(chǔ)及可視化

高性能計(jì)算應(yīng)用在大規(guī)模模擬中,往往會(huì)遇到大量數(shù)據(jù)存儲(chǔ)的問題.比如分子模擬常規(guī)100 萬原子體系中,每個(gè)原子單精度浮點(diǎn)數(shù)僅存儲(chǔ)三維的速度和位移數(shù)據(jù),大約需要24MB 存儲(chǔ)空間,如果需要1 萬幀軌跡數(shù)據(jù)則需要240GB 文件,如果再需要不同溫度和壓強(qiáng)的組合條件下的軌跡數(shù)據(jù),對(duì)存儲(chǔ)空間需求劇增.另外,數(shù)據(jù)格式在異構(gòu)系統(tǒng)上不兼容,直接數(shù)據(jù)傳輸存儲(chǔ)會(huì)產(chǎn)生問題.當(dāng)前可以借助專門的HDF(hierarchical data format)文件格式接口來避免,它可存儲(chǔ)不同類型的圖像和數(shù)碼數(shù)據(jù)的文件格式,并可在不同類型的機(jī)器上傳輸,同時(shí)還有統(tǒng)一處理這種文件格式的函數(shù)庫(kù).

數(shù)值模擬數(shù)據(jù)的可視化特點(diǎn)見表2,數(shù)據(jù)源可以按網(wǎng)格特征進(jìn)行劃分,數(shù)據(jù)的可視化也具有一般的流程.大規(guī)模數(shù)值模擬數(shù)據(jù)對(duì)其可視化也帶來了困難,主要體現(xiàn)在空間表示和空間存儲(chǔ)方面.比如粒子數(shù)超過了32 位整數(shù)可表示的范圍可能導(dǎo)致軟件崩潰,渲染時(shí)GPU 顯存一次性無法容納一個(gè)數(shù)據(jù)文件時(shí),導(dǎo)致數(shù)據(jù)分批導(dǎo)入顯存,增加了算法的復(fù)雜度也往往降低了計(jì)算性能.

Table 2 Visualization features in numerical simulation表2 數(shù)值模擬數(shù)據(jù)可視化特點(diǎn)

2.4 軟件框架設(shè)計(jì)

軟件框架決定了軟件的生命力,其最主要特征是軟件功能的適應(yīng)性,具體包括理論建模支撐、計(jì)算模塊、邏輯流程控制和數(shù)據(jù)分析等模塊.隨著計(jì)算機(jī)硬件特性和軟件編程語(yǔ)言的發(fā)展,許多科學(xué)計(jì)算軟件已不能很好地適應(yīng)應(yīng)用需求,大部分問題是軟件的原框架受到?jīng)_擊導(dǎo)致的.如原有程序采用Fortran 語(yǔ)言編寫,后來為了適應(yīng)新的運(yùn)行環(huán)境改用C++實(shí)現(xiàn)大規(guī)模設(shè)計(jì),又需要新的C++11 改進(jìn)編程效率,導(dǎo)致軟件框架不穩(wěn)定.其次早期采用OpenMP 為主的多CPU 核計(jì)算模型,需要轉(zhuǎn)為異構(gòu)OpenCL 計(jì)算模式;另外由于新的加速設(shè)備,原來的CUDA 程序可能需要改造為OpenACC 并行編程方式.以上問題在高性能計(jì)算軟件的設(shè)計(jì)中普遍存在,內(nèi)在原因?yàn)檐浖⑿杏?jì)算模式較多,硬件體系結(jié)構(gòu)差異大.

對(duì)于高性能計(jì)算軟件來說,優(yōu)秀的軟件框架可提高編程開發(fā)效率,適應(yīng)新的軟件需求和用戶服務(wù)推廣.如何以不變的設(shè)計(jì)原則應(yīng)對(duì)多變的軟件需求是一個(gè)關(guān)鍵問題.對(duì)于軟件框架來說,分層式設(shè)計(jì)常用方式如通用并行軟件JASMIN 和格點(diǎn)QCD 應(yīng)用軟件Chroma.但從異構(gòu)計(jì)算這方面考慮,面向異構(gòu)計(jì)算環(huán)境有機(jī)統(tǒng)一計(jì)算模擬的各方面是根本性解決問題,如統(tǒng)一處理器核計(jì)算、主存和加速內(nèi)存管理和網(wǎng)絡(luò)通信等三大模塊,屏蔽異構(gòu)計(jì)算帶來的差異性可以使上層軟件設(shè)計(jì)保持不變性;當(dāng)前Khronos 組織發(fā)表的SYCL 規(guī)范[36]和Intel oneAPI[37]正在建立面向異構(gòu)計(jì)算環(huán)境的統(tǒng)一編程接口,SYCL 底層基于OpenCL 跨平臺(tái)異構(gòu)API 規(guī)范,基于SYCL 的DPC++作為Intel oneAPI 的主體,已經(jīng)實(shí)現(xiàn)了Intel GPU 和Nvidia GPU 的異構(gòu)加速計(jì)算;針對(duì)AMD GPU 常用選擇ROCm 框架下的HIP 編程方式,而HIP 實(shí)際上能夠兼容CUDA 在Nvidia GPU 上編譯運(yùn)行,另外,值得注意的是,基于HIP 也可以實(shí)現(xiàn)SYCL 統(tǒng)一規(guī)范[38].

3 總結(jié)

本文基于高性能計(jì)算應(yīng)用軟件的現(xiàn)狀趨勢(shì),對(duì)宇宙N 體模擬、大氣與地球系統(tǒng)模式、計(jì)算材料相場(chǎng)動(dòng)力學(xué)、分子動(dòng)力學(xué)、量子計(jì)算化學(xué)和格點(diǎn)量子色力學(xué)等應(yīng)用開展了并行計(jì)算特征分析,提煉了應(yīng)用的典型算法和軟件的共性問題,并開展了解決的對(duì)策討論.

在國(guó)產(chǎn)異構(gòu)計(jì)算高性能計(jì)算環(huán)境下,由于國(guó)產(chǎn)架構(gòu)不同于主流的體系架構(gòu),其算法實(shí)現(xiàn)和軟件設(shè)計(jì)有很大差異.另外我國(guó)自主的硬件體系架構(gòu)設(shè)計(jì)起步較晚,所構(gòu)建的軟件生態(tài)環(huán)境包括系統(tǒng)工具和程序調(diào)試工具等等還遠(yuǎn)不成熟,給應(yīng)用程序開發(fā)團(tuán)隊(duì)開發(fā)程序增加了現(xiàn)實(shí)難度.為了有效開展基于國(guó)產(chǎn)高性能異構(gòu)計(jì)算算法和軟件設(shè)計(jì),現(xiàn)給出兩點(diǎn)發(fā)展建議:

(1) 盡快建立我國(guó)高性能計(jì)算的統(tǒng)一軟硬件標(biāo)準(zhǔn).可從指令集、API 接口、軟件架構(gòu)、網(wǎng)格環(huán)境對(duì)接等方面建立統(tǒng)一規(guī)范.

(2) 兼容已有國(guó)際主流標(biāo)準(zhǔn)并積極參與國(guó)際規(guī)范標(biāo)準(zhǔn)的制定,提高國(guó)產(chǎn)超級(jí)計(jì)算平臺(tái)與開源軟件的兼容性.