郝 赫，司雨蒙，韋建文，文敏華，林新華,2+

1.上海交通大學(xué) 高性能計算中心，上海 200240

2.NVIDIA Technology Center Asia Pacific，Singapore 999002

天體物理成團研究中的非規(guī)則訪存優(yōu)化*

郝 赫1，司雨蒙1，韋建文1，文敏華1，林新華1,2+

1.上海交通大學(xué) 高性能計算中心，上海 200240

2.NVIDIA Technology Center Asia Pacific，Singapore 999002

HAO He,SI Yumeng,WEI Jianwen,et al.Optimizing irregular memory access in astrophysical clustering studies.Journal of Frontiers of Computer Science and Technology,2017,11(1)：80-90.

HGGF（halo-based galaxy group finder）算法實現(xiàn)了基于暗物質(zhì)暈的星系找群，在研究宇宙大尺度結(jié)構(gòu)及宇宙的演化等領(lǐng)域中占有至關(guān)重要的地位。但由于數(shù)據(jù)規(guī)模的增長，急需對HGGF算法進(jìn)行優(yōu)化，以縮短運行時間。經(jīng)分析，算法的熱點部分耗時受到非規(guī)則訪存的嚴(yán)重影響，因此針對算法的結(jié)構(gòu)和非規(guī)則訪存模型，提出了數(shù)據(jù)預(yù)排序方法，并分析了該方法如何影響訪存過程。在此基礎(chǔ)上，利用數(shù)據(jù)對齊、循環(huán)分解進(jìn)一步優(yōu)化訪存效率，利用負(fù)載均衡和互斥變量私有化的方法提高了OpenMP的并行效率，最終將HGGF應(yīng)用使用12線程加速11.6倍，同時取得了更好的可擴展性。主要有三點貢獻(xiàn)：（1）分析了HGGF算法的非規(guī)則訪存問題；（2）提出并分析了數(shù)據(jù)預(yù)排序方法；（3）使用數(shù)據(jù)對齊、循環(huán)分解、負(fù)載均衡、互斥變量私有化方法提高了HGGF應(yīng)用的并行性能。

天體物理成團；非規(guī)則訪存優(yōu)化；數(shù)據(jù)預(yù)排序；并行計算

1 引言

隨著天文觀測能力的不斷提高，天文數(shù)據(jù)急劇增加，其中星系的觀測總量已經(jīng)達(dá)到109量級。面對來自大型數(shù)字巡天計劃的海量數(shù)據(jù)，如何從數(shù)據(jù)中迅速準(zhǔn)確地提取所需要的內(nèi)容，直接影響著天文學(xué)的發(fā)展和研究進(jìn)程。其中，如何將可觀測的星系和理論中不可見的暗物質(zhì)暈聯(lián)系起來，是天體物理研究中的重要難題之一[1]。天文學(xué)家利用計算機設(shè)計出多種Group-Finding算法來尋找同一暗物質(zhì)暈中的星系。FOF（friends-of-friends）是Group-Finding算法中最早的算法之一[2]，對于理解宇宙的結(jié)構(gòu)和演化具有重大意義。HGGF（halo-based galaxy group finder）算法由FOF算法演變而來，雖然同樣是對星系進(jìn)行分類，但是HGGF算法是依賴暗物質(zhì)暈的性質(zhì)，能實現(xiàn)對某些孤立的星系和暗物質(zhì)暈的聯(lián)接[3]。HGGF算法還將亮度、紅移距離、暗暈質(zhì)量等性質(zhì)共同作為篩選的依據(jù)，相比僅依靠距離的算法更有實際意義和準(zhǔn)確性。

星系成團類算法基本步驟是預(yù)先指定星系中心，根據(jù)連接長度[4]逐個將星系樣本分配到某一星系群，再反復(fù)迭代直到收斂，使用的數(shù)據(jù)結(jié)構(gòu)多為基于指針的圖或樹，如FOF算法即基于K-D樹[5]實現(xiàn)。雖然整體框架基本相同，但是HGGF算法使用數(shù)組替代指針作為主要數(shù)據(jù)結(jié)構(gòu)，而在天文學(xué)或分子動力學(xué)的N體模擬問題中，粒子間的位置和連接關(guān)系往往非常復(fù)雜，使用數(shù)組的存儲效率更高[6]。在HGGF算法的核心部分，由于連接體積和索引數(shù)組[7]的存在，主數(shù)組的訪問順序被徹底打亂，因此將HGGF算法的訪存方式視為非規(guī)則訪存。非規(guī)則訪存算法通常表現(xiàn)出極少的數(shù)據(jù)局部性和數(shù)據(jù)重用，大量的動態(tài)非連續(xù)存儲訪問以及細(xì)粒度的同步或通訊[8]，在現(xiàn)代處理器基于緩存的架構(gòu)下，緩存的利用效率會對程序的性能起決定性作用。而非規(guī)則訪存算法由于訪問地址不連續(xù)，甚至是隨機訪問，緩存內(nèi)的數(shù)據(jù)不能得到高效利用，加上CPU和內(nèi)存之間的速度差異，即使CPU計算能力很強，也不能獲得好的性能。

隨著異構(gòu)計算的發(fā)展，在科學(xué)計算或是工程計算中加速卡的使用已經(jīng)非常普遍。人們嘗試將HGGF應(yīng)用移植到Intel公司的MIC（many integrated core）上，但是由于處理器的性能相比CPU差距較大[9]，HGGF算法沒有在MIC上取得理想的性能，需要進(jìn)一步優(yōu)化。

