張仁強(qiáng) 蔣翔宇 俞炯弛 曾充 宮明 徐順

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

2) (浙江大學(xué)高分子科學(xué)與工程學(xué)系, 杭州 310058)

3) (浙江大學(xué)計算機(jī)科學(xué)與技術(shù)學(xué)院, 杭州 310058)

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

5) (中國科學(xué)院大學(xué)物理科學(xué)學(xué)院, 北京 100049)

格點量子色動力學(xué)(格點QCD)是一種以量子色動力學(xué)為基礎(chǔ), 被廣泛應(yīng)用于強(qiáng)相互作用相關(guān)計算的理論, 作為一種可以給出精確可靠理論結(jié)果的研究方法, 近年來隨著計算機(jī)能力的提升, 正在發(fā)揮著越來越重要的作用.蒸餾算法是格點QCD中計算強(qiáng)子關(guān)聯(lián)函數(shù)的一種重要數(shù)值方法, 可以提高所計算物理量的信噪比.但用它來構(gòu)造關(guān)聯(lián)函數(shù)時, 同樣面臨著數(shù)據(jù)量大和數(shù)據(jù)維數(shù)多的問題, 需要進(jìn)一步提升計算效率.本文開發(fā)了一套利用蒸餾算法產(chǎn)生夸克雙線性算符的關(guān)聯(lián)函數(shù)的程序, 利用MPI (message passing interface,消息傳遞接口, https://www.open-mpi.org), OpenMP (open multi-processing, 共享存儲并行) 和SIMD (single instruction multiple data, 單指令多數(shù)據(jù)流)多級別優(yōu)化技術(shù)解決其中計算性能瓶頸問題.對程序進(jìn)行了多方面的測試, 結(jié)果表明本文的設(shè)計方案能夠支持大規(guī)模的計算, 在強(qiáng)擴(kuò)展性測試下512個進(jìn)程并行計算仍能達(dá)到70%左右的效率, 大大提升了計算關(guān)聯(lián)函數(shù)的能力.

1 引 言

在自然界中有四種基本相互作用: 電磁相互作用、弱相互作用、強(qiáng)相互作用、引力相互作用, 它們決定了目前所知世界的物質(zhì)運(yùn)動規(guī)律.比如摩擦力實際上是電磁相互作用的宏觀效果; 而強(qiáng)相互作用決定了強(qiáng)子的基本性質(zhì), 如強(qiáng)子的質(zhì)量[1]、衰變寬度[2]、形狀因子[3]等, 同時將質(zhì)子和中子結(jié)合在一起組成原子核的核力, 實際上是剩余強(qiáng)相互作用,這種力決定了原子核的聚變和裂變以及質(zhì)量等性質(zhì)[4,5].因此研究強(qiáng)相互作用具有非常重要的意義.描述強(qiáng)相互作用的理論是量子色動力學(xué)(quantum chromodynamics, QCD)[6], 這種理論與電磁相互作用不同, 它在低能標(biāo)是強(qiáng)耦合的, 這意味著無法用微擾展開的方式去近似計算強(qiáng)子的低能性質(zhì).目前也發(fā)展出了一些研究低能QCD的方法, 如手征微擾論[7]、大Nc展開法[8]等, 這些方法都對理論做了一定程度的修改, 都依賴于模型.而格點QCD則完全不需要任何的模型假設(shè), 它直接從QCD出發(fā), 利用蒙特卡羅方法進(jìn)行可靠的數(shù)值模擬[9].

格點QCD計算應(yīng)用需要借助大規(guī)模計算來實現(xiàn)組態(tài)的計算模擬, 需要巨大的高性能計算資源支撐, 其計算密集型和數(shù)據(jù)密集型特征對高性能計算技術(shù)的發(fā)展提出了不少挑戰(zhàn).格點QCD是未來美國E級計算機(jī)的主要應(yīng)用之一[10], 格點QCD在我國E級超算的應(yīng)用還處于探索階段, 有些算法方面的移植優(yōu)化工作[11].

在格點QCD中, 通常需要計算如下形式的量[12]:

其中, S [U] 表 示QCD的作用量, U (t,x) 是規(guī)范場,Oobserve代表任意可觀測量.在歐幾里得時空中, 可以將(1)式理解為簡單的抽樣模型

然后 用 蒙特 卡羅 方 法 去 模擬, 將 ρ (x) 理解 為概 率幅.只要產(chǎn)生一組滿足 e-S[U]分布的規(guī)范場 U (t,x) ,任意一個可觀測量 Oobserve就是在這一組規(guī)范場上的統(tǒng)計平均值, 并且可以有效地估計其誤差.本工作研究的可觀測量為夸克雙線性算符的兩點關(guān)聯(lián)函數(shù):

為了得到更好的統(tǒng)計精度, 近似地計算全傳播子, 有人提出了蒸餾算法[16,17].這種算法在近似計算全傳播子時需要的存儲空間比傳統(tǒng)方法小很多,還可以很方便地進(jìn)行算符構(gòu)造, 即便在結(jié)束傳播子的計算之后, 仍然可以改變涂摩(smear)[18,19]波函數(shù), 提供了算符構(gòu)造[20]更多的自由空間.蒸餾算法把傳播子的計算轉(zhuǎn)換到拉普拉斯算符的本征空間.本征方程表示為其中

Nv表示采用本征矢量的個數(shù), Nv越大, 對拉普拉斯算符的近似越好, 也就越接近于計算全傳播子.f(λc) 是一個關(guān)于本征值的函數(shù).在蒸餾算法中只需要保存本征值、本征矢量和約化傳播子(perambulator):

另外, 蒸餾算法在進(jìn)行算符構(gòu)造時, 只需要對本征矢量的內(nèi)積進(jìn)行操作, 不需要去更改約化傳播子.這意味著, 只需要計算一次約化傳播子, 便可以重復(fù)地利用; 配合上對本征矢量的各種操作, 就可以構(gòu)造各種各樣的算符.這種靈活性和復(fù)用性是傳統(tǒng)方法不具備的; 而且非連通部分的計算可以直接利用連通部分的中間結(jié)果, 不需要進(jìn)行額外的計算.因此, 蒸餾算法具有很多傳統(tǒng)方法無法企及的優(yōu)勢.

雖然蒸餾算法有很多優(yōu)點, 但是組合約化傳播子和拉普拉斯算符本征系統(tǒng)來構(gòu)建關(guān)聯(lián)函數(shù)過程的計算復(fù)雜度正比于將 Nv=10 和Nv=100 進(jìn) 行 對 比, Nv=100 的 計 算 復(fù) 雜 度 是Nv=10 的10000倍.而越是精確的計算越是要求Nv要足夠大, 因此并行計算以及相應(yīng)的優(yōu)化是很有必要的.我們根據(jù)蒸餾算法的特性, 開發(fā)實現(xiàn)了一套利用約化傳播子和拉普拉斯算符本征系統(tǒng)計算關(guān)聯(lián)函數(shù)的程序, 并進(jìn)行了多方面的計算優(yōu)化和測試.

本文將探討蒸餾算法中利用約化傳播子和拉普拉斯算符本征系統(tǒng)計算夸克雙線性算符的關(guān)聯(lián)函數(shù)的實現(xiàn)方法, 結(jié)合算法的計算特點提出了代碼優(yōu)化的設(shè)計方法, 以提高關(guān)聯(lián)函數(shù)的計算效率.首先介紹通過蒸餾算法構(gòu)造關(guān)聯(lián)函數(shù)的基本原理; 其次介紹具體的程序優(yōu)化方案; 在第3部分, 給出了測試結(jié)果和分析; 最后在第4部分進(jìn)行總結(jié).

2 核心算法

為了降低問題的復(fù)雜性, 以關(guān)聯(lián)函數(shù)的連通部分為例, 解釋如何利用蒸餾算法計算關(guān)聯(lián)函數(shù)并且優(yōu)化程序.非連通部分的關(guān)聯(lián)函數(shù)的計算可以利用連通部分的計算中間結(jié)果得到, 因此只對連通部分進(jìn)行討論.

在格點QCD計算中, 通過蒸餾算法構(gòu)造介子連通部分關(guān)聯(lián)函數(shù)的形式[16]為

其中 c1,c2,c3,c4均為對拉普拉斯算符本征值和本征矢量序號的標(biāo)記. Nv表示本征矢量的個數(shù).α,β,γ,κ 均為物理學(xué)中的旋量指標(biāo), 它們的取值范圍為 { 0,1,2,3} . A ,B 標(biāo)記著物理上的不同算符, 每兩個具有相同量子數(shù)的算符就可以用來構(gòu)造關(guān)聯(lián)函數(shù), 它們會在不同的哈密頓量本征態(tài)上有不同的投影, 對應(yīng)著不同的能量.記關(guān)聯(lián)函數(shù)矩陣[21]為

Nop表示關(guān)聯(lián)函數(shù)矩陣的維數(shù).(6)式中Φ和τ的表達(dá)式為

x ,ω 均為空間指標(biāo).通常τ被稱為約化傳播子, 它是拉普拉斯算符本征空間的傳播子, 對應(yīng)傳統(tǒng)算法的傳播子. a ,b,c,d 標(biāo)記的是物理學(xué)中的色指標(biāo), 它的取值范圍是 { 0,1,2} .Φ的不同定義決定了算符的量子數(shù), 對應(yīng)到真實世界中的不同粒子.實際上,(7)式中的G就是傳統(tǒng)方法中的傳播子, 可以看到由于本征矢量與傳播子的空間指標(biāo)和色指標(biāo)進(jìn)行了求和收縮合并, 使得τ和Φ只剩下了表示本征矢量的指標(biāo)、時間的指標(biāo)和旋量的指標(biāo).注意e-σλc1這一項由格點QCD中的涂摩算法給出. λc1表示第 c1個本征矢量的本征值, 其值越大, 該本征矢量對關(guān)聯(lián)函數(shù)的貢獻(xiàn)就越小, 因此通常選取 Nv個最小的本征值對應(yīng)子空間來近似計算全傳播子.

通過蒸餾算法計算關(guān)聯(lián)函數(shù)的步驟如圖1, 首先解拉普拉斯算符的本征值和本征矢量; 再利用本征矢量, 計算τ, 并保存以便重復(fù)使用; 根據(jù)研究的需要利用本征矢量和本征值構(gòu)造Φ; 最后按照(6)式, 利用Φ和τ計算相應(yīng)的關(guān)聯(lián)函數(shù).本文工作主要涉及后兩步.計算λ, v, τ的過程只需進(jìn)行一次.在后續(xù)的計算中根據(jù)不同的研究需要, 對本征矢v進(jìn)行相應(yīng)的操作再與約化傳播子τ進(jìn)行縮并可以構(gòu)造出相應(yīng)的關(guān)聯(lián)函數(shù).重復(fù)利用λ, v, τ.相較于傳統(tǒng)方法在每次進(jìn)行算符構(gòu)造時都需要進(jìn)行傳播子求解, 蒸餾算法構(gòu)造算符更加方便, 在有了約化傳播子和本征值、本征矢量后, 計算關(guān)聯(lián)函數(shù)更經(jīng)濟(jì).同時, 由于產(chǎn)生約化傳播子的過程與構(gòu)造算符的過程相對獨立, 蒸餾算法可以在研究課題沒有確定的情況下, 不斷利用可用的計算資源產(chǎn)生約化傳播子和本征值、本征矢量, 合理利用資源.

圖1 利用蒸餾算法計算關(guān)聯(lián)函數(shù)的流程Fig.1.The procedure of computing correlators via distillation method.

(8)式中Φ的計算實際上可以分解為兩部分,Φ=ΦC?ΦG. ΦC代表拉普拉斯本征空間的計算,ΦG代表旋量空間相關(guān)的計算. ΦC和 ΦG的計算是互相獨立的, 并且 ΦG是一些常數(shù)矩陣的組合, 不需要額外進(jìn)行計算.在程序計算時, 將 ΦG固定在程序內(nèi)部, 可變的 ΦC由外部文件讀入, 然后在根據(jù)需要將 ΦC和 ΦG通過一定方式組合起來得到所需的具有確定量子數(shù)的Φ.

本征矢量的個數(shù) Nv的取值影響著物理信號,通常需要 Nv足夠大以便于對盡可能多的情況保持良好的信號, 而蒸餾算法從約化傳播子計算關(guān)聯(lián)函數(shù)的復(fù)雜度雖然相比傳統(tǒng)方法計算一個全傳播子的關(guān)聯(lián)函數(shù)的計算量要小很多, 但仍然需要處理次乘法.通常 Nop的取值是 O (10) 級別, Nv是 O (100) 級 別, Nt也是 O(100)級別.

為加速計算, 在進(jìn)程級別、線程級別和指令集級別分別進(jìn)行了算法的計算優(yōu)化.

3 優(yōu)化方案

對(6)式進(jìn)行計算約化, 以減少不必要的計算開銷.結(jié)合蒸餾算法計算的特性, 選擇時間維度上的切片操作, 切割后的數(shù)據(jù)處理相對獨立, 使用MPI多線程方式實現(xiàn)時, 避免了數(shù)據(jù)通信成為瓶頸的可能; 另外高維約化傳播子和Φ的數(shù)據(jù)被劃分到多個計算節(jié)點上處理, 也避免了單個計算節(jié)點的內(nèi)存容量限制.為了進(jìn)一步提高計算并行度, 在MPI進(jìn)程中, 實現(xiàn)了基于共享內(nèi)存特性的線程級并行計算, 最大可能地提升了計算節(jié)點中資源的利用率.最后, 由于涉及到大量復(fù)雜的復(fù)數(shù)運(yùn)算操作,開展了數(shù)據(jù)結(jié)構(gòu)的矢量化適配, 實現(xiàn)了SIMD指令級的并行計算.

3.1 計算約化

關(guān)聯(lián)函數(shù)計算表達(dá)式(6)式變形成

注意到括號中的 c2和β指標(biāo)下的部分計算用T矩陣表示, 那么關(guān)聯(lián)函數(shù)可以寫成

在計算關(guān)聯(lián)函數(shù)之前, 可以先計算出T矩陣, 然后利用它來計算不同的關(guān)聯(lián)函數(shù).因為T的計算量∝Nop×Nv, 于是關(guān)聯(lián)函數(shù)的計算量可見計算復(fù)雜度比原來小很多, 在 Nv和 Nop越大時越明顯.實現(xiàn)計算約化后的程序流程如圖2.

圖2 含計算約化的關(guān)聯(lián)函數(shù)計算的流程.其中 T A 和TB 表示兩個中間計算量.利用中間量的計算減少了總體的計算量, 讓計算量從 變成 極大地減少了計算量Fig.2.The flowchart of computing correlation function.TA and T B are two intermidiate quantities.After introducted intermediate quantities, the computation consumption is highly reducted to

3.2 進(jìn)程和線程級優(yōu)化

注意到關(guān)聯(lián)函數(shù)的計算關(guān)于各個時間片t和 t′獨立, 因此選擇對時間維度進(jìn)行切分, 這樣就能夠在各個進(jìn)程進(jìn)行獨立的計算, 完成各個t, t′所負(fù)責(zé)的那部分?jǐn)?shù)據(jù)進(jìn)行收縮合并, 無需考慮通信開銷.程序?qū)崿F(xiàn)是選擇MPI并行計算方式, 每個MPI進(jìn)程負(fù)責(zé)各個獨立的時間片, 在MPI進(jìn)程中進(jìn)一步啟用OpenMP線程共享內(nèi)存的并行計算, OpenMP多線程沒有跨計算節(jié)點的通信, 因此不會增加節(jié)點間通信.

