畢玉江,周 超,吳郁非,黎睿翔,劉朝峰,3,陳建海,徐 順

1(中國科學院 高能物理研究所,北京 100049)

2(中國科學院計算機網(wǎng)絡信息中心,北京 100190)

3(中國科學院大學,北京 100049)

4(浙江大學,杭州 310058)

在微觀粒子世界中起主要作用的是強相互作用、電磁相互作用和弱相互作用,描述強相互作用的量子色動力學(QCD)和描述電磁、弱相互作用的電弱統(tǒng)一理論(EW)統(tǒng)稱為粒子物理中的標準模型.就強相互作用而言,其低能特性如夸克禁閉問題,仍是需要解決的世紀難題.實際上,對于夸克禁閉以及與之密切相關(guān)的強子結(jié)構(gòu)和強子性質(zhì)等低能區(qū)強相互作用現(xiàn)象需要有效的非微擾方法.格點QCD[1]是目前公認有效的從第一原理出發(fā)的求解QCD的非微擾方法,除標準模型基本參數(shù)之外沒有任何額外模型參數(shù),因而其計算結(jié)果被認為是對強相互作用現(xiàn)象的可靠描述.該格點QCD理論對不少重大物理研究有著不可替代的作用,主要體現(xiàn)在用高能物理實驗結(jié)果對現(xiàn)有理論的驗證、用理論分析結(jié)果指導實驗設計、以及進一步地尋找標準模型以外的新物理等方面.

格點QCD的研究手段是通過計算機進行大規(guī)模的矩陣運算和蒙特卡羅數(shù)值模擬,其算法特點具有內(nèi)在的高并行度和高可擴展性.然而,隨著計算機硬件體系結(jié)構(gòu)的日益復雜,研發(fā)計算高效的格點QCD應用程序也變得愈加挑戰(zhàn)性.大規(guī)模格點QCD計算過程對高性能計算資源需求巨大,其計算過程體現(xiàn)出計算密集型特征,而其數(shù)據(jù)分析體現(xiàn)出數(shù)據(jù)密集型特征.計算精度和數(shù)據(jù)分析效率直接反映了格點QCD計算方法的應用水平.格點QCD已經(jīng)逐漸成為與高能物理實驗研究和理論解析研究并列的高能物理第三大分支,研究范圍已經(jīng)拓展到強相互作用現(xiàn)象的各個方面.美國、日本、歐洲等國家無一例外將格點QCD應用作為高性能計算的重要應用.美國更是將其列為未來E級計算重要的應用之一,凸顯了格點QCD數(shù)值模擬內(nèi)在的科學計算應用的價值[2].近年來,格點QCD的數(shù)值模擬方法發(fā)展出了大量理論和算法,形成了美國USQCD[3]和MILC[4]、英國Grid[5]和日本Bridge++軟件[6]等諸多優(yōu)秀的格點QCD模擬軟件.

從格點QCD數(shù)值計算的特點來分析,其主要過程為矩陣向量間的運算,并行化時涉及格點體系的區(qū)域分解并行劃分,邊界格子進程間通信還有組態(tài)存儲分析等任務[7].如何實現(xiàn)高效的大規(guī)模格點QCD數(shù)值計算一直以來都是一個現(xiàn)實問題.但我們研究發(fā)現(xiàn),格點QCD抽象分層軟件架構(gòu)的設計方法得到了重要應用,軟件框架的這種分層抽象的設計方法極大地提高了代碼的開發(fā)效率和運行效率.比如美國主導的USQCD軟件框架分為4層:最上層應用層包括Chroma、CPS和MILC 等,下面3、4層為支撐庫層,分別實現(xiàn)消息通訊、數(shù)據(jù)結(jié)構(gòu)定義和基本算法實現(xiàn),而其中的QDP++支撐庫是用來實現(xiàn)數(shù)據(jù)并行應用接口的,提供邏輯上的數(shù)據(jù)并行操作,隱藏底層數(shù)據(jù)依賴體系結(jié)構(gòu)的操作,提供上層統(tǒng)一的數(shù)據(jù)通信函數(shù)接口.同樣在這方面,英國Peter Boyle 等人2015年研發(fā)出了Grid軟件[5],它抽象出了一個數(shù)據(jù)并行計算模式,針對格點QCD的熱點計算抽象了統(tǒng)一的數(shù)據(jù)并行接口,支持多種平臺的SIMD矢量化計算,如AVX2、AVX512、NEON等形式,同時也支持組合OpenMP和MPI 等多種并行計算編程模式.在實現(xiàn)線程化、SIMD操作可以實現(xiàn)比QDP++更有效.