數(shù)據(jù)預(yù)取是通過優(yōu)化局部性改善訪存性能的常見方法。數(shù)據(jù)預(yù)取包括硬件預(yù)取和軟件預(yù)取，硬件預(yù)取一般由編譯器自動執(zhí)行，而軟件預(yù)取需要在代碼中插入明顯的預(yù)取指令[10]。我們使用了兩種軟件預(yù)取方法，但效果一般。經(jīng)過進(jìn)一步分析算法結(jié)構(gòu)，發(fā)現(xiàn)主數(shù)組訪存順序同時受到連接體積和索引數(shù)組影響，因此采用了按轉(zhuǎn)換坐標(biāo)對輸入數(shù)據(jù)預(yù)排序的方法，明顯改善了數(shù)據(jù)局部性，減少了緩存未命中的次數(shù)，在單核上加速超過2倍。同時，利用OpenMP在單節(jié)點上實現(xiàn)了HGGF算法中篩選星系過程的多線程并行，并且從負(fù)載均衡、同步開銷等方面對算法進(jìn)一步優(yōu)化，最終在12線程上獲得11.6倍的加速，在將數(shù)據(jù)規(guī)模擴大16倍后，獲得15倍加速。

綜上，本文主要有以下三點貢獻(xiàn)：

（1）分析了HGGF算法的非規(guī)則訪存問題；

（2）提出并分析了數(shù)據(jù)預(yù)排序方法；

（3）使用傳統(tǒng)方法提高了HGGF應(yīng)用的并行性能。

本文組織結(jié)構(gòu)如下：第2章介紹相關(guān)工作；第3章介紹HGGF算法及其訪存模型；第4章詳細(xì)介紹優(yōu)化方法；第5章是實驗環(huán)境及結(jié)果分析；最后總結(jié)全文，并提出下一步工作。

2 相關(guān)工作

目前，星系成團研究中對并行計算的使用還比較有限，現(xiàn)有工作局限在對FOF算法[11]、HOP算法[12]的研究中。文獻(xiàn)[11]介紹了卡內(nèi)基梅隆大學(xué)的Gibson等人利用Hadoop集群對超大天文數(shù)據(jù)集分解處理的DiscFinder方法。DiscFinder基于初始版本的FOF算法，通過Hadoop和MapReduce實現(xiàn)數(shù)據(jù)并行。文獻(xiàn)[12]介紹了西北大學(xué)的Liu等人為HOP算法設(shè)計的并行方案。他們通過對K-D樹進(jìn)行域分解，利用MPI（message passing interface）實現(xiàn)了多節(jié)點多核的并行加速。HGGF算法與FOF算法和HOP算法相比計算精度更高，與實際觀測數(shù)據(jù)更接近，在星系成團研究中的地位已得到廣泛認(rèn)可[13-14]，但計算過程更復(fù)雜，數(shù)據(jù)依賴更多，尚未有文獻(xiàn)對其模擬運算性能進(jìn)行研究。HGGF算法利用OpenMP共享內(nèi)存實現(xiàn)的多任務(wù)并行效率較低，并且算法中數(shù)據(jù)互相依賴造成外層循環(huán)難以拆分，本文并未如文獻(xiàn)[12]采用多節(jié)點MPI并行，而是在深入分析算法結(jié)構(gòu)后，針對訪存問題和共享變量等問題給出多種優(yōu)化方法，在單節(jié)點多核環(huán)境下獲得很好的效果。

由于訪存延遲常常嚴(yán)重影響程序性能，至今已有多種方法通過改善數(shù)據(jù)局部性優(yōu)化訪存性能，如緩存分塊[15]、單指令多數(shù)據(jù)（single instruction multiple data，SIMD）[16]、預(yù)取技術(shù)等。本文所述方法既包括如數(shù)據(jù)對齊、循環(huán)分解等在文獻(xiàn)[9,17]中記載的方法，也有針對算法特定結(jié)構(gòu)提出的數(shù)據(jù)預(yù)排序方法。主數(shù)組-索引數(shù)組結(jié)構(gòu)在N體模擬問題中應(yīng)用廣泛。文獻(xiàn)[7]介紹的非規(guī)則訪存模型中包含主數(shù)組-索引數(shù)組結(jié)構(gòu)，其中提出的數(shù)據(jù)重排方法基于Inspector-Executor機制[18]，通過在程序內(nèi)插入預(yù)處理函數(shù)動態(tài)調(diào)整數(shù)據(jù)順序。這種方法雖然更靈活，但由于大規(guī)模數(shù)值模擬計算中循環(huán)迭代次數(shù)太多，節(jié)省的時間不足以抵消Inspector函數(shù)的耗時。

本文所述數(shù)據(jù)預(yù)排序方法在計算前的預(yù)處理階段根據(jù)索引數(shù)組特點重排輸入數(shù)據(jù)，與文獻(xiàn)[7]不同，對迭代次數(shù)多的大規(guī)模計算加速更明顯，同時還區(qū)別于文獻(xiàn)[19-20]中記載的直接排序方法，可應(yīng)用于結(jié)構(gòu)更復(fù)雜的訪存模式。

3 HGGF算法及其非規(guī)則訪存模型

3.1HGGF算法

本節(jié)將對算法的核心部分，同時也是程序熱點部分（占總運行時間超過90%）進(jìn)行詳細(xì)介紹。

