田 卓,陳一峯

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

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

時間演化類問題從時間軸來模擬體系演化.時間被離散化為一定步長為間隔的時間點,下一個時間點的狀態(tài)依賴上一個時間點的狀態(tài),兩個狀態(tài)之間需要一次通信以交換數(shù)據(jù).典型的應用是分子動力學模擬[5,6],分子沿著這些時間點運動,每個時間點對應分子的一個狀態(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)間的關系是相互依賴的,需要通信.而神威上的通信延遲較長[7],包括網(wǎng)卡上的通信延遲以及較長的數(shù)據(jù)通路所帶來的通信延遲都太長了.此類計算模式的時間演化類程序在高延遲的神威集群上如何優(yōu)化,規(guī)避通信延遲,提高性能,是本文的宗旨.

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

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

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

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

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

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

1 神威結構及計算實例

1.1 神威結構

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

Fig.4 Architecture of the SW26010 processor圖4 神威體系結構

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

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

1.2 計算實例

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

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

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

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

考慮常溫下,硅晶體保持規(guī)則的金剛石晶胞結構,每個原子周圍有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)實的.為了提高硅原子模擬的迭代頻率,有效的做法是減少迭代步間通信的延遲.由神威體系高性能集群的特點可知,其上的通信延遲較長.那么針對諸如分子動力學模擬的時間演化類程序,如何在神威上減少通信延遲、提高迭代頻率,本文給出了一系列優(yōu)化技術.

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

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

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)卡上和每臺服務器上的延遲都太長了.對通信受限的時間演化類程序而言,通信延遲將極大地制約神威機器性能的發(fā)揮.因此,我們希望在神威上優(yōu)化通信延遲.

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

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

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

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

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

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

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

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

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

異構體系結構下代碼分為兩部分:主核代碼和從核代碼,如圖9、圖10所示.

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

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

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

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

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

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

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

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

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

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

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

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

2.2.1 基本的并行模式

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

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

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

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

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

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

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

多體作用模型中,如圖16所示,計算i粒子受到鄰居j的作用力時,其他3 個鄰居k0,k1,k2也對i粒子產生了作用力.這意味著j,k0,k1,k2分別受到了來自i粒子的反作用力,需要從它們的受力上減去.考慮i粒子的全部受力,不僅需要計算它受到其鄰居粒子的作用力,還需要計算i粒子作為鄰居粒子時,它所受到的反作用力.如圖17所示,當i粒子作為j的鄰居粒子時,它會受到j粒子對它的反作用力?func3.此外,如圖18所示,j粒子在計算它的鄰居k對其作用力時,會對i產生一個反作用力?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粒子的力不同,與位置有關;

(2)計算i粒子受到的反作用力.i粒子所受到的反作用力來自于以它的鄰居粒子為主的計算過程中,如i粒子作為j粒子的鄰居時,i粒子所受到反作用力.計算j粒子受力時,遍歷j的所有鄰居.假設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粒子受到的反作用力在原算法中是由負責計算它的鄰居粒子受力的線程完成,如func3 和func4 是由負責計算j粒子的線程完成,那么線程間會產生相互依賴.在新的算法中,i粒子受到的反作用力仍由負責計算i粒子受力的線程計算.即將互相寫同步的緊耦合計算模式改為獨立計算的松耦合模式.每個線程在計算粒子受力時獨立計算,通信僅發(fā)生在迭代計算完成后,減少了通信次數(shù)和核間等待,優(yōu)化了通信延遲.

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

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

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

2.3.1 高效訪存

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

每個線程負責N個粒子的計算,線程i將所負責的所有原子和其鄰居原子的速度、位移等信息從主存加載到SPM 中,如圖20所示.迭代計算前,將粒子信息通過同步DMA 加載到SPM 中,m_memcpy(MEM_TO_CACHE),將所有粒子的鄰居粒子信息通過異步DMA 加載到SPM 中,m_memcpy_async(MEM_TO_CACHE);迭代計算完成后,通過異步DMA 方式,將本線程所負責粒子的更新信息從SPM 拷貝回系統(tǒng)內存,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ī)則化

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

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

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

由硅粒子的規(guī)則晶胞結構可知,a粒子的受力需計算4 個鄰居粒子與a粒子的距離.我們利用規(guī)則化的數(shù)據(jù)結構,首先將粒子三維坐標的數(shù)據(jù)結構規(guī)則化.神威的并行C 語言上擴展的數(shù)據(jù)類型為doublev4,由4 個雙精度浮點數(shù)構成.粒子坐標即定義為doublev4 類型的數(shù)據(jù)結構,如doublev4a.s;表示a粒子的位置,它的4 個浮點數(shù)分別代表三維坐標和0.規(guī)則化的數(shù)據(jù)結構去計算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 個向量進行轉置,實現(xiàn)規(guī)則化計算.修改后的形式如圖24所示.

Fig.24 Distance representation after transposition圖24 轉置后的距離表示

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

3 實驗結果

3.1 單節(jié)點性能

我們的實例是基于硅結晶的模擬,實驗設計上,我們逐步對比了不同的優(yōu)化方法對性能的影響,以及不同問題規(guī)模下性能的差異.在這一部分,我們主要闡述了在第2 節(jié)中所提出的3 種優(yōu)化技術對性能的提升效果.為了說明可擴展性,我們測試了兩種不同的問題規(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)化方法最終能夠實現(xiàn)18x 的加速比.共享內存等待與從核同步相結合的方式對加速比的提升影響較大,緊耦合的計算模式、訪存優(yōu)化及規(guī)則化對性能提升也有一定的影響.從性能提升的角度,我們減少了消息個數(shù),優(yōu)化了通信以及訪存延遲,對時間演化類問題迭代頻率的提升做出了較大貢獻.

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

3.2 多節(jié)點性能

我們的實驗平臺是基于神威太湖之光高性能集群.為測試程序的強可擴展性,我們首先測試了單個CPU 的性能.單個SW26010 處理器的性能取決于每個從核上所負責計算的原子個數(shù),我們測試了每個從核上處理不同數(shù)目的原子時性能的變化,見表2.模擬速度最大可以達到120KSteps/s.每個CPU 上所負責的原子個數(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 個進程,但負責通信的進程只有0號主進程.為了驗證程序的強可擴展性,我們考察了不同規(guī)模的節(jié)點個數(shù),從8 個節(jié)點到32 768 個節(jié)點.

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

我們的實驗為考察程序的強可擴展性,隨著進程數(shù)的增加,粒子數(shù)同比例增長,即每個進程上所負責計算的粒子數(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所示.縱坐標是每個迭代步所消耗的時間,橫坐標表示了進程數(shù)的不同,在每個點上的數(shù)據(jù),對應了不同的粒子數(shù).

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

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

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

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

performance=natoms×nsteps/sec.

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

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

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

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

4 相關工作

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

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

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

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

5 總結

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