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

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

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

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

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

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

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

1 引 言

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

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

在格點(diǎn)QCD中, 通常需要計(jì)算如下形式的量[12]:

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

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

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

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

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

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

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

2 核心算法

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

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

其中 c1,c2,c3,c4均為對拉普拉斯算符本征值和本征矢量序號的標(biāo)記. Nv表示本征矢量的個(gè)數(shù).α,β,γ,κ 均為物理學(xué)中的旋量指標(biāo), 它們的取值范圍為 { 0,1,2,3} . A ,B 標(biāo)記著物理上的不同算符, 每兩個(gè)具有相同量子數(shù)的算符就可以用來構(gòu)造關(guān)聯(lián)函數(shù), 它們會(huì)在不同的哈密頓量本征態(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)到真實(shí)世界中的不同粒子.實(shí)際上,(7)式中的G就是傳統(tǒng)方法中的傳播子, 可以看到由于本征矢量與傳播子的空間指標(biāo)和色指標(biāo)進(jìn)行了求和收縮合并, 使得τ和Φ只剩下了表示本征矢量的指標(biāo)、時(shí)間的指標(biāo)和旋量的指標(biāo).注意e-σλc1這一項(xiàng)由格點(diǎn)QCD中的涂摩算法給出. λc1表示第 c1個(gè)本征矢量的本征值, 其值越大, 該本征矢量對關(guān)聯(lián)函數(shù)的貢獻(xiàn)就越小, 因此通常選取 Nv個(gè)最小的本征值對應(yīng)子空間來近似計(jì)算全傳播子.

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

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

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

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

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

3 優(yōu)化方案

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

3.1 計(jì)算約化

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

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

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

圖2 含計(jì)算約化的關(guān)聯(lián)函數(shù)計(jì)算的流程.其中 T A 和TB 表示兩個(gè)中間計(jì)算量.利用中間量的計(jì)算減少了總體的計(jì)算量, 讓計(jì)算量從 變成 極大地減少了計(jì)算量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ù)的計(jì)算關(guān)于各個(gè)時(shí)間片t和 t′獨(dú)立, 因此選擇對時(shí)間維度進(jìn)行切分, 這樣就能夠在各個(gè)進(jìn)程進(jìn)行獨(dú)立的計(jì)算, 完成各個(gè)t, t′所負(fù)責(zé)的那部分?jǐn)?shù)據(jù)進(jìn)行收縮合并, 無需考慮通信開銷.程序?qū)崿F(xiàn)是選擇MPI并行計(jì)算方式, 每個(gè)MPI進(jìn)程負(fù)責(zé)各個(gè)獨(dú)立的時(shí)間片, 在MPI進(jìn)程中進(jìn)一步啟用OpenMP線程共享內(nèi)存的并行計(jì)算, OpenMP多線程沒有跨計(jì)算節(jié)點(diǎn)的通信, 因此不會(huì)增加節(jié)點(diǎn)間通信.

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

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

圖3 按照時(shí)間切分實(shí)現(xiàn)并行計(jì)算的方式, 根據(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)化

在計(jì)算中, 考慮到形如

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

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

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

在計(jì)算中, 每個(gè)數(shù)字實(shí)際上是復(fù)數(shù), 而由于復(fù)數(shù)乘法規(guī)則復(fù)雜, 在不更改數(shù)據(jù)排列的情況下,編譯器無法保證利用SIMD矢量運(yùn)算提升復(fù)數(shù)計(jì)算性能.因此將矩陣的實(shí)部和虛部分離存儲(chǔ),利用運(yùn)算法則將1個(gè)復(fù)矩陣的乘法分解為4個(gè)實(shí)矩陣乘法和2個(gè)實(shí)矩陣加法, 形式為 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ù)性,有助于編譯器自動(dòng)完成SIMD矢量計(jì)算優(yōu)化.通過確認(rèn)編譯得到的匯編代碼, 編譯器確實(shí)啟用SIMD矢量運(yùn)算指令.當(dāng)然, 也根據(jù)測試平臺(tái)手動(dòng)調(diào)用SIMD相關(guān)函數(shù)和指令實(shí)現(xiàn)了計(jì)算, 經(jīng)過測試相比編譯器自動(dòng)優(yōu)化提升不到10%, 考慮到硬編碼的繁瑣之處和兼容問題, 我們最終采用的方式是修改數(shù)據(jù)格式和實(shí)現(xiàn)復(fù)數(shù)計(jì)算后, 通過編譯器自動(dòng)優(yōu)化完成SIMD向量化.

4 測試結(jié)果

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

4.1 SIMD計(jì)算加速測試

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

圖4 使用SIMD優(yōu)化前后各階段耗時(shí)對比.I/O代表圖1中第一步和第二步的時(shí)間, Calc.prepare代表圖1中第三步的時(shí)間, Calc.result代表圖1中第四步的時(shí)間, Others代表圖1中第五步的時(shí)間, Init代表程序初始化的時(shí)間.圖例SIMD表示啟用了AVX形式的SIMD計(jì)算性能, 而Complex表示程序直接調(diào)用標(biāo)準(zhǔn)庫中的復(fù)數(shù)計(jì)算函數(shù)(此處未使用SIMD計(jì)算).其中16個(gè)MPI進(jìn)程并行計(jì)算的結(jié)果是在超線程計(jì)算狀態(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計(jì)算性能, 而Complex表示程序直接調(diào)用標(biāo)準(zhǔn)庫中的復(fù)數(shù)計(jì)算函數(shù)(此處未使用SIMD計(jì)算).其中, 在SIMD啟用時(shí)16個(gè)超線程計(jì)算結(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指令后, 在不考慮超線程的時(shí), 運(yùn)行時(shí)間變?yōu)樵瓉淼?/4.Calc.prepare的時(shí)間變?yōu)樵瓉淼?/6.Calc.prepare是程序的熱點(diǎn), 此處計(jì)算程序中間量, 計(jì)算時(shí)間隨Nv的變化最明顯.

4.2 OpenMP計(jì)算加速測試

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

圖6 使用OpenMP優(yōu)化前后耗時(shí)對比.圖例如圖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個(gè)超線程加速相對單核獲得了約10倍的加速.同時(shí)采用OpenMP線程的計(jì)算效率與采用MPI進(jìn)程的計(jì)算效率幾乎相同.

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

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

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

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

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

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

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

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

圖10 做變分前的結(jié)果.在時(shí)間比較小時(shí), 三個(gè)算符得到的有效質(zhì)量的行為有很大差別, 證明在不同態(tài)的投影是不同的, 意味著變分會(huì)有一定的效果.在時(shí)間較大時(shí), 三個(gè)算符的有效質(zhì)量趨于同一平臺(tái), 說明它們的量子數(shù)是相同的,可以用來變分.traditional method表示第一個(gè)算符通過傳統(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.

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

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

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

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

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

5 結(jié)論與展望

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