3.1.1 連接體積

文獻(xiàn)[4]首次提出了連接體積的概念。每個星系在篩選時，本應(yīng)遍歷空間中所有其他星系，但若兩星系間距超過特定閾值，根據(jù)經(jīng)驗它們不可能屬于同一類，沒有互相計算的必要。圖1為星系分布的二維空間抽象圖，每個點代表一個星系。圖中有效半徑外與圓心位置的星系肯定不同群，因此在后面的篩選過程中，只需計算球體（擴展到正方體）內(nèi)的樣本。

初始化連接體積的偽代碼如下：

前面說的有效半徑不是實際的空間距離，而是將三維空間轉(zhuǎn)換成以L為邊長的正方體（簡稱L空間）后按比例換算出的數(shù)值。式（1）中，x、y、z分別代表利用輸入數(shù)據(jù)中星系的極坐標(biāo)求得的直角坐標(biāo)；x′、y′、z′為直角坐標(biāo)經(jīng)過放縮平移到L空間（0＜x′,y′,z′＜L）內(nèi)的值，L′為經(jīng)驗常數(shù)。式（2）中，x″、y″、z″為x′、y′、z′的取整，函數(shù)Func的值為Volume數(shù)組的下標(biāo)（后面統(tǒng)稱為轉(zhuǎn)換坐標(biāo)），并且可以保證僅當(dāng)x″、y″、z″都相等時，F(xiàn)unc的值才會相等。由于x″、y″、z″為原浮點數(shù)坐標(biāo)的取整，誤差可能導(dǎo)致不同星系轉(zhuǎn)換坐標(biāo)相同，此時，Linklist數(shù)組會以類似鏈表的方式逐一記錄下轉(zhuǎn)換坐標(biāo)相等的星系編號，若未與其他星系相同，由Volume數(shù)組記錄星系編號。

Fig.1 Valid radius ofLspace圖1 L空間內(nèi)的有效半徑

經(jīng)過以上步驟，Volume和Linklist兩個數(shù)組記錄了所有的星系編號。進(jìn)入篩選過程前，需要先確定有效半徑內(nèi)包含了哪些星系。本文將以當(dāng)前星系為中心（多次迭代后為星系群的中心），二倍有效半徑為邊長（如圖1虛線框）的正方體稱為連接體積，則連接體積必然包含所有有效半徑內(nèi)的星系。建立索引數(shù)組時，按遞增順序篩選連接體積內(nèi)的所有整數(shù)點，若該點的Func函數(shù)值剛好等于某一星系的轉(zhuǎn)換坐標(biāo)，則將Volume中的星系編號加入索引數(shù)組，偽代碼如下：

掃描結(jié)束后，連接體積內(nèi)所有星系的編號都存入了索引數(shù)組。

3.1.2 遍歷篩選星系

建立索引數(shù)組后進(jìn)入篩選星系的過程，篩選的對象是連接體積內(nèi)的星系，即索引數(shù)組和對應(yīng)Linklist中存儲的星系編號。判斷星系是否同類的條件為式（3），Dist1為紅移距離，Dist2為空間距離，p1、p2、p3、p4和r、x、y、z分別代表被篩選星系和中心星系的坐標(biāo)。每次篩選過程需要兩次判斷，若紅移距離不符合要求，則跳轉(zhuǎn)到下一個星系，若滿足條件，再判斷空間距離，如果兩個條件同時滿足，則記錄下該星系編號。此外，在每次篩選結(jié)束后會查看對應(yīng)Linklist數(shù)組中是否存有其他星系。

3.2 非規(guī)則訪存模型

HGGF算法的非規(guī)則訪存模型為二重循環(huán)，外層循環(huán)為當(dāng)前分類數(shù)，隨著星系成團的過程，星團數(shù)目會不斷減少，兩個內(nèi)層循環(huán)分別對應(yīng)上一小節(jié)建立索引數(shù)組和遍歷篩選星系的過程。第一個內(nèi)層循環(huán)建立的索引數(shù)組用作第二個內(nèi)層循環(huán)中主數(shù)組的下標(biāo)，每個外層循環(huán)索引數(shù)組根據(jù)當(dāng)前星系編號對應(yīng)更新，包含新的要篩選的星系。訪存模型偽代碼如下：

因為編譯時主數(shù)組下標(biāo)即索引數(shù)組的值不能立即獲取，所以訪存的順序完全由索引數(shù)組決定，且連接體積內(nèi)的星系離散分布在所有樣本星系中。以最小輸入數(shù)據(jù)為例，樣本數(shù)為639 359個星系，則索引數(shù)組index的值不規(guī)則分布于區(qū)間1～639 359。主數(shù)組讀取數(shù)據(jù)時，會按一個緩存行的大?。?4 B）讀入，即16個相鄰的數(shù)組元素，但受索引數(shù)組的影響，訪問地址不連續(xù)，需要的數(shù)組元素很難位于同一數(shù)據(jù)塊中，這會嚴(yán)重影響緩存性能。因此，可將這種訪存模式歸為非規(guī)則訪存。

4 優(yōu)化方法

優(yōu)化之前，使用VTune[21]測試出熱點函數(shù)每行代碼的具體耗時以及硬件性能計數(shù)器的數(shù)據(jù)。結(jié)果顯示篩選星系過程耗時占熱點函數(shù)的95%以上，因此優(yōu)化工作主要針對這部分展開。

