神威太湖之光上分子動力學(xué)模擬的性能優(yōu)化?

田 卓,陳一峯

(北京大學(xué) 信息科學(xué)與技術(shù)學(xué)院,北京 100 871)

分子動力學(xué)模擬是考察系統(tǒng)隨時間演化的行為[1,2].在一定時間內(nèi),靠計算機模擬分子和原子體系運動,廣泛應(yīng)用于物理、化學(xué)、生物等領(lǐng)域.用分子動力學(xué)模擬微觀生物現(xiàn)象如蛋白質(zhì)折疊[3,4]所需時間尺度為毫秒,為捕獲原子震蕩現(xiàn)象,每個迭代步設(shè)為1fs~2.5fs,模擬總時長若為毫秒,需要1012個迭代步.相鄰兩個迭代步間具有相關(guān)性,需要一次通信,而頻繁通信成為此類問題計算性能的主要瓶頸.

時間演化類問題從時間軸來模擬體系演化.時間被離散化為一定步長為間隔的時間點,下一個時間點的狀態(tài)依賴上一個時間點的狀態(tài),兩個狀態(tài)之間需要一次通信以交換數(shù)據(jù).典型的應(yīng)用是分子動力學(xué)模擬[5,6],分子沿著這些時間點運動,每個時間點對應(yīng)分子的一個狀態(tài).分子在t時刻的狀態(tài)用向量X(t)表示,記錄t時刻所有粒子的位置,速度等狀態(tài)信息.t+1 時刻的狀態(tài)X(t+1)由t時刻狀態(tài)X(t)演算而來,即X(t+1)=f(X(t)).因而,不同時刻粒子狀態(tài)間的關(guān)系是相互依賴的,需要通信.而神威上的通信延遲較長[7],包括網(wǎng)卡上的通信延遲以及較長的數(shù)據(jù)通路所帶來的通信延遲都太長了.此類計算模式的時間演化類程序在高延遲的神威集群上如何優(yōu)化,規(guī)避通信延遲,提高性能,是本文的宗旨.

神威太湖之光超級計算機以理論峰值125.4PFLOPS/s 位居世界第一,其SW26010 處理器是針對神威定制的多核處理器,每個處理器由4 個核組組成,每個核組包含1 個主核和64 個從核.神威超級計算機上共安裝了40 96 0 個SW26010 處理器,共10 64 9 600 個核[8].在神威上,大部分應(yīng)用是以一個CPU 作為單位編程[9],只用4個核組中的一個核組進行計算和通信,如圖1所示.CPU 上的4 個核組共享一個網(wǎng)絡(luò)端口,若問題規(guī)模不變,這就迫使我們充分利用所有核組的計算能力,以提高性能,如圖2,4 個主核同時計算和通信.由圖可見,兩種模式的區(qū)別是:(1)進程數(shù)不同;(2)消息個數(shù)不同.設(shè)問題規(guī)模為N,進程個數(shù)為P,模式1,進程上任務(wù)量為N/P,消息個數(shù)m.模式2,每進程上任務(wù)量N/4P,消息個數(shù)4×m.消息個數(shù)增加了4 倍,因而帶來了通信延遲.如何在以核組為單位的編程模式中減少消息個數(shù),優(yōu)化通信延遲.

如圖3所示,每個CPU 只選一個進程用于與外界通信[10],片上所有核間的數(shù)據(jù)同步是通過共享片外存儲來實現(xiàn),減少了用于通信的進程數(shù),即實現(xiàn)了消息個數(shù)的減少,優(yōu)化了通信延遲.這是神威上較快的共享數(shù)據(jù)的方式[7],但由于多核訪存需要經(jīng)過訪存隊列,這種方式仍然較慢,帶來訪存延遲.

Fig.1 Node master process mode圖1 節(jié)點主進程方式

Fig.2 Core groups programming圖2 核組為單位編程

Fig.3 Reducing the number of processes圖3 減少進程數(shù)方式

針對諸如分子動力學(xué)模擬等延遲敏感的時間演化類應(yīng)用,本文的主要貢獻是,給出了完整的主從核異構(gòu)的國產(chǎn)神威處理器上的一系列優(yōu)化技術(shù),為類似的通信受限類程序提供了參考.

(1)異構(gòu)體系結(jié)構(gòu)下核間同步的優(yōu)化,采用共享內(nèi)存等待與從核同步相結(jié)合的方式,優(yōu)化訪存延遲,進而優(yōu)化通信延遲;

(2)改變計算模式,采取有利于減少核間數(shù)據(jù)關(guān)聯(lián)和依賴關(guān)系的計算模式,減少相鄰迭代步間的通信次數(shù)和數(shù)據(jù)等待,優(yōu)化通信延遲;

(3)高效訪存及規(guī)則化以提升性能.數(shù)據(jù)傳輸與計算重疊,以優(yōu)化從核訪主存效率.規(guī)則化數(shù)據(jù)結(jié)構(gòu)以提升性能,使得向量長度對齊、連續(xù)、大小規(guī)整.

本文第1 節(jié)對神威太湖之光的結(jié)構(gòu)及實例的計算模式進行描述.第2 節(jié)給出神威上的幾種優(yōu)化技術(shù),并進行詳細闡述.第3 節(jié)對實現(xiàn)結(jié)果進行描述.第4 節(jié)對已有研究工作進行分析和總結(jié).第5 節(jié)總結(jié)全文.

1 神威結(jié)構(gòu)及計算實例

1.1 神威結(jié)構(gòu)