在程序?qū)崿F(xiàn)時, 選擇依次對t和 t′切分.先對t進(jìn)行切分, 當(dāng)進(jìn)程數(shù)大于 Nt時再考慮對 t′進(jìn)行切分.以 Nt=128 為例, 如果進(jìn)程數(shù)等于64, 那么第一個進(jìn)程就要計算 t =0,1 , t′=0,1,···,127 的關(guān)聯(lián)函數(shù), 第二個進(jìn)程就要計算 t =2,3 , t′=0,1,···,127 的關(guān)聯(lián)函數(shù)的數(shù)據(jù).如果進(jìn)程數(shù)等于256, 那么第一個進(jìn)程就要計算 t =0 , t′=0,1,···,63 的關(guān)聯(lián)函數(shù), 第二個進(jìn)程就要計算 t =0 , t′=64,65,···,127 的數(shù)據(jù), 以此類推.這樣進(jìn)行劃分的另外一個優(yōu)勢是, 由于關(guān)聯(lián)函數(shù)具有隨 t′指數(shù)衰減的性質(zhì),以及對于t具有平移不變性, 可以利用這兩個性質(zhì)來檢查計算結(jié)果的正確性.

由于約化傳播子τ以及Φ的維度 ∝ Nv2, 當(dāng) Nv的取值接近1000, Nt=128 時, 按雙精度存儲,τ的數(shù)據(jù)量大約1 TB, 因此在MPI實現(xiàn)時使用并行IO的數(shù)據(jù)讀入方式, 并行IO也是按時間片劃分方式, 和MPI并行計算方案保持一致.注意到Φ只有一個時間指標(biāo), 而τ有兩個時間指標(biāo).在實際計算中, 當(dāng)并行進(jìn)程數(shù) Np≤Nt時, 只對τ其中一個時間指標(biāo)進(jìn)行切分, 分配到各個進(jìn)程上去, 而Φ不切分; 當(dāng)并行進(jìn)程數(shù) Np＞Nt時, 同時對τ的兩個時間維度進(jìn)行切分, 而切片τ的計算只依賴部分Φ, 即每個并行進(jìn)程上都只需Φ的部分矩陣元,因此Φ切片依賴于τ的切分方式.其流程如圖3所示.