對于非規(guī)則訪存問題，軟件預(yù)取是優(yōu)化局部性最常見的方法。由于沒有明顯性能提升，只簡要介紹軟件預(yù)取方法：通過插入具體預(yù)取指令，編譯器可提前將數(shù)據(jù)讀入緩存，并利用流水線上其他任務(wù)覆蓋這部分開銷。本文嘗試了兩種：

（1）使用#pragma prefetch指示編譯器進(jìn)行預(yù)??；

（2）使用intrinsic代碼_mm_prefetch將數(shù)組元素預(yù)取到指定緩存。

從profile結(jié)果看，數(shù)據(jù)被成功預(yù)取，但開銷不僅沒有被覆蓋，還和優(yōu)化前一樣。分析原因是需要預(yù)取的計算任務(wù)簡單，流水線不能覆蓋時間更長的訪存開銷[22]。

4.1 非規(guī)則訪存優(yōu)化

4.1.1 數(shù)據(jù)預(yù)排序

通過分析連接體積算法的原理可以發(fā)現(xiàn)，循環(huán)內(nèi)主數(shù)組元素的訪問順序由索引數(shù)組決定，而索引數(shù)組作為Volume數(shù)組的子集，存儲了當(dāng)前星系有效半徑內(nèi)的所有星系編號（除轉(zhuǎn)換坐標(biāo)相等的星系由Linklist數(shù)組存儲）。因為遍歷連接體積時是按轉(zhuǎn)換坐標(biāo)從小到大的順序，所以索引數(shù)組中星系編號的存儲位置同樣應(yīng)該按轉(zhuǎn)換坐標(biāo)從小到大排列。設(shè)想一下，如果能讓輸入數(shù)據(jù)在程序開始前按同樣順序排好，則通過遍歷連接體積內(nèi)整數(shù)點得到的索引數(shù)組，其中的星系編號很可能已經(jīng)是連續(xù)的，即1～639 359中的任意一個連續(xù)區(qū)間，但不包含轉(zhuǎn)換坐標(biāo)相等的星系。

原輸入數(shù)據(jù)中星系的性質(zhì)與轉(zhuǎn)換坐標(biāo)并沒有直接關(guān)系，也沒有與轉(zhuǎn)換坐標(biāo)順序一致或相似的一列數(shù)據(jù)，因此只有通過對輸入數(shù)據(jù)進(jìn)行預(yù)處理，才能達(dá)到使輸入數(shù)據(jù)按轉(zhuǎn)換坐標(biāo)有序的目的，這個預(yù)處理過程稱之為數(shù)據(jù)預(yù)排序。從優(yōu)化的思路來看，數(shù)據(jù)預(yù)排序與Inspector-Executor有一定相似之處，兩種方法都是通過改變數(shù)據(jù)的分布來優(yōu)化訪存性能，但前者是在程序運行前通過預(yù)處理輸入數(shù)據(jù)改變數(shù)據(jù)局部性，而Inspector-Executor是通過在程序中插入數(shù)據(jù)重排函數(shù)，動態(tài)改變數(shù)據(jù)的分布。

下面介紹數(shù)據(jù)預(yù)排序的具體操作方法。如圖2，首先利用與原程序相同的算法求出所有星系的轉(zhuǎn)換坐標(biāo)，然后將得到的值按原輸入數(shù)據(jù)順序添加到文件中（原輸入數(shù)據(jù)中沒有轉(zhuǎn)換坐標(biāo)那一列，相當(dāng)于給每個星系增加一個新的性質(zhì)），再按轉(zhuǎn)換坐標(biāo)從小到大對輸入文件按行排序，并作為新的輸入文件，整個過程并不需要修改源程序的代碼。

Fig.2 Procedure of data pre-sorting圖2 數(shù)據(jù)預(yù)排序步驟

圖3是數(shù)據(jù)預(yù)排序前后索引數(shù)組分配的抽象圖，可看出預(yù)排序前，各星系編號被隨機分配到索引數(shù)組中，這樣讀取數(shù)據(jù)時，由于星系編號差別很大，導(dǎo)致每次讀入的一行數(shù)據(jù)中可能只有一個有效數(shù)據(jù)。預(yù)排序后，由于星系編號已經(jīng)按轉(zhuǎn)換坐標(biāo)排序，會按順序被分配到索引數(shù)組中，使得每次讀入緩存的數(shù)據(jù)塊中的全部數(shù)組元素都是本次循環(huán)或下個循環(huán)要被使用的數(shù)據(jù)。圖4是數(shù)據(jù)預(yù)排序前后索引數(shù)組（陰影部分）及Linklist數(shù)組（突出部分）存儲的星系編號，因為遍歷的順序類似于深度優(yōu)先搜索，所以預(yù)排序后數(shù)組內(nèi)星系編號可視為完全連續(xù)。雖然主數(shù)組-索引數(shù)組的結(jié)構(gòu)未變，但是數(shù)據(jù)局部性已經(jīng)徹底得到改善。

Fig.3 Distribution of index array圖3 索引數(shù)組分配

Fig.4 Change of index and Linklist array圖4 索引數(shù)組及Linklist數(shù)組變化

相比訪存延遲的巨大額外開銷，數(shù)據(jù)預(yù)排序本身的開銷可以被抵消掉，而且僅通過一次預(yù)排序，可以為之后多次模擬實驗節(jié)省時間，顯然是可以接受的。

4.1.2 數(shù)據(jù)對齊