神威太湖之光超級計算機是基于主從核異構(gòu)的體系結(jié)構(gòu)[7],如圖4所示.處理器芯片SW26010 由4 個核組(core gro ups,簡稱CGs)構(gòu)成,每個核組包含一個主核MPE(運算控制核心,management pro cessing e lement,簡稱MPE)和一個CPE 集群(computing p rocessing elements,簡稱CPEs).CPE 集群由64 個從核按照8×8 的網(wǎng)狀格式組成,64 個從核陣列之間可采用低延遲的寄存器方式通信.4 個核組間通過片上網(wǎng)絡(luò)(network on chip,簡稱Noc)連接,每個核組有自己的內(nèi)存空間,通過內(nèi)存控制器(memory controller,簡稱MC)連接到MPE 和CPEs 上.協(xié)議處理單元(protocol processing unit,簡稱PPU)負責(zé)處理來自CPEs 和MPE 不同類型的請求.所有核組通過系統(tǒng)接口(system interface,簡稱SI)與外界相連.

Fig.4 Architecture of the SW26010 processor圖4 神威體系結(jié)構(gòu)

每個從核有64KB 的SPM(scratch pad memory).為了提高內(nèi)存帶寬的使用率,從核可通過DMA 方式批量訪問主存(main memory),DMA 通道能夠高效地在主存和SPM 間傳遞數(shù)據(jù).國產(chǎn)眾核服務(wù)器的高性能集群上,每個節(jié)點大小為4,節(jié)點內(nèi)相對進程號為0 的進程是節(jié)點主進程.節(jié)點內(nèi)4 個進程共享同一個網(wǎng)絡(luò)端口,由系統(tǒng)接口SI 連接.在優(yōu)化算法中,我們選用0 號主進程用于通信,其余3 個進程與0 號進程間通過片上共享數(shù)據(jù)交互信息.

申威異構(gòu)的體系架構(gòu)下,若問題本身的計算具有高頻迭代性質(zhì),同時不同處理單元上的所分配的計算任務(wù)之間具有相互依賴性,由此產(chǎn)生的相互依賴的計算需要在多個異構(gòu)單元下進行高頻通信以及同步訪存.本文提出的優(yōu)化策略能夠為此類延遲敏感類應(yīng)用的性能優(yōu)化提供一定的參考.本文的算例是基于時間演化類程序的典型代表分子動力學(xué)模擬為例,此外,本文提出的優(yōu)化策略同樣普適于類似的高頻迭代類應(yīng)用如stencil 計算,它是很多科學(xué)計算類應(yīng)用最為重要和耗時最多的計算核心.此類應(yīng)用的特征是問題本身的求解需要大量高頻迭代,迭代算法耗時較多.此類應(yīng)用在異構(gòu)體系結(jié)構(gòu)下的性能瓶頸是系統(tǒng)較長的延遲,性能優(yōu)化的重點在于提高迭代頻率.

1.2 計算實例

分子動力學(xué)模擬是指用計算機模擬大量粒子系統(tǒng)中每個粒子的運動過程.系統(tǒng)中每個粒子在其他粒子所形成的勢場下運動,運動方程遵循牛頓定律[11].牛頓方程如下:

其中:N表示系統(tǒng)中的原子數(shù);Mi,ri,Fi分別表示原子i的質(zhì)量、位移以及受力;Φ表示勢函數(shù),勢函數(shù)需能夠描述系統(tǒng)在時間演化的各個階段,即原子各種狀態(tài)下的相互作用.復(fù)雜系統(tǒng)的描述多采用多體勢,本文研究的硅材料,原子間的作用勢取決于原子間的距離以及原子間的成鍵方向[12,13].硅原子間的多體作用模型采用典型的共價鍵勢Tersoff 勢描述.Tersoff 勢函數(shù)的數(shù)學(xué)表示為

fR為排斥勢,fA為吸引勢,fC為截斷函數(shù),bij為鍵級項:

本文以分子動力學(xué)模擬中多體作用模型實現(xiàn)的單晶硅體系原子模擬為計算實例,說明這種計算模式的應(yīng)用在神威集群上運行時所遇到的性能瓶頸,并給出幾種優(yōu)化技術(shù).硅原子之間的多體作用模型采用Tersoff 勢函數(shù)[14],該模型的特點是原子之間的相互作用依賴于原子所處的局域環(huán)境.每個迭代步將原子t時刻的狀態(tài)推演到t+1 時刻,更新原子狀態(tài)時,需獲得該原子的所有鄰居原子的狀態(tài)信息,即每個迭代步需要一次同步.原子之間的依賴關(guān)系具有局部相關(guān)性,因而導(dǎo)致一定的異步性.此類時間演化類問題,硅原子的狀態(tài)是由時間來推進演化的,因此,不同迭代步之間是相互依賴的,需要一次通信.

考慮常溫下,硅晶體保持規(guī)則的金剛石晶胞結(jié)構(gòu),每個原子周圍有4 個鄰近原子.如圖5所示,硅原子a有4個鄰居原子b,c,d,e,每個原子由不同線程計算.t時刻,每個線程更新了5 個粒子的位置信息,進入下一個迭代步.a原子的受力受到它的所有鄰居原子的影響,需要得到4 個鄰居粒子的位置信息,來計算不同鄰居對a的分子間作用力.因此,需要一次通信來獲得所有鄰居在t時刻的位置信息.

Fig.5 Schematic of interatomic interaction圖5 硅原子間相互作用示意圖

有效的原子模擬現(xiàn)象,需要1012個迭代步,相鄰迭代步間需要一次通信,那么在實驗室環(huán)境下,觀察到有效的原子模擬現(xiàn)象是不現(xiàn)實的.為了提高硅原子模擬的迭代頻率,有效的做法是減少迭代步間通信的延遲.由神威體系高性能集群的特點可知,其上的通信延遲較長.那么針對諸如分子動力學(xué)模擬的時間演化類程序,如何在神威上減少通信延遲、提高迭代頻率,本文給出了一系列優(yōu)化技術(shù).