圖3 按照時間切分實現(xiàn)并行計算的方式, 根據(jù) N p 與 Nt的相對大小, 由于數(shù)據(jù)的特性, 對τ和Φ按情況采用不同的切分方法.Fig.3.Data sgemented according to time.Two conditions are considered which decided how τ and Φ are treated because of the feature of data.

3.3 指令級優(yōu)化

在計算中, 考慮到形如

的張量縮并過程, 可以將其簡化為與縮并指標(biāo)相對應(yīng)的矩陣乘法.上述的兩個表達(dá)式都可以看作形如(4Nv,4Nv) 的矩陣乘法.如將原本按照 (4,4,Nv,Nv)行優(yōu)先存儲的數(shù)據(jù)轉(zhuǎn)換為 ( 4Nv,4Nv) 行優(yōu)先存儲的格式, 那么可得如下的計算方式:

行優(yōu)先存儲格式下按(13)式作矩陣乘法時, B矩陣無法實現(xiàn)連續(xù)的內(nèi)存訪問, 因此將B矩陣轉(zhuǎn)置, 計算變?yōu)?/p>

則計算中在循環(huán)k時可以內(nèi)存連續(xù)地訪問矩陣的元素.

在計算中, 每個數(shù)字實際上是復(fù)數(shù), 而由于復(fù)數(shù)乘法規(guī)則復(fù)雜, 在不更改數(shù)據(jù)排列的情況下,編譯器無法保證利用SIMD矢量運(yùn)算提升復(fù)數(shù)計算性能.因此將矩陣的實部和虛部分離存儲,利用運(yùn)算法則將1個復(fù)矩陣的乘法分解為4個實矩陣乘法和2個實矩陣加法, 形式為 Re(AB)=Re(A)Re(B)-Im(A)Im(B) , Im(AB)=Re(A)×Im(B)+Im(A)Re(B) .