Grid軟件包近年來受到越來越多的應用,利用它開展的物理研究包括K介子衰變中直接電荷共軛-宇稱(CP)聯(lián)合破壞參數(shù)的格點QCD計算[8],這能幫助物理學家理解宇宙中觀測到的正反物質(zhì)不對稱性.文獻[8]利用Grid軟件改進了他們以前的物理測量代碼,性能得到顯著提升.另外核子結(jié)構(gòu)的研究[9]以及超出標準模型新物理的研究[10]等也使用了Grid軟件.文獻[11]中作者們初步嘗試將Grid軟件移植到GPU上,他們比較了幾種策略:OpenACC、OpenMP 4.5、Just-In-Time編譯和CUDA 并討論了各自的優(yōu)缺點.

從現(xiàn)有工作來看,格點QCD計算軟件是以分層軟件架構(gòu)為主,以數(shù)據(jù)并行計算模式為核心,如何提升格點QCD中數(shù)據(jù)并行計算效率和傳輸效率成為關(guān)鍵點,這是現(xiàn)有USQCD和Grid 等軟件包重要的研發(fā)目標.本文基于格點QCD計算Grid軟件包,開展格點QCD數(shù)據(jù)并行計算的理論分析,研究數(shù)據(jù)并行計算模式在實現(xiàn)上的關(guān)鍵問題,分析典型應用軟件Grid和QDP++的軟件架構(gòu)實現(xiàn),最后結(jié)合具體的計算平臺,測試相關(guān)的應用軟件性能,針對實驗數(shù)據(jù)開展分析,驗證相應的理論分析的有效性.

1 格點QCD模擬數(shù)據(jù)并行計算分析

1.1 格點QCD計算算法原理

格點QCD數(shù)值計算的基本思想為:1)將四維連續(xù)時空離散化為一個四維超立方時空格子Lx×Ly×Lz×Lt,這樣體系的自由度變?yōu)橛邢蘅蓴?shù);2)在該時空格子上定義基本自由度夸克場和膠子場,并在此基礎上定義QCD的格點作用量形式;3)采用路徑積分量子化的方式;4)物理觀測量的測量類似統(tǒng)計物理中的系綜平均值.于是問題就轉(zhuǎn)化為一個包含近鄰相互作用(因果律要求物理體系的相互作用必須是局域的)、多自由度體系的統(tǒng)計物理問題,可以采用蒙特卡羅數(shù)值模擬方法對物理觀測量進行計算.格點QCD數(shù)值模擬計算的簡單流程如下:

1)按照Boltzmann分布因子進行蒙特卡羅隨機抽樣,產(chǎn)生QCD物質(zhì)狀態(tài)系綜(規(guī)范場組態(tài)),這個過程中會涉及大量的線性系統(tǒng)求解 Dx=b 來用于組態(tài)更新,其中涉及到維數(shù)達到Lx×Ly×Lz×Lt×4×3 (一般超過千萬數(shù)量級)的稀疏矩陣(以典型的威爾遜格點作用量為例):

這里n和m的取值為整個格子Lx×Ly×Lz×Lt,α和β為夸克自旋指標分別有4個取值,下標a和b為夸克的顏色自由度分別有3個取值.

2)物理量測量,這個過程涉及在第1)步生成的每一個組態(tài)上做大量的線性系統(tǒng)求解 D x=b.與1)的不同在于向量b的設置不一樣以及解向量x 這里是用于物理量的計算.

3)對物理量做組態(tài)平均,通過數(shù)據(jù)擬合與分析獲得物理結(jié)果.

總體而言,格點QCD的主要計算熱點就是求解大規(guī)模線性方程組D x=b.人們一般采用CG、BiCGstab、MinRes、GCR 等迭代算法,并且配合Even-odd、Multi-grid 等適當?shù)念A處理算法,在{b,Db,D2b,D3b,···} 張開的Krylov 子空間里迭代求解所需的大規(guī)模線性方程組 D x=b.在求解過程中,這些算法會涉及大規(guī)模數(shù)據(jù)處理操作,因而“矩陣D 乘以向量”操作是計算熱點的熱點.另外,在超大規(guī)模計算中,“計算兩個向量內(nèi)積”由于涉及到全局歸約,也會產(chǎn)生額外開銷.