2 神威上的優(yōu)化

CPU 的從核線程負責(zé)計算,兩個從核線程位于CPU0 和CPU1 上,如圖6.CPU0 上的從核線程產(chǎn)生的數(shù)據(jù),是CPU1 上的從核線程下一個迭代步的輸入,輸出又作為CPU0 上從核線程下個迭代步的輸入.這樣產(chǎn)生了數(shù)據(jù)依賴關(guān)系.

Fig.6 Data communication loop between different CPUs on Sunway圖6 神威上不同CPU 間的數(shù)據(jù)通信環(huán)路

若不做任何計算,數(shù)據(jù)在兩個CPU 之間傳遞一個來回也需要花費時間.CPU 計算能力即使再快,這個通信環(huán)路花費的時間長的話,也沒有任何意義.不考慮任何計算時間,系統(tǒng)的速度極限即每秒能做多少次迭代,與這個環(huán)路每秒能做多少次通信是一樣的.經(jīng)測試,神威上的速度極限是每秒能做3 00 0 次通信,但這個速度實際上太慢了.環(huán)路上的通信延遲,包括發(fā)生在網(wǎng)卡上和每臺服務(wù)器上的延遲都太長了.對通信受限的時間演化類程序而言,通信延遲將極大地制約神威機器性能的發(fā)揮.因此,我們希望在神威上優(yōu)化通信延遲.

SW26010 眾核處理器采用片上融合的異構(gòu)體系結(jié)構(gòu)[15],包含4 個核組,如圖7所示.4 個核組之間通過片上互聯(lián)網(wǎng)絡(luò)(network on chip,簡稱Noc)連接,片上互聯(lián)網(wǎng)絡(luò)與系統(tǒng)接口(system i nterface,簡稱SI)相連,系統(tǒng)接口連接到以太網(wǎng).由申威處理器4 個核組與以太網(wǎng)連接的方式可以看出:4 個核組在硬件上共享同一個網(wǎng)絡(luò)接口,系統(tǒng)的通信能力固定.如果4 個核組的260 個核并發(fā)通信,多核并發(fā)通信將造成進程間競爭網(wǎng)絡(luò)資源;同時,較多的消息個數(shù)將導(dǎo)致網(wǎng)絡(luò)阻塞,通信性能降低.因此,本文考慮選用一個主核負責(zé)與外界通信,芯片內(nèi)部采用核間同步策略來降低競爭成本,實現(xiàn)網(wǎng)絡(luò)通信消息個數(shù)的減少,優(yōu)化通信延遲.

Fig.7 Connection of network on SW26010圖7 申威芯片網(wǎng)絡(luò)連接示意圖

分子動力學(xué)模擬中,粒子的狀態(tài)迭代函數(shù)需要獲取其遠程鄰居粒子的狀態(tài)信息,如更新后的位置,以進行新的受力計算,進而更新加速度、速度和位移.迭代函數(shù)的性質(zhì)決定了不同處理器之間的計算是相互依賴的,需要大量頻繁的數(shù)據(jù)同步.本文采用的單進程通信能夠降低處理器內(nèi)部的通信開銷,以及維護大量通信鏈接對資源的占用而帶來性能優(yōu)化.這部分的優(yōu)勢是因為分子動力學(xué)模擬程序是時間演化類以及高頻迭代類程序的典型代表,不同處理器核間頻繁大量的數(shù)據(jù)交換導(dǎo)致此類程序?qū)ρ舆t較敏感,網(wǎng)絡(luò)延遲對算法的性能影響較大.因此,神威系統(tǒng)上單進程通信的優(yōu)化算法主要針對對異構(gòu)架構(gòu)下類似的延遲敏感及通信受限類應(yīng)用的優(yōu)化提供了參考.

2.1 減少消息個數(shù),優(yōu)化通信延遲

2.1.1 片上同步減少消息個數(shù),但帶來訪存延遲

SW26010 芯片由4 個核組構(gòu)成,以核組為單位的編程模式中,每個核組的64 個從核在計算結(jié)束后,需要通過其所在主核與其他進程通信,交換粒子更新后的狀態(tài)信息,以進入下一個迭代步.由SW26010 芯片的結(jié)構(gòu)特點可知:4 個核組通過系統(tǒng)接口SI 共享同一個網(wǎng)絡(luò)端口,系統(tǒng)的通信能力固定.那么,如何優(yōu)化通信延遲成為首要解決的問題.

本文通過減少用于通信的進程數(shù)以減少消息個數(shù),優(yōu)化通信延遲.每個節(jié)點選取0 號核組即0 號進程作為通信進程,CPU 上其余3 個進程與0 號進程間通過片上同步進行通信.與以核組為單位的編程模式相比,通信進程數(shù)減少到原來的四分之一.

具體做法如圖8所示,每個節(jié)點的0 號主核首先進入通信等待狀態(tài),等待接收其他節(jié)點在上一時刻所更新的粒子信息.接收消息后,0 號主核需要將消息同步給本節(jié)點的所有從核以開始計算.由于神威處理器芯片是主從核異構(gòu)模式,該同步過程包括兩部分:主核同步和從核同步.0 號主核將消息同步給其他3 個主核,此時,4 個主核完成同步.每個主核需要將消息同步給該核組的64 個從核,實現(xiàn)主從核的同步.至此,所有負責(zé)計算的256 個從核得到了其他節(jié)點的粒子信息,開始計算.由于主從核的異步執(zhí)行,此時4 個核組的0 號主核等待本核組的所有從核計算結(jié)束.核組內(nèi)的從核計算結(jié)束后,通過與該核組的主核執(zhí)行一次主從核同步,使得該核組的主核獲得本核組的所有從核已計算完成的消息.4 個核組執(zhí)行一次主核同步,同步所有核組計算完成的狀態(tài).0 號主核將計算后的粒子狀態(tài)通信給其他節(jié)點.至此,一個完整的迭代步結(jié)束.