在上述調(diào)整下, 變量的內(nèi)存訪問具有連續(xù)性,有助于編譯器自動完成SIMD矢量計算優(yōu)化.通過確認(rèn)編譯得到的匯編代碼, 編譯器確實啟用SIMD矢量運(yùn)算指令.當(dāng)然, 也根據(jù)測試平臺手動調(diào)用SIMD相關(guān)函數(shù)和指令實現(xiàn)了計算, 經(jīng)過測試相比編譯器自動優(yōu)化提升不到10%, 考慮到硬編碼的繁瑣之處和兼容問題, 我們最終采用的方式是修改數(shù)據(jù)格式和實現(xiàn)復(fù)數(shù)計算后, 通過編譯器自動優(yōu)化完成SIMD向量化.

4 測試結(jié)果

為測試 加速效果, 在 Nt=128 , Ns=16 的格子上進(jìn)行測試.組態(tài)的海里只包含兩味粲夸克.計算過程中采用雙精度. a-t1=9.6Gev .分別對SIMD,MPI, OpenMP的加速效果進(jìn)行測試, 并對物理計算的結(jié)果進(jìn)行展示.

4.1 SIMD計算加速測試

為了驗證SIMD計算加速效果, 選擇一個含有Nv=10 個本征向量的算例, 按照(9)式來計算關(guān)聯(lián)函數(shù)矩陣( 5 ×5 維).測試的處理器支持FMA和AVX SIMD指令集.該CPU處理器包含8個物理核心, 可通過超線程技術(shù)實現(xiàn)的16進(jìn)程并行.計算時對整個流程按步驟分別計時, 具體步驟包括: 初始化、數(shù)據(jù)IO讀取、計算T中間量、利用T矩陣計算關(guān)聯(lián)函數(shù).為了實現(xiàn)使用SIMD優(yōu)化效果, 還對照測試了直接調(diào)用標(biāo)準(zhǔn)庫中的復(fù)數(shù)計算函數(shù)的程序設(shè)計.按步驟計時的結(jié)果如圖4, 使用SIMD前后計算性能如圖5.