1.2 格點QCD計算的數(shù)據(jù)并行計算特征

格點QCD數(shù)值模擬時,四維超立方格子上的點被分配到多個處理核做并行計算,三維中的示意圖如圖1所示.

圖1 格點QCD計算中的空間分解并行劃分方法

由于QCD 具有近鄰相互作用特征,每個處理器核心在涉及小立方邊界處就會向其他相關(guān)核心發(fā)起通信以獲得相鄰立方邊界的數(shù)據(jù).因此切分格子的方法會對計算效率產(chǎn)生明顯的影響,需考慮計算與通信的平衡.

1.3 數(shù)據(jù)并行計算的典型設計

Grid 最重要也是其突出的特點就體現(xiàn)在數(shù)據(jù)處理的矢量化設計上.前面提到過,在格點QCD計算中,為了能夠適應多節(jié)點并行計算,通常是將格子在時空維度進行劃分,然后在單個節(jié)點上進行各種矩陣的計算,這是很自然的想法.在當時處理器沒有向量化或向量化能力比較弱的情況下,這種處理方法沒有什么問題的.早期的格點QCD 程序如QDP++在考慮處理器矢

圖2 QDP++矩陣向量乘法示意圖

量化時,一般考慮的就是在單個矩陣計算時的矢量化(如圖2,其中等號左邊列向量中實點的值等于等號右邊矩陣中實點所在行與列向量的內(nèi)積(圖引自文獻[5]),其并行方法是在單個矩陣計算內(nèi),將多個元素的內(nèi)積及規(guī)約進行并行),這需要MPI進程之間進行很多的歸約,形成通訊上的瓶頸,對性能造成影響.但是數(shù)據(jù)的劃分方式除了時空劃分這種方式,還可以通過每個格點上的自旋和洛侖茲指標進行劃分,Grid在設計時,就部分采用了這種劃分方式,采用了與矩陣內(nèi)矢量化不同的方式,如圖3,其考慮兩個或4個矩陣一起對每個元素進行計算,在SIMD線上,沒有MPI歸約,減少了通訊開銷.圖3中,SIMD可實現(xiàn)多個矩陣向量乘法并行操作,多個向量對等位上的數(shù)存儲在一個SIMD向量去計算,避免了SIMD位長限制和規(guī)約操作(圖3引自文獻[5]).

圖3 Grid矩陣向量乘法的SIMD 實現(xiàn)示意圖

這種方式很適合格點 QCD計算,考慮到AVX512有8個double 寬度,由于要計算的矩陣有4 維旋量指標,自然就適合4個矩陣一起進行計算,而且矢量化也非常簡便.Grid軟件中采用的這種SIMD計算方式充分發(fā)揮了數(shù)據(jù)并行計算的潛在能力.此外,這種SIMD矢量化方式對各種數(shù)值運算都進行了矢量化抽象,只需要完成相應的數(shù)據(jù)矢量化代碼,Grid就能夠很容易地擴展到新的硬件架構(gòu)中.我們在開發(fā)格點QCD代碼時,也可以借鑒Grid優(yōu)秀的設計思想.

2 實驗結(jié)果與分析

2.1 AVX512與AVX2 性能對比

選擇Grid軟件Wilson 實例作為測試對象,我們比較了在Intel Xeon Gold 5120 機器在采用AVX2和AVX512向量化版本的性能差異,性能結(jié)果數(shù)據(jù)如表1和圖4.

表1 AVX512與AVX2 性能對比

圖4 Intel Xeon Gold 5120上AVX512與AVX2性能加速對比圖

從圖4中,我們可以看出浮點計算性能在AVX512版本下優(yōu)于AVX2,但隨著線程數(shù)(Threads)增加,越來越接近總核數(shù)(28),會發(fā)現(xiàn)其性能對比差異開始下降.

2.2 OpenMP 共享內(nèi)存測試

接下來分析線程綁定對數(shù)據(jù)并行計算的影響.我們選擇在2路2核的Intel Xeon Silver 4114機器上做線程綁定測試,采用Intel編譯器來編譯Grid軟件,通過設置環(huán)境變量KMP_AFFINITY 實現(xiàn)線程綁定,主要分為兩種模式,即compcat和scatter.compact模式下,線程將優(yōu)先集中在一個CPU上進行計算;scatter模式下,線程將均勻分布在多個CPU上進行計算,訪存變大.