篩選星系過程的主數(shù)組是二維數(shù)組，記錄了星系的全部16個性質(zhì)，并且聲明在內(nèi)存的相鄰位置，但是計算過程中只訪問到16個元素中的前4個（三維坐標(biāo)、紅移距離），這會影響數(shù)據(jù)讀取的效率。如圖5（a）是數(shù)據(jù)未對齊時最好的情況，即一次將4個數(shù)組元素全部讀入，但是浪費了后面的12個單元。圖5（b）則是最差的情況，4個元素由于前端未對齊被分開，此時需要訪問兩次內(nèi)存才能取出全部4個元素。圖5（c）是使用數(shù)據(jù)對齊，并且將數(shù)組分開聲明的情況，再配合數(shù)據(jù)預(yù)排序，一次可同時讀入4組有效數(shù)據(jù)，比之前最差的情況效率提高8倍。

Fig.5 Data alignment圖5 數(shù)據(jù)對齊

數(shù)據(jù)對齊方法：首先確保數(shù)組的起始地址對齊，使用可以將普通的內(nèi)存分配和釋放改寫為對齊的intrinsic原語，即將malloc()改為_mm_malloc()，free()改為_mm_free()，由于剛好被16整除，不需要考慮后端對齊。

4.1.3 循環(huán)分解

由于CPU緩存空間有限，循環(huán)中無關(guān)數(shù)據(jù)過多可能使原本應(yīng)保留在緩存中的數(shù)據(jù)被擠出，而換入的數(shù)據(jù)若再被換出，會形成惡性循環(huán)，導(dǎo)致有效數(shù)據(jù)總是不能保留在緩存中。如果將沒有依賴關(guān)系的計算拆開，分成兩個或多個循環(huán)，不僅可以提高各部分對緩存的利用效率，還可以分別展開進(jìn)一步優(yōu)化。

4.2 OpenMP并行效率優(yōu)化

本文使用OpenMP對篩選星系過程實現(xiàn)并行，但在12線程上只加速不到4倍，因此使用了以下方法進(jìn)一步優(yōu)化算法。

4.2.1 負(fù)載均衡

為避免各線程任務(wù)分配不均，可利用OpenMP制導(dǎo)語句選項nowait或schedule實現(xiàn)任務(wù)的動態(tài)分配。加入選項后，編譯器會根據(jù)各線程的執(zhí)行情況動態(tài)分配任務(wù)，先執(zhí)行完的線程會被分配新的任務(wù)。兩種方法對HGGF算法的效果基本相同。

4.2.2 互斥變量私有化

因為程序中的部分操作和變量在迭代過程中只能串行執(zhí)行或由多線程共享，所以需要使用鎖機制避免該部分程序因為多線程并行出現(xiàn)讀寫沖突。OpenMP中實現(xiàn)互斥最常見的方法是使用指導(dǎo)選項#ompcritical或#ompatomic。添加指導(dǎo)選項后，同一時間只能有一個線程對互斥區(qū)進(jìn)行訪問。需要特別指出，atomic比critical雖然在性能上會有小幅提升，但由于atomic選項只能作用于簡單語句，當(dāng)互斥區(qū)內(nèi)指令比較復(fù)雜時，可能造成結(jié)果錯誤。

雖然指導(dǎo)選項使用便捷，但線程間只能互相等待，嚴(yán)重影響程序并行效率。HGGF算法中互斥區(qū)內(nèi)指令的作用是計數(shù)及記錄星系編號，線程執(zhí)行的先后順序不影響最終結(jié)果。利用這個特點，將線程序號與數(shù)組下標(biāo)結(jié)合，為每個線程單獨聲明一組變量，使得原來需要互斥訪問的共享變量可以被各線程私有和同步讀寫，實現(xiàn)了線程間的完全并行。偽代碼如下：

4.3 MIC移植介紹

MIC采用傳統(tǒng)x86架構(gòu)，擁有數(shù)十個精簡核心，提供高度并行的計算能力，支持現(xiàn)有程序的平滑移植。MIC有靈活的編程模式，既可以用作協(xié)處理器，也可以作為獨立的節(jié)點，Native和Offload是兩種比較常見的模式[23]。但是，MIC并不總是能提供好的性能，其CPU核心不支持亂序執(zhí)行，復(fù)雜算法在MIC上的表現(xiàn)一般。與傳統(tǒng)多核CPU相比，MIC的SIMD和緩存對算法更敏感，訪存密集或訪存不規(guī)則的程序受到的影響更加嚴(yán)重。

為了進(jìn)一步驗證文中優(yōu)化方法的有效性，以及探索未來程序向眾核平臺移植的可能性，在MIC的Natvie模式下測試了HGGF應(yīng)用，僅比較并行部分的耗時，即使開啟全部60個核上的60/120/180/240個線程，并行性能仍落后于CPU。

要在MIC上獲得好的性能，往往需要從多方面進(jìn)行優(yōu)化。受到算法限制，緩存分塊、單指令多數(shù)據(jù)等技術(shù)并不適用于HGGF算法。但考慮到同樣是x86架構(gòu)，前面介紹的優(yōu)化方法應(yīng)該也適用于MIC。因此本文將同樣的優(yōu)化方法應(yīng)用在MIC的程序中，取得了與CPU上接近的加速比。

5 實驗環(huán)境及結(jié)果分析

5.1 實驗環(huán)境

5.1.1 數(shù)據(jù)集