Fig.8 Process of synchronization on the chip圖8 片上同步的實現(xiàn)流程

異構(gòu)體系結(jié)構(gòu)下代碼分為兩部分:主核代碼和從核代碼,如圖9、圖10所示.

圖8 表示一個完整的迭代步,包括計算部分和通信部分.一個迭代步內(nèi),需要兩次主核同步、兩次主核與從核同步.由于神威處理器沒有cache,同步均需通過共享內(nèi)存數(shù)據(jù)的方式來實現(xiàn)[16],即同步需訪存.掣肘于神威處理器訪存長延遲的特點,節(jié)點同步雖減少了網(wǎng)絡(luò)通信的消息個數(shù),但同時也帶來了訪存延遲,對優(yōu)化消息個數(shù)帶來了挑戰(zhàn).那么,如何優(yōu)化訪存延遲?

申威處理器上內(nèi)存總大小為32GB,每個核組通過內(nèi)存控制器與8GB 的本核組內(nèi)存相連,如圖7所示,不同核組之間如需同步,需要訪問片外存儲器.那么,多核組之間共享內(nèi)存空間在內(nèi)存分配上將帶來一定的性能損失[16].該部分的性能損失在本文的優(yōu)化算法中,采用節(jié)點共享變量的方式規(guī)避.多核組間通過檢測節(jié)點共享變量實現(xiàn)同步,而非顯式同步的方式優(yōu)化訪存延遲.在節(jié)點主進程上設(shè)置節(jié)點共享變量,不同核組間通過共享變量以獲取所有核組間的同步狀態(tài),以此規(guī)避了對大內(nèi)存的頻繁訪問而帶來的訪存延遲.

Fig.9 Code on MPE for node-sync algorithm圖9 片上同步的主核代碼

Fig.10 Code on CPEs for node-sync algorithm圖10 片上同步的從核代碼

2.1.2 共享內(nèi)存等待與核間同步相結(jié)合,優(yōu)化訪存延遲

神威處理器主從核異構(gòu)的體系結(jié)構(gòu)特點,導(dǎo)致主從核之間的同步需要通過訪問片外存儲器實現(xiàn).多核訪存需要經(jīng)過訪存隊列,這種同步方式較慢,帶來訪存延遲.可見,片上同步的主要瓶頸在于同步需訪存.神威體系結(jié)構(gòu)下,主核間同步可通過節(jié)點共享變量實現(xiàn),單個核組的從核間同步速度較快,無需訪存.那么,如何實現(xiàn)所有主從核間的快速同步?我們采用共享內(nèi)存等待與核間同步相結(jié)合的方式,具體做法如圖11所示.

Fig.11 Process of combining shared memory waiting with inter-core synchronization圖11 共享內(nèi)存等待與核間同步相結(jié)合實現(xiàn)流程

流程如下:每個節(jié)點的0 號主核負責(zé)通信,將通信完成指示變量定義為節(jié)點共享變量,節(jié)點內(nèi)所有主核同時可見.那么節(jié)點內(nèi)的主核間可通過節(jié)點共享變量進行同步,無需顯式同步.此時,每個核組的0 號從核等待通信完成指示變量的變化.當(dāng)該變量加1 后,說明主核已完成通信,接收到其他節(jié)點在t時刻的狀態(tài)信息.0 號從核此時需要將此狀態(tài)同步給所在核組的其余63 個從核.同一核組的從核間同步速度較快,無需訪存.每個核組的從核完成組內(nèi)同步后,從核開始并行加速計算.核組內(nèi)所有從核完成計算后,該核組內(nèi)的64 個從核經(jīng)過一次從核同步,以保證64 個從核狀態(tài)統(tǒng)一.從核同步后,每個核組的0 號從核將計算完成的節(jié)點共享變量加1.該共享變量對所有核組的主核同時可見.所有從核加速計算的同時,0 號通信主核等待4 個0 號從核計算完成的節(jié)點共享變量的變化.當(dāng)0 號主核檢測到4 個0 號從核均將計算完成指示變量加1 后,0 號主核開始與其他節(jié)點進行通信.至此,完成一個完整的迭代步.主從核的代碼實現(xiàn)參見圖12、圖13 的虛代碼所示.

Fig.12 Code on MPE for shared memory waiting圖12 共享內(nèi)存等待算法的主核代碼

Fig.13 Code on CPEs for shared memory waiting圖13 共享內(nèi)存等待算法的從核代碼

神威異構(gòu)體系結(jié)構(gòu)下,單核組通信的模式需利用核間同步,以同步通信數(shù)據(jù).SW26010 芯片的硬件結(jié)構(gòu)決定了主從核之間同步需要訪問片外存儲器,申威處理器較長的訪存延遲導(dǎo)致該同步模式效率較低,本文采用共享內(nèi)存等待的方式規(guī)避訪存延遲.共享內(nèi)存等待方式的優(yōu)勢是:核組間無需顯式同步,只需要設(shè)置共享內(nèi)存等待變量;該變量對節(jié)點內(nèi)所有主核同時可見,只需要檢測共享內(nèi)存指示變量的變化來實現(xiàn)不同核組間的隱式同步,避免了核組同步時的訪存操作,規(guī)避了申威處理器較長的訪存時延.

2.2 改變計算模式,減少核間等待,優(yōu)化通信延遲