圖4 使用SIMD優(yōu)化前后各階段耗時對比.I/O代表圖1中第一步和第二步的時間, Calc.prepare代表圖1中第三步的時間, Calc.result代表圖1中第四步的時間, Others代表圖1中第五步的時間, Init代表程序初始化的時間.圖例SIMD表示啟用了AVX形式的SIMD計算性能, 而Complex表示程序直接調(diào)用標(biāo)準(zhǔn)庫中的復(fù)數(shù)計算函數(shù)(此處未使用SIMD計算).其中16個MPI進(jìn)程并行計算的結(jié)果是在超線程計算狀態(tài)下獲得Fig.4.The cost of time of program's each part to see the effects of SIMD.I/O labels the time of first step and second step in Fig.1, Calc.prepare labels the time of the third step in Fig.1, Calc.result labels the time of the fourth step in Fig.1, Others labels the time of the fifth step in Fig.1.Init labels the time of initialization.SIMD in the picture means SIMD optimization was adopted and Complex in the picture means the stdandard library of complex computation was used.And hyper-threading technology was used for 16 MPI process.

圖5 使用SIMD優(yōu)化前后性能對比.圖例SIMD表示啟用了AVX形式的SIMD計算性能, 而Complex表示程序直接調(diào)用標(biāo)準(zhǔn)庫中的復(fù)數(shù)計算函數(shù)(此處未使用SIMD計算).其中, 在SIMD啟用時16個超線程計算結(jié)果未參與數(shù)據(jù)擬合Fig.5.The cost of time of program's each part to see the effects of SIMD.SIMD in the picture means SIMD optimization was adopted and Complex in the picture means the stdandard library of complex computation was used.And hyper-threading technology was used for 16 MPI process.