原始HGGF算法已經(jīng)在單線程機器上驗證并使用過。其使用的輸入數(shù)據(jù)是基于P3M算法得到的仿真數(shù)據(jù)[24]，輸出數(shù)據(jù)正確性可由2度視場星系紅移巡天計劃（2-degree Field Galaxy Redshift Survey，2dFGRS）和斯隆巡天計劃（Sloan Digital Sky Survey，SDSS）數(shù)據(jù)驗證。在本次并行化的升級算法中，將和原來的單線程結(jié)果相互比較驗證。

5.1.2 軟硬件環(huán)境

實驗的硬件環(huán)境如表1所示，軟件環(huán)境如表2所示。

Table 1 Hardware表1 硬件環(huán)境

Table 2 Software表2 軟件環(huán)境

5.2 結(jié)果分析

5.2.1 非規(guī)則訪存優(yōu)化

為了驗證數(shù)據(jù)預(yù)排序?qū)υL存性能的影響，利用VTune測出了程序單核運行時的各級緩存丟失率。從表3中可以看出，使用該方法后，各級緩存丟失率均有明顯下降。這說明數(shù)據(jù)預(yù)排序已經(jīng)改善了局部性，使得CPU每次都能讀取更多的有效數(shù)據(jù)，與本文的分析基本一致。

Table 3 Cache miss rate during pre-sorting表3 預(yù)排序前后緩存丟失對比 %

確認(rèn)了預(yù)排序方法的作用后，又在多線程上測試了程序的運行時間。如圖6，單獨使用預(yù)排序方法時（紅色線），單線程上的加速效果最明顯，有1.9倍加速，而使用全部12個線程后，僅有1.2倍加速。原因是使用多核并行后，每個核有獨立的一、二級緩存，等效于每次循環(huán)可將更多數(shù)據(jù)讀入緩存，削弱了數(shù)據(jù)預(yù)排序的效果。此外，數(shù)據(jù)對齊和緩存分解方法是已有記載的優(yōu)化方法，不僅針對單一訪存模型，無需在此詳細(xì)分析。使用全部3種非規(guī)則訪存優(yōu)化方法（藍(lán)色線）的運行時間，相比僅使用預(yù)排序方法有1.2倍加速，在不同核數(shù)上加速效果基本相同。因為核數(shù)增多后，數(shù)據(jù)預(yù)排序加速效果減弱，所以設(shè)計了使用更大數(shù)據(jù)集的實驗。如圖7，數(shù)據(jù)規(guī)模從1到16倍，預(yù)排序的加速比從1.2倍增加到1.6倍。原因是隨著數(shù)據(jù)規(guī)模擴大，緩存容量再次不足，數(shù)據(jù)預(yù)排序的效果逐漸回升，證明其弱可擴展性較好。

Fig.6 Strong scalability of memory optimization圖6 訪存優(yōu)化強可擴展性

Fig.7 Weak scalability of data pre-sorting圖7 數(shù)據(jù)預(yù)排序弱可擴展性

Fig.8 Strong scalability of OpenMP optimization圖8 OpenMP并行優(yōu)化強可擴展性

5.2.2 OpenMP并行效率優(yōu)化

如圖8，分別測試了負(fù)載均衡（紅色線）和互斥變量私有化（藍(lán)色線）兩種方法的運行時間。可以看出負(fù)載均衡對性能的提升并不明顯，只有5%左右，說明原算法線程任務(wù)分配已經(jīng)比較均勻。而互斥變量私有化隨著使用核數(shù)的增加加速效果越來越明顯，12線程直接加速2倍，這也說明原算法中互斥區(qū)的存在嚴(yán)重影響了線程的并行。

最后來看一下綜合所有優(yōu)化方法的加速效果及可擴展性。如表4，優(yōu)化后12線程的加速比為11.6倍，相比優(yōu)化前僅有3.9倍的加速有了大幅提高。而優(yōu)化后整個程序（不僅Kernel部分）的加速比也達(dá)到8倍。圖9顯示了程序經(jīng)整體優(yōu)化后弱可擴展性良好。

5.2.3 MIC優(yōu)化結(jié)果

本文分別測試了使用數(shù)據(jù)預(yù)排序方法和全部優(yōu)化方法在MIC上的運行時間。如圖10，僅使用數(shù)據(jù)預(yù)排序時，最多加速1.4倍，但不同線程數(shù)的加速效果相差不大，主要是因為在MIC上每個核最多可開啟4個線程，雖然線程數(shù)不同，但是使用的核數(shù)都是60個，所以并未改變緩存的容量，就不會對數(shù)據(jù)預(yù)排序的效果產(chǎn)生影響。使用全部優(yōu)化方法后，最多在240個線程上加速2.2倍，主要原因是線程越多，互斥變量私有化方法的效果越明顯。

Table 4 Speedup of all optimizations表4 所有優(yōu)化方法加速情況

Fig.9 Weak scalability of final optimization圖9 最終優(yōu)化弱可擴展性

Fig.10 Result on MIC圖10 MIC測試結(jié)果

6 總結(jié)與下一步工作