2.2.1 基本的并行模式

共價鍵分子的原子間相互作用勢不僅取決于原子間的距離,與原子間的成鍵方向有密切關(guān)系,因此需要考慮周圍原子的影響.硅原子之間的多體作用模型采用Tersoff 勢函數(shù)[14],計算涉及兩個以上的原子.計算原子i受到它的鄰居j的作用力時,需要考慮i的其他3 個鄰居此時對i的影響.見圖14:由Tersoff 勢函數(shù)可得到鄰居j對i的作用力func1,力是相互作用力,原子j也受到了一個大小相等方向相反的作用力?func1,并且需要加到j(luò)原子上[17].在神威上是多線程計算,j粒子同時在由另外一個線程計算j粒子與它的鄰居粒子之間的受力.這個方向相反的作用力?func1 在同一時刻不能同時加到j(luò)粒子身上.基本算法是:需要存儲本線程所計算的對應(yīng)原子受到其所有鄰居原子的全部作用力,通過同步,傳遞給對應(yīng)原子所在線程.所有鄰居原子線程的受力計算結(jié)束后,接收消息,將原子受到的鄰居粒子的反作用力疊加起來,保證所有原子作用力的計算是完整的、準(zhǔn)確的.

對每個原子i而言,以計算它的鄰居粒子j=0 對它的作用力func1 為例,說明受力計算的過程.首先,計算鄰居粒子j對i的作用力時,需要考慮i粒子其余3 個鄰居粒子k0,k1,k2此時對i的影響.func1 是i和j之間大小相等、方向相反的相互作用力,即i粒子對j粒子的反作用力?func1,需同時累加到j(luò)粒子上.若為并行程序,此時j粒子正在由其他線程計算j粒子與它的鄰居粒子之間的受力.不能同時把?func1 累加到j(luò)粒子上,需要交互粒子信息.最后,其他的3 個鄰居粒子也受到了反作用力?func2,同樣需要在3 個粒子的受力上減去.受力計算結(jié)束后,算法采用蛙跳格式更新位移.

Fig.14 Inter-atomic forces圖14 原子間作用力

在分子動力學(xué)模擬中,分子間作用力的計算時間占80%以上[18,19],是最重要的部分,也是發(fā)揮神威機器計算性能的關(guān)鍵所在.在迭代過程中,計算核心需要頻繁地交互粒子信息,以保證分子間作用力計算的準(zhǔn)確性.神威機器高延遲的特點,將使粒子信息的交互遇到瓶頸,此類應(yīng)用會限制神威性能的發(fā)揮.程序特點是典型的BSP 模式,如圖15所示.每個迭代步將粒子的位置,速度等信息存儲在一個向量中,該向量包含多個分量.粒子間作用力的計算只涉及到其鄰居粒子,向量的不同分量間具有局部相關(guān)性.即:后面的分量是由前面幾個分量計算而來,僅與前面的幾個分量有關(guān).數(shù)據(jù)的局部相關(guān)性和依賴性,導(dǎo)致具有一定的異步性.

Fig.15 Implementation process of parallel algorithm圖15 并行實現(xiàn)算法流程圖

2.2.2 改變計算模式,減少通信次數(shù),優(yōu)化消息個數(shù)

在原來的計算模式中,向量的一個分量在計算過程中會產(chǎn)生一些中間數(shù)據(jù),另外的分量可以利用.向量的不同分量之間是相互依賴的關(guān)系,是緊耦合的,導(dǎo)致芯片內(nèi)部通過主存的這種數(shù)據(jù)傳輸和數(shù)據(jù)交換會大量增加.我們打破了這種依賴,變?yōu)樗神詈?把原來互相寫同步的計算模式變成了每個線程自己去做多步的計算,是獨立計算,減少迭代步中間的同步次數(shù),更流水化.代價是會造成計算量少量增加,但是整個吞吐率反而提高了.

多體作用模型中,如圖16所示,計算i粒子受到鄰居j的作用力時,其他3 個鄰居k0,k1,k2也對i粒子產(chǎn)生了作用力.這意味著j,k0,k1,k2分別受到了來自i粒子的反作用力,需要從它們的受力上減去.考慮i粒子的全部受力,不僅需要計算它受到其鄰居粒子的作用力,還需要計算i粒子作為鄰居粒子時,它所受到的反作用力.如圖17所示,當(dāng)i粒子作為j的鄰居粒子時,它會受到j(luò)粒子對它的反作用力?func3.此外,如圖18所示,j粒子在計算它的鄰居k對其作用力時,會對i產(chǎn)生一個反作用力?func4,也需要從i粒子的受力上減去.

具體做法如下:

(1)計算i粒子受到的作用力.i粒子所受到的作用力來自于它的所有鄰居粒子.虛代碼如圖19所示.變量j循環(huán)i粒子的所有鄰居,例如,i粒子的第0 個鄰居j=0 時,計算j對i的作用力func1.變量k循環(huán)i粒子除了j以外的所有鄰居,它們對i粒子的力為func2,不同的鄰居,如k0,k1,k2對i粒子的力不同,與位置有關(guān);

(2)計算i粒子受到的反作用力.i粒子所受到的反作用力來自于以它的鄰居粒子為主的計算過程中,如i粒子作為j粒子的鄰居時,i粒子所受到反作用力.計算j粒子受力時,遍歷j的所有鄰居.假設(shè)j的0 號鄰居為i,計算得到i作為j的鄰居時,i受到的反作用力func3,需要從i粒子的受力上減去.atom[i].f?=func3.若0 號鄰居不為i,j對i的反作用力func4,atom[i].f?=func4.

Fig.16 Force of i atom圖16 i 粒子受到的作用力