具體測試方案:分別在scatter模式下和compact模式下進行線程綁定的測試,分別測試線程數(shù)為1,2,3,4,8,10,12,16,20,30,32,36和40 下的計算性能,得到結(jié)果如表2,同時圖示化對比結(jié)果如圖5所示.

表2 線程綁定scatter與compact模式下性能對比

圖5 不同線程綁定方式性能對比圖.

由于程序存在大量的cache 訪問,在線程數(shù)分別為2,3,4,8,10,12,16和20的情況下,線程綁定導致共享相同cache,造成沖突,競爭使用cache,造成cache 頻繁讀寫,降低訪存速度,相比分散線程,性能要差.如果應用程序需要大量的訪存操作,需要大量的cache,將OpenMP線程綁定核時,盡量將線程分布在每個CPU,可以很好的提升性能,在本次測試中,性能較compact 下最高提升1.7倍.如果程序的每個線程需要共享大量數(shù)據(jù),則將所有線程綁定在一個CPU時,性能更好.

2.3 ARM架構(gòu)向量化性能對比

如今越來越多的服務器開始采用ARM架構(gòu)的處理器,比如華為的鯤鵬處理器,亞馬遜的Graviton處理器,飛騰的FT系列處理器等等,“天河三號”原型系統(tǒng)便使用了64位ARM架構(gòu)的64核心處理器 FT-2000plus和Matrix-2000加速卡,其與Intel的Xeon架構(gòu)類似.ARM架構(gòu)的CPU向量化操作 (asimd)有NEON技術(shù).我們在“天河三號”原型系統(tǒng)上,測試了向量化前后Grid的性能對比情況.

GCC編譯器可以通過-march 選項設置相關(guān)參數(shù)來啟用向量化編譯.首先我們測試在沒有開啟向量化編譯情況下Grid的性能.在每個節(jié)點上,我們使用 1個MPI 進程和64個OpenMP 線程進行測試,每個節(jié)點負責的子格子大小是固定的164.我們記錄程序的計算性能和計算時間,如表3所示,其中,第2列和第3列是計算 Dslash和Wilson fermion時單個節(jié)點的總計算性能(Dslash是格點QCD 中的最重要的算符,Wilson fermion是格點QCD 中的一種費米子),第4列和第5列則是計算所消耗的總時間和通訊時間.

表3 ARM平臺未開啟矢量化優(yōu)化測試

在開啟向量化(asimd)后,Grid的性能測試結(jié)果如表4,各列含義與表3相同.

表4 ARM平臺下Grid開啟ASIMD向量化測試結(jié)果

作為對比,結(jié)合表3和表4,在圖6中畫出了表示性能的Wilson列數(shù)值隨節(jié)點數(shù)變化的曲線圖,其中w/o asmid代表沒有開啟asimd;而w asimd表示開啟了asimd.

圖6 Grid弱可擴展性測試

從圖6可以看出,一方面,隨著節(jié)點數(shù)增加,單個節(jié)點計算能力慢慢下降,但總計算能力呈線性增加,這表明Grid 有很好的弱擴展性;另一方面,目前我們只是采取了基本的自動化向量化編譯,可以看出矢量化后程序總體計算性能有一定提升,而且加入矢量化之后,同樣具有非常好的弱擴展性,可以預測在更大節(jié)點數(shù)時,也能能夠保持較好的計算效率.

矢量化前后的計算時間、通訊時間和總時間隨節(jié)點數(shù)的變化情況如圖7所示.

圖7 ARM下Grid開啟ASIMD向量化前后Wilson算符的計算、通訊時間和總時間對比圖.

可以看到的是,兩種情況下計算時間基本一樣,減少明顯的是通訊時間.這是可以理解的,在開啟向量化之前,MPI間進行同步和歸約時會有瓶頸,每次迭代計算時有大量的同步;而在向量化后,MPI會一次同步和歸約多個元素,因此同步和歸約的次數(shù)會大大減少,這一部分不再是瓶頸,因此時間也相對減少了很多.

3 結(jié)論與展望

本文圍繞格點QCD計算中采用數(shù)據(jù)并行模型提升上層應用計算效率的關(guān)鍵問題,開展了計算算法的理論分析,針對典型應用軟件Grid 開展了面向數(shù)據(jù)并行計算模型的軟件架構(gòu)分析,并且結(jié)合具體計算平臺,測試了Grid應用軟件的性能,包括向量化特征和大規(guī)?？蓴U展性分析.最終闡述了數(shù)據(jù)并行計算模式,對格點QCD計算的有效性和重要性.