本文通過對天文應(yīng)用Halo-based Galaxy Group Finder進(jìn)行分析與優(yōu)化，歸納出算法核心部分的非規(guī)則訪存模型，并對該類型問題使用數(shù)據(jù)預(yù)排序方法優(yōu)化訪存性能，取得了不錯的優(yōu)化結(jié)果。還分析了數(shù)據(jù)預(yù)排序方法對該非規(guī)則訪存模型產(chǎn)生影響的原因，為類似的訪存問題提供了新的優(yōu)化思路和方法。在數(shù)據(jù)預(yù)排序的基礎(chǔ)上，使用了幾種傳統(tǒng)優(yōu)化方法：（1）使用數(shù)據(jù)對齊、循環(huán)分解優(yōu)化訪存性能；（2）使用負(fù)載均衡、互斥變量私有化提高算法的并行效率。綜合上述方法，在單節(jié)點12個線程上對HGGF應(yīng)用加速11.6倍，并且隨著數(shù)據(jù)規(guī)模的擴大，加速效果更加明顯。將HGGF應(yīng)用以本地模式移植到MIC卡上，并測試了性能，驗證了本文優(yōu)化方法的有效性。綜上，本文的優(yōu)化成果對于HGGF應(yīng)用涉及的天文學(xué)問題的研究將有很大幫助，也為相似訪存類型應(yīng)用提供了更多優(yōu)化途徑。

在下一步工作中，將繼續(xù)通過優(yōu)化算法改善其數(shù)據(jù)局部性和可擴展性，并利用多節(jié)點（MPI）或異構(gòu)平臺（CPU+GPU，CPU+MIC）進(jìn)行加速。

致謝感謝上海交通大學(xué)高性能計算中心提供的實驗環(huán)境。感謝上海交通大學(xué)物理與天文系楊小虎教授在程序優(yōu)化過程中的悉心指導(dǎo)，以及對本文的修改意見。

[1]Yang Xiaohu,Mo H J,van den Bosch F C,et al.A halobased galaxy group finder:calibration and application to the 2dFGRS[J].Monthly Notices of the Royal Astronomical Society,2005,356(4):1293-1307.

[2]Huchra J P,Geller M J.Groups of galaxies:I-nearby groups [J].TheAstrophysical Journal,1982,257:423-437.

[3]Yang Xiaohu,Mo H J,van den Bosch F C,et al.Galaxy groups in the SDSS DR4.I.the catalog and basic properties [J].TheAstrophysical Journal,2007,671(1):153.

[4]Eke V R,Baugh C M,Cole S,et al.Galaxy groups in the2dFGRS:the group-finding algorithm and the 2PIGG catalogue[J].Monthly Notices of the Royal Astronomical Society,2004,348(3):866-878.

[5]Shevtov M,Soupikov A,Kapustin A.Highly parallel fast KD-tree construction for interactive ray tracing of dynamic scenes[J].Computer Graphics Forum,2007,26(3):395-404.

[6]Badawy A H,Aggarwal A,Yeung D,et al.The efficacy of software prefetching and locality optimizations on future memory systems[J].Journal of Instruction-Level Parallelism,2004,6(7).

[7]Lin Yuan,Padua D.Analysis of irregular single-indexed array accesses and its applications in compiler optimizations [C]//LNCS 1781:Proceedings of the 2000 Joint European Conferences on Theory and Practice of Software,Berlin, Mar 25-Apr 2,2000.Berlin,Heidelberg:Springer,2000: 202-218.

[8]Hennessy J L,Patterson D A.Computer architecture:a quantitative approach[M].San Francisco,USA:Morgan Kaufmann Publishers Inc,2011.

[9]Lin Xinhua,Li Shuo,Zhao Jiaming,et al.Node-level memory access optimization on Intel Knights Corner[J].Computer Science,2015,42(11):37-42.

[10]Gornish E H,Granston E D,Veidenbaum A V.Compilerdirected data prefetching in multiprocessors with memory hierarchies[C]//Proceedings of the 25th International Conference on Supercomputing Anniversary Volume,Tucson, USA,May 31-Jun 4,2011.New York:ACM,2014:128-142.

[11]Fu Bin,Ren Kai,Lopez J,et al.DiscFinder:a data-intensive scalable cluster finder for astrophysics[C]//Proceedings of the 19th International Conference on High Performance Distributed Computing,Chicago,USA,Jun 20-25,2010. New York:ACM,2010:348-351.

[12]Liu Ying,Liao Weikeng,Choudhary A.Design and evaluation of a parallel HOP clustering algorithm for cosmological simulation[C]//Proceedings of the 17th International Parallel and Distributed Processing Symposium,Nice,France,Apr 22-26,2003.Piscataway,USA:IEEE,2003:82.

[13]Faltenbacher A,Li Cheng,Mao Shude,et al.Three different types of galaxy alignment within dark matter halos[J]. TheAstrophysical Journal Letters,2007,662(2):L71-L74.

[14]Yang Xiaohu,Mo H J,van den Bosch F C.Galaxy groups in the SDSS DR4.II.halo occupation statistics[J].The Astrophysical Journal,2008,676(1):248-261.

[15]Nishtala R,Vuduc R W,Demmel J W,et al.When cache blocking of sparse matrix vector multiply works and why [J].Applicable Algebra in Engineering,Communication and Computing,2007,18(3):297-311.

[16]Knobe K,Lukas J D,Steele G L.Data optimization:allocation of arrays to reduce communication on SIMD machines [J].Journal of Parallel and Distributed Computing,1990,8 (2):102-118.

[17]Hager G,Wellein G.Introduction to high performance computing for scientists and engineers[M].Boca Raton,USA: CRC Press,Inc,2010.

[18]Arenaz M,Tourino J,Doallo R.An inspector-executor algorithm for irregular assignment parallelization[C]//LNCS 3358:Proceedings of the 2nd International Symposium on Parallel and Distributed Processing and Applications,Hong Kong,China,Dec 13-15,2004.Berlin,Heidelberg:Springer, 2005:4-15.