程序在啟用SIMD指令后, 在不考慮超線程的時, 運(yùn)行時間變?yōu)樵瓉淼?/4.Calc.prepare的時間變?yōu)樵瓉淼?/6.Calc.prepare是程序的熱點, 此處計算程序中間量, 計算時間隨Nv的變化最明顯.

4.2 OpenMP計算加速測試

現(xiàn)測試OpenMP計算加速效果, 沿用了上一節(jié)測試的算例, 處理器開啟CPU的超線程, 并關(guān)閉SIMD加速.同樣按步驟計時的測試結(jié)果如圖6.

圖6 使用OpenMP優(yōu)化前后耗時對比.圖例如圖4.圖例Serial表示串行版本, 即未開啟OpenMP多線程和MPI多進(jìn)程Fig.6.The effects of OpenMP optimization was showed.Legends are the same as 4.Serial lables the results of serial program which means no OpenMP and MPI was adopted.

從結(jié)果可見, 程序采用OpenMP線程加速方式有性能提升, 使用16個超線程加速相對單核獲得了約10倍的加速.同時采用OpenMP線程的計算效率與采用MPI進(jìn)程的計算效率幾乎相同.

4.3 MPI強(qiáng)擴(kuò)展性測試

為了測試MPI并行計算的強(qiáng)擴(kuò)展性, 選擇一個含有 Nv=100 個本征向量的大算例, 此時輸入文件體積達(dá)到40 GB, 采用不同的計算核數(shù)分別按照(9)式來計算關(guān)聯(lián)函數(shù)矩陣( 5 ×5 維).測試的處理器支持FMA和AVX SIMD指令集, 最大MPI進(jìn)程達(dá)到512個.同樣按初始化、數(shù)據(jù)IO讀取、計算T中間量、利用T矩陣計算關(guān)聯(lián)函數(shù)等步驟計時.結(jié)果如圖7.MPI并行的效果明顯, 適合大規(guī)模并行.MPI進(jìn)程數(shù)的增加, 效果體現(xiàn)在Calc.prepare部分, 也就是(11)式中間量的計算加速效果最明顯.由于各部分加速的效果不同, 當(dāng)512個進(jìn)程時,Calc.prepare的計算時間與I/O的時間差別不大.I/O隨進(jìn)程數(shù)的變化時間改變并不明顯, 耗時可以忽略.