Fig.17 Reaction force func3圖17 i 粒子受到的反作用力func3

Fig.18 Reaction force func4圖18 i 粒子受到的反作用力func4

Fig.19 Algorithm for changing the calculation mode圖19 改變計算模式的算法

i粒子受到的反作用力在原算法中是由負責(zé)計算它的鄰居粒子受力的線程完成,如func3 和func4 是由負責(zé)計算j粒子的線程完成,那么線程間會產(chǎn)生相互依賴.在新的算法中,i粒子受到的反作用力仍由負責(zé)計算i粒子受力的線程計算.即將互相寫同步的緊耦合計算模式改為獨立計算的松耦合模式.每個線程在計算粒子受力時獨立計算,通信僅發(fā)生在迭代計算完成后,減少了通信次數(shù)和核間等待,優(yōu)化了通信延遲.

硅粒子結(jié)晶的多體作用模型中,原計算模式下,粒子i只需要負責(zé)計算它的4 個鄰居粒子對它的作用力,其他粒子對i粒子的反作用力通過進程間同步而得.優(yōu)化后的計算模式為:減少粒子i與其鄰居粒子間的同步次數(shù),由i粒子完成反作用力部分的計算.所增加的計算量是i粒子需計算其鄰居粒子對它的所有反作用力.以i粒子受到其鄰居粒子j的反作用力為例,當(dāng)i粒子作為j粒子的第1 個鄰居粒子時,i粒子受到j(luò)對i的反作用力func3.當(dāng)i粒子作為j粒子的第2 個、第3 個或第4 個鄰居粒子時,i粒子受到j(luò)對i的反作用力為func4.那么func3和func4 的值需要從i粒子的受力計算中減去.該部分的計算與原計算模式相比為增加的部分.計算模式的改變打破了數(shù)據(jù)依賴,將緊耦合變?yōu)樗神詈?雖然帶來了計算量的增加,但是提高了吞吐.

本文提出的計算模式的優(yōu)化策略,適用于類似的時間演化類或高頻迭代程序,當(dāng)問題本身的結(jié)構(gòu)特性為數(shù)據(jù)具有局部相關(guān)性和依賴性,由此可產(chǎn)生一定的異步性的性質(zhì),可利用本文提出的優(yōu)化計算模式的方式,將不同數(shù)據(jù)之間緊耦合的依賴關(guān)系變?yōu)樗神詈?將原計算模式中互相寫同步的計算模式改為更流水化的獨立計算.該部分的計算對計算強度影響的主要原因是,由于減少了進程間同步操作而避免了受力計算的同步更新,因此需要在本地進程額外替其他進程執(zhí)行相應(yīng)的計算,雖然計算強度少量增加,但該算法利用了處理器較強的計算性能,掩蓋了因計算強度增加而帶來的性能損耗,提高了吞吐,優(yōu)化了性能.

2.3 高效訪存及規(guī)則化

2.3.1 高效訪存

在分子間作用力的計算過程中,從核需要頻繁地訪問粒子信息,從核訪問主存的延遲極大地限制了神威性能的發(fā)揮.每個從核上提供了64KB 的高速緩存空間SPM,在算法實現(xiàn)上,需要頻繁訪問的粒子信息通過DMA方式批量傳送至SPM,以優(yōu)化從核訪主存效率.

每個線程負責(zé)N個粒子的計算,線程i將所負責(zé)的所有原子和其鄰居原子的速度、位移等信息從主存加載到SPM 中,如圖20所示.迭代計算前,將粒子信息通過同步DMA 加載到SPM 中,m_memcpy(MEM_TO_CACHE),將所有粒子的鄰居粒子信息通過異步DMA 加載到SPM 中,m_memcpy_async(MEM_TO_CACHE);迭代計算完成后,通過異步DMA 方式,將本線程所負責(zé)粒子的更新信息從SPM 拷貝回系統(tǒng)內(nèi)存,m_memcpy_async(CACHE_TO_MEM).異步DMA 方式將計算與主存訪問相重疊.虛代碼如圖21所示.

Fig.20 SPM and main memory exchange data in DMA mode圖20 SP M 與主存間以DMA 方式交互數(shù)據(jù)

Fig.21 Algorithm for efficient accessing memory圖21 高效訪存算法

2.3.2 規(guī)則化

分子動力學(xué)模擬中,以硅結(jié)晶過程為例.問題本身不一定是規(guī)則的,但是可以以盡量規(guī)則化的方法進行計算,在數(shù)據(jù)結(jié)構(gòu)上盡量規(guī)則化.所謂規(guī)則化,就是按照向量長度是對齊,連續(xù),大小規(guī)整[20].規(guī)整是按照硬件結(jié)構(gòu)來規(guī)整的,具體做法是:將數(shù)據(jù)結(jié)構(gòu)規(guī)則化,向量寬度規(guī)則化.程序中的主要數(shù)據(jù)結(jié)構(gòu)是存儲例子位移、速度、加速度等信息,有一些非規(guī)則的部分.由于神威芯片的核是向量化的、規(guī)則化的,即硬件是規(guī)則化的,我們在設(shè)計類型和數(shù)據(jù)結(jié)構(gòu)的時候,將數(shù)據(jù)規(guī)則存儲.由于芯片的向量化、規(guī)則化,因而規(guī)則化數(shù)據(jù)結(jié)構(gòu)以提升性能.本文考慮常溫下,晶體硅的模擬.硅晶體為規(guī)則的晶胞結(jié)構(gòu),有4 個鄰居粒子.受力計算中,需分別計算4 個鄰居粒子b,c,d,e與a粒子的距離,如圖22所示.

Fig.22 Distance between atoms圖22 原子間的距離