[19]Oh G H,Kim J M,Kang W H,et al.Reducing cache misses in hash join probing phase by pre-sorting strategy[C]//Proceedings of the 2012 SIGMOD International Conference on Management of Data,Scottsdale,USA,May 20-24,2012. New York:ACM,2012:864-864.

[20]Abouzied A,Abadi D J,Silberschatz A.Invisible loading: access-driven data transfer from raw files into database systems[C]//Proceedings of the 16th International Conference on Extending Database Technology,Genoa,Italy,Mar 18-22,2013.New York:ACM,2013:1-10.

[21]Wang P.How to use Vtune amplifier XE 2013 on Intel Xeon Phi coprocessor[R/OL].Intel Corporation(2012)[2015-10-19].https://software.intel.com/zh-cn/blogs/2012/11/05/vtuneamplifier-xe-2013-intel-xeon-phi-coprocessor?language=it. [22]Lee J,Kim H,Vuduc R W.When prefetching works,when it doesn’t,and why[J].ACM Transactions on Architecture and Code Optimization,2012,9(1):2.

[23]Jeffers J,Reinders J.Intel Xeon Phi coprocessor high performance programming[M].Beijing:Posts and Telecom Press,2014.

[24]Jing Y P,Suto Y.Triaxial modeling of halo density profiles with high-resolutionN-body simulations[J].The Astrophysical Journal,2002,574(2):538-553.

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

[9]林新華,李碩,趙嘉明,等.Intel Knights Corner的結(jié)點級內(nèi)存訪問優(yōu)化[J].計算機科學(xué),2015,42(11):37-42.

[23]Jeffers J,Reinders J.Intel Xeon Phi協(xié)處理器高性能編程指南[M].北京:人民郵電出版社,2014.

HAO He was born in 1989.He is an M.S.candidate at Shanghai Jiao Tong University.His research interest is high performance computing.

郝赫（1989—），男，遼寧錦州人，上海交通大學(xué)碩士研究生，主要研究領(lǐng)域為高性能計算。

SI Yumeng was born in 1992.He is an M.S.candidate at Shanghai Jiao Tong University.His research interest is high performance computing.

司雨蒙（1992—），男，山西晉城人，上海交通大學(xué)碩士研究生，主要研究領(lǐng)域為高性能計算。

WEI Jianwen was born in 1986.He is an assistant engineer at the Center for High Performance Computing of Shanghai Jiao Tong University.His research interest is high performance computing.

韋建文（1986—），男，廣西百色人，上海交通大學(xué)高性能計算中心助理工程師，主要研究領(lǐng)域為高性能計算。

WEN Minhua was born in 1988.He is an assistant engineer at the Center for High Performance Computing of Shanghai Jiao Tong University.His research interest is high performance computing.

文敏華（1988—），男，江西會昌人，上海交通大學(xué)高性能計算中心助理工程師，主要研究領(lǐng)域為高性能計算。

LIN Xinhua was born in 1979.He is the vice director at the Center for High Performance Computing of Shanghai Jiao Tong University.His research interests include computer architecture and code optimization.

林新華（1979—），男，浙江紹興人，上海交通大學(xué)高性能計算中心副主任，主要研究領(lǐng)域為體系結(jié)構(gòu)，代碼優(yōu)化。

Optimizing Irregular MemoryAccess inAstrophysical Clustering Studies*

HAO He1,SI Yumeng1,WEI Jianwen1,WEN Minhua1,LIN Xinhua1,2+
1.Center for High Performance Computing,Shanghai Jiao Tong University,Shanghai 200240,China
2.NVIDIATechnology CenterAsia Pacific,Singapore 999002
+Corresponding author:E-mail:james@sjtu.edu.cn

Halo-based galaxy group finder(HGGF)tries to find galaxies in the same dark matter halo which is not directly visible.It plays a very important role in the research of large-scale structure of the universe.However,because of the growth of data scale,it’s extremely necessary to increase the running speed by optimizing the group finder coding algorithm.After a thorough investigation on the original HGGF code,it is found that the kernel part of the algorithm is seriously affected by the irregular memory access.This paper proposes a specific data pre-sorting approach and analyzes how it affects the process of memory access according to the structure of the algorithm and the irregular memory access pattern.Moreover,this paper uses data alignment and loop fission to optimize the memory access as well as improving the efficiency of OpenMP with load balance and mutex privatization.Eventually the HGGF application gets 11.6 times speedup on 12 threads,and gets better weak scalability.The following is the original contributions:(1)Analyze the irregular memory access of the HGGF application;(2)Propose and analyze the data pre-sorting;(3)Improve the parallel performance of HGGF application with another four approaches including data alignment,loop fission,load balance and mutex privatization.

A

：TP391

10.3778/j.issn.1673-9418.1512078

*The National High Technology Research and Development Program of China under Grant No.2014AA01A302(國家高技術(shù)研究發(fā)展計劃(863計劃));the Project of Japan Society for the Promotion of Science RONPAKU Fellowship(日本學(xué)術(shù)振興會項目).

Received 2015-11,Accepted 2016-01.

CNKI網(wǎng)絡(luò)優(yōu)先出版:2016-01-07,http://www.cnki.net/kcms/detail/11.5602.TP.20160107.1540.008.html

Key words:astrophysical clustering studies;optimizing irregular memory access;data pre-sorting;parallel computing