圖7 MPI并行強(qiáng)擴(kuò)展性測試.隨著MPI進(jìn)程數(shù)增加, 計算時間成比例減少.圖例如圖4Fig.7.MPI parallelism in strong scale tests.The cost time decrease with MPI process numbers.Legends are the same as Fig.4.

計算程序在強(qiáng)擴(kuò)展性測試中的并行計算效率如圖8, 從結(jié)果可見, 總體并行計算效率很高, 在512個進(jìn)程的情況下仍能達(dá)到70%左右的效率.因此程序適合大規(guī)模計算加速.

圖8 MPI并行強(qiáng)擴(kuò)展性測試，不同MPI進(jìn)程數(shù)的測試相對于16個進(jìn)程的計算效率Fig.8.MPI parallelism in strong scale tests.The efficiency of strong scale tests compared with 16 MPI processes.

4.4 MPI弱擴(kuò)展性測試

為了測試MPI并行計算的弱擴(kuò)展性, 選擇幾個不同規(guī)模的算例, 其進(jìn)程數(shù)分別為 Np=2,16,128,1024 , 關(guān)聯(lián)函數(shù)矩陣維度仍然為 5 ×5 .測試的處理器支持FMA和AVX SIMD指令集, 對不同算例等比例使用不同規(guī)模的并行進(jìn)程數(shù), 使得理論上每個進(jìn)程的負(fù)載相同.按步驟計時的結(jié)果如圖9,忽略I/O部分, 其他部分表現(xiàn)良好.但I(xiàn)/O部分當(dāng)進(jìn)程數(shù)等于1024時, 耗時出現(xiàn)了明顯的增加, 并與熱點部分(Calc.prepare)的耗時可比.這受限于硬件的性能.在當(dāng)前測試中, 當(dāng) Np≤128 時程序在弱擴(kuò)展性測試中的并行效率較高, 但 Np=1024 的大規(guī)模測試中展現(xiàn)出了一個因系統(tǒng)資源限制產(chǎn)生的問題, 由于磁盤帶寬是一定的, I/O請求數(shù)激增, I/O部分時間大大增加, 甚至超過計算中間變量T矩陣(Calc.prepare)部分成為了新的熱點; 但并行I/O能夠解決串行讀入數(shù)據(jù)時內(nèi)存不足的問題.由于計算量是隨著 Nv的三次方增加的, 而數(shù)據(jù)量大小是隨著 Nv平方增加的, 這會導(dǎo)致如果要維持每個進(jìn)程的計算量不變相應(yīng)的數(shù)據(jù)量就會變少, 讀入的數(shù)據(jù)會隨著 Nv的增加而變得零碎, 最終受限于I/O.