4 次距離的計算是不凝聚的,我們針對神威核向量化的特點,考慮指令優(yōu)化,進行規(guī)則化計算.三維距離是計算兩點間3 個維度上的差的平方和再開方,即:

由硅粒子的規(guī)則晶胞結(jié)構(gòu)可知,a粒子的受力需計算4 個鄰居粒子與a粒子的距離.我們利用規(guī)則化的數(shù)據(jù)結(jié)構(gòu),首先將粒子三維坐標(biāo)的數(shù)據(jù)結(jié)構(gòu)規(guī)則化.神威的并行C 語言上擴展的數(shù)據(jù)類型為doublev4,由4 個雙精度浮點數(shù)構(gòu)成.粒子坐標(biāo)即定義為doublev4 類型的數(shù)據(jù)結(jié)構(gòu),如doublev4a.s;表示a粒子的位置,它的4 個浮點數(shù)分別代表三維坐標(biāo)和0.規(guī)則化的數(shù)據(jù)結(jié)構(gòu)去計算a與b,c,d,e之間距離時,每兩點間需進行差平方后的求和操作.如圖23所示,a與b間的距離首先需計算公式(12):

Fig.23 Vector representation of the distance between atoms圖23 原子間的距離的向量表示

若順序執(zhí)行4 次這樣的求和指令,那么這4 次計算是不凝聚的.考慮規(guī)則化計算,以提高訪存凝聚性.規(guī)則化計算的方式進行指令優(yōu)化,通過一個求和指令,實現(xiàn)4 次求和計算.規(guī)則化的向量化求和是計算4 個向量之和,而不是向量的4 個分量之和.因而考慮將原始的4 個向量進行轉(zhuǎn)置,實現(xiàn)規(guī)則化計算.修改后的形式如圖24所示.

Fig.24 Distance representation after transposition圖24 轉(zhuǎn)置后的距離表示

由此可見:對轉(zhuǎn)置后的4 個向量進行向量化求和后再開平方,得到的向量的4 個分量分別表示4 個粒子與a粒子的歐幾里得距離.由此可見:向量化的SIMD 指令等價于一個4 次循環(huán)操作,減少了指令數(shù),降低了對指令訪問帶寬的要求,而且減少了循環(huán)引起的控制相關(guān),提升了訪存凝聚性,提高了效率.

3 實驗結(jié)果

3.1 單節(jié)點性能

我們的實例是基于硅結(jié)晶的模擬,實驗設(shè)計上,我們逐步對比了不同的優(yōu)化方法對性能的影響,以及不同問題規(guī)模下性能的差異.在這一部分,我們主要闡述了在第2 節(jié)中所提出的3 種優(yōu)化技術(shù)對性能的提升效果.為了說明可擴展性,我們測試了兩種不同的問題規(guī)模,分別是4 096 個粒子和32 768 個粒子.我們的測試是運行在單節(jié)點上,啟用4 個核組.表1 給出了神威的系統(tǒng)配置.

Table 1 Su nway TaihuLight supercomputer system configuration表1 神威太湖之光超級計算機的系統(tǒng)配置

圖25 描述了每一步采用不同的優(yōu)化方法帶來的加速比的提升,我們的優(yōu)化方法最終能夠?qū)崿F(xiàn)18x 的加速比.共享內(nèi)存等待與從核同步相結(jié)合的方式對加速比的提升影響較大,緊耦合的計算模式、訪存優(yōu)化及規(guī)則化對性能提升也有一定的影響.從性能提升的角度,我們減少了消息個數(shù),優(yōu)化了通信以及訪存延遲,對時間演化類問題迭代頻率的提升做出了較大貢獻.

Fig.25 Speedup of different optimization methods圖25 不同優(yōu)化方法的加速比

3.2 多節(jié)點性能

我們的實驗平臺是基于神威太湖之光高性能集群.為測試程序的強可擴展性,我們首先測試了單個CPU 的性能.單個SW26010 處理器的性能取決于每個從核上所負責(zé)計算的原子個數(shù),我們測試了每個從核上處理不同數(shù)目的原子時性能的變化,見表2.模擬速度最大可以達到120KSteps/s.每個CPU 上所負責(zé)的原子個數(shù)最大為98K,浮點利用率達到15%.

Table 2 Single SW26010 processor performance表2 單SW26010 處理器的性能

程序的弱可擴展性測試是基于相同的數(shù)據(jù)規(guī)模,節(jié)點個數(shù)不同時,考察系統(tǒng)的性能.強可擴展性的測試是基于不同的粒子數(shù),擴展系統(tǒng)的節(jié)點數(shù),節(jié)點數(shù)目隨著粒子數(shù)的增加而增長.我們的測試是針對典型的細粒度分布,來考察典型計算模式下,系統(tǒng)的速率及強可擴展性.

我們的測試利用了SW26010 處理器的所有4 個核組,每個處理器運行4 個進程,但負責(zé)通信的進程只有0號主進程.為了驗證程序的強可擴展性,我們考察了不同規(guī)模的節(jié)點個數(shù),從8 個節(jié)點到32 768 個節(jié)點.

每個節(jié)點4 個進程,最多可達到131K 個進程.不同的系統(tǒng)規(guī)模下,每個進程上負責(zé)處理的原子個數(shù)固定不變,為512 個原子.我們考察這種典型計算模式下系統(tǒng)的運算速率.

我們的實驗為考察程序的強可擴展性,隨著進程數(shù)的增加,粒子數(shù)同比例增長,即每個進程上所負責(zé)計算的粒子數(shù)目保持不變.在上述細粒度分布的條件下,我們得到這種典型計算模式下的速率,即每秒迭代的步數(shù),見表3 第4 列.

節(jié)點數(shù) 進程數(shù) 粒子數(shù) Steps/s us/step Time