圖9 MPI并行弱擴(kuò)展性測試. N p 表示使用的并行進(jìn)程數(shù), N v 表示本征向量數(shù).圖例說明同圖4Fig.9.MPI parallelism in weak scale tests. N p represents the number of process, N v represents the number of eigenvectors.Legends are the same as Fig.4.

4.5 物理結(jié)果

為了檢驗結(jié)果的正確性, 對物理結(jié)果進(jìn)行了測試.構(gòu)造了三個JPC量子數(shù)為 0-+的算符, 分別計算它們的關(guān)聯(lián)函數(shù)并分析有效質(zhì)量,這里 C (t) 表示某種量子數(shù)算符的關(guān)聯(lián)函數(shù).通過比較這三個算符所得的關(guān)聯(lián)函數(shù)的有效質(zhì)量, 在時間片比較大的時候趨于同一個平臺, 與傳統(tǒng)方法所得結(jié)果一致, 見圖10.

圖10 做變分前的結(jié)果.在時間比較小時, 三個算符得到的有效質(zhì)量的行為有很大差別, 證明在不同態(tài)的投影是不同的, 意味著變分會有一定的效果.在時間較大時, 三個算符的有效質(zhì)量趨于同一平臺, 說明它們的量子數(shù)是相同的,可以用來變分.traditional method表示第一個算符通過傳統(tǒng)方法所得到的有效質(zhì)量, 作為蒸餾算法的參照Fig.10.Results before variation.The behaviors of the effective mass of these three operators are very different and it means variational analysis would give good results.When time is large enough, these three operators approach to the same plateau so that they should have the same quantum numbers.traditial method label the effective mass of first operator throgh traditional mehtod, which can is matched with distillation method.

因為三個算符在不同本征能量態(tài)上的投影有很大不同, 如果中間計算沒有錯誤, 通過變分分析[21]可以近似得到三個主要投影到最低的三個能量本征態(tài)的算符, 變分方程如下:

變分后, 可以得到本征矢 Vi, 對應(yīng)第i個態(tài) mi,將 CA,B與 Vi組合, 可以得到對應(yīng)到相應(yīng)態(tài)的關(guān)聯(lián)函數(shù)

并且利用這三個算符所計算出的有效質(zhì)量要能夠相互區(qū)別開, 在時間片比較大的時候趨于不同的平臺.因此可以此來驗證計算的正確性.

利用上述三個算符構(gòu)造關(guān)聯(lián)函數(shù)矩陣, 做變分分析, 得到如圖11所示的結(jié)果.變分后得到三個不同的態(tài), 分別對應(yīng)著基態(tài)、第一激發(fā)態(tài)、第二激發(fā)態(tài).利用變分本征矢量組合后的關(guān)聯(lián)函數(shù)所得的有效質(zhì)量在時間較大時, 分別趨于不同的平臺, 證明變分將原來的三個算符分別組合成了三個主要投影到最低三個態(tài)的優(yōu)化算符.第二激發(fā)態(tài)的信號略差.變分效果明顯, 再次驗證了結(jié)果的正確性.

圖11 做變分后的結(jié)果Fig.11.Results after variation.

5 結(jié)論與展望

本文圍繞格點QCD蒸餾算法中關(guān)聯(lián)函數(shù)的計算優(yōu)化問題, 提出了一套加速程序計算的有效解決方案, 結(jié)合理論分析采用了MPI, OpenMP和SIMD多級別的計算優(yōu)化技術(shù).文中首先介紹了蒸餾算法中關(guān)聯(lián)函數(shù)計算中各個模塊的優(yōu)化設(shè)計, 并對程序的性能進(jìn)行了多方面的分析和測試.實驗結(jié)果表明蒸餾算法計算關(guān)聯(lián)函數(shù)效率的顯著提升, 并且驗證了物理結(jié)果的正確性.