Table 3 Comparison of software solution on Sunway表3 神威上不同系統(tǒng)規(guī)模下的對比

我們的測試通過改變系統(tǒng)規(guī)模,即進程數(shù)來觀察系統(tǒng)性能的變化,如圖26所示.縱坐標(biāo)是每個迭代步所消耗的時間,橫坐標(biāo)表示了進程數(shù)的不同,在每個點上的數(shù)據(jù),對應(yīng)了不同的粒子數(shù).

Fig.26 Time spent on each iteration step at different system sizes圖26 不同系統(tǒng)規(guī)模下每個迭代步的耗時

同時,我們還測試了在不同系統(tǒng)規(guī)模下,程序總耗時的變化,如圖27所示.在圖27 的折線圖中,標(biāo)出了每個進程點所對應(yīng)的粒子個數(shù),可觀察到時間的變化.

Fig.27 Time changes with different scales圖27 不同系統(tǒng)規(guī)模下時間的變化

我們的目標(biāo)是提高通信受限類程序的迭代頻率,為了比較系統(tǒng)的實際性能,我們采用的度量公式如下:

performance=natoms×nsteps/sec.

在不同的系統(tǒng)規(guī)模下,對應(yīng)的粒子個數(shù)也不同,綜合對比系統(tǒng)性能.如圖28所示,系統(tǒng)性能隨著進程數(shù)的增加而穩(wěn)步提升.

Fig.28 Performance with different system sizes圖28 不同系統(tǒng)規(guī)模下的性能

由表4 可知,本文給出的優(yōu)化策略高于其他軟件解決方案.當(dāng)原子個數(shù)小于10million 時,迭代速度大于10Ksteps/s,該迭代速度高于現(xiàn)有的絕大部分軟件解決方案.同時,較Anton 硬件解決方案的優(yōu)勢是它的可擴展性,系統(tǒng)規(guī)?？蛇_到幾萬個處理器并行執(zhí)行,同時可模擬的原子個數(shù)可達到5 千萬以上.這對大規(guī)模的分子動力學(xué)模擬提供了有價值的參考.

Table 4 Performance comparisons of hardware and software solutions for MD表4 分子動力學(xué)模擬軟硬件解決方案的性能對比

4 相關(guān)工作

計算密集型和訪存密集型程序在神威“太湖之光”上得以優(yōu)化.Stencil 問題具有較高的計算吞吐,在神威上實現(xiàn)計算-通信重疊,優(yōu)化通信開銷[25].優(yōu)化神威上HPCG 算法中的有效內(nèi)存帶寬以及增強算法的可擴展性[9].GTC-P 大規(guī)模并行模擬的高性能計算程序針對神威的訪存帶寬進行優(yōu)化[26].

神威“太湖之光”超級計算機強大的運算能力,使其能夠處理多種大規(guī)模的應(yīng)用.在神威高性能集群上實現(xiàn)了超大規(guī)模的氣象模擬[27].大規(guī)模非線性地震模擬[28]針對神威體系結(jié)構(gòu)特點給出并行化解決方案.此類應(yīng)用的特點是數(shù)據(jù)規(guī)模龐大,針對內(nèi)存空間和帶寬給出了優(yōu)化方案.

時間演化類應(yīng)用旨在解決的問題是提高迭代頻率,加速時間演化過程.Anton[21,22]機器是針對分子動力學(xué)模擬設(shè)計的一款專用目的計算機,硬件上設(shè)計的低延遲、高帶寬特點的網(wǎng)絡(luò)以用于快速同步,但是限制了系統(tǒng)的物理規(guī)模.在神威上,從數(shù)據(jù)的預(yù)取和向量化角度優(yōu)化[29]LAMMPS 中對內(nèi)存數(shù)據(jù)的訪問.針對計算密集型的GROMACS 程序,在神威“太湖之光”上解決內(nèi)存帶寬限制的問題[30].

由于時間演化類應(yīng)用本身數(shù)據(jù)依賴性的特點,不同處理器間的頻繁通信將極大地制約迭代頻率的提高.本文以減少延遲敏感的時間演化類程序的通信延遲為主要目標(biāo),優(yōu)化通信,并提出幾種有效的并行化策略,為類似的通信受限類程序在異構(gòu)的國產(chǎn)化神威機器上的應(yīng)用提供了藍本.

5 總結(jié)

在本文中,我們實現(xiàn)了分子動力學(xué)模擬程序在神威太湖之光超級計算機上的優(yōu)化.我們的實現(xiàn)是基于以核組為單位的編程模式,在系統(tǒng)規(guī)模和網(wǎng)絡(luò)通信能力不變的前提下,利用片上同步,減少了消息個數(shù),優(yōu)化了通信延遲.通過共享內(nèi)存等待與從核同步相結(jié)合的方式,進一步優(yōu)化了片上同步帶來的訪存延遲.同時,我們針對分子間多體作用力的計算模式進行修改,將互相寫同步的緊耦合計算模式改為松耦合.減少了迭代步中間的同步次數(shù),打破了不同線程間的依賴關(guān)系,提高了吞吐.此外,進行了訪存優(yōu)化以及規(guī)則化數(shù)據(jù)結(jié)構(gòu)以提高訪存凝聚性.我們的工作是針對諸如分子動力學(xué)模擬等延遲敏感的時間演化類應(yīng)用如何提高迭代頻率,給出的一系列優(yōu)化技術(shù),為類似的通信受限類程序在主從核異構(gòu)的國產(chǎn)神威處理器上的優(yōu)化提供了參考.今后的工作中,我們將進一步探索神威上的優(yōu)化技術(shù),對時間演化類程序進行高效模擬.