周生昌劉衛(wèi)國宋振亞楊曉丹

(1.山東大學 軟件學院,山東 濟南250101;2.自然資源部 第一海洋研究所,山東 青島266061;3.青島海洋科學與技術(shù)試點國家實驗室 區(qū)域海洋動力學與數(shù)值模擬功能實驗室,山東 青島266237;4.海洋環(huán)境科學和數(shù)值模擬自然資源部重點實驗室,山東 青島266061)

研究海洋的手段主要有理論研究、觀測試驗和數(shù)值模擬,其中觀測試驗是數(shù)值模擬和理論研究的基礎(chǔ),理論研究為數(shù)值模擬和觀測試驗設(shè)計提供指導,而數(shù)值模式則集成了人們從觀測數(shù)據(jù)和理論研究中所獲得的知識和理解,不僅可以用來驗證理論假說和彌補觀測資料的不足,而且是能夠最終用于預測未來變化的唯一工具[1]。隨著海洋和氣候變化研究的不斷深入,高性能計算機的發(fā)展,海洋數(shù)值模式正逐步朝著更高分辨率、更多物理過程和更快計算速度的方向發(fā)展[2]。一般來說,水平分辨率每提高1倍,計算量會增大到原來的8～10倍[3]。同時,隨著分辨率的提高,會有越來越多的原次網(wǎng)格不能描述的物理過程被包含進來,而新的次網(wǎng)格過程的參數(shù)化也需要通過將多組數(shù)值試驗結(jié)果與實測數(shù)據(jù)進行比較來進行評估和調(diào)整,這顯然也提高了對計算量和計算能力的要求[4]。因此,充分利用和發(fā)揮高性能計算機的性能來提高海洋數(shù)值模式計算速度已成為模式發(fā)展的一個必要條件。

目前高性能海洋數(shù)值模擬的研究主要采用以下幾種方案提高模式計算能力:1)基于分布式內(nèi)存并行化程序。主要采用MPI(Message Passing Interface)消息傳遞接口使計算任務(wù)合理分布在多個獨立進程的計算單元內(nèi),并通過高速網(wǎng)絡(luò)進行數(shù)據(jù)通信。該方法實現(xiàn)復雜,需考慮邊界交換、負載均衡等問題,雖然可擴展性好,但節(jié)點之間交換會帶來通訊開銷。2)基于共享式內(nèi)存并行化程序。主要采用Open MP等接口將計算任務(wù)分配給多個線程,并通過共享內(nèi)存地址的方式實現(xiàn)數(shù)據(jù)交換。該方法實現(xiàn)相對容易,無通訊開銷,但一般來說無法跨節(jié)點,這限制了其并行規(guī)模。3)基于進程/線程并行化程序。該方法是以上2種并行方式的結(jié)合,實現(xiàn)較為復雜,在異構(gòu)計算機(如CPU+GPU、國產(chǎn)申威處理器等)上效果特別顯著。

上述3種方案主要從并行擴展方面提高模式計算性能,這需要大量的計算資源,且程序在設(shè)計和實現(xiàn)中并未考慮代碼現(xiàn)代化。特別是隨著高性能計算機的發(fā)展,這些方案已經(jīng)無法充分利用和發(fā)揮現(xiàn)代計算機的優(yōu)勢。楊曉丹等[3]在海浪模式MASNUM上開展了代碼現(xiàn)代化優(yōu)化方面的嘗試,通過檢查代碼的規(guī)范性、合理性,并針對語言特性對代碼進行改進,使之更好地適應(yīng)計算機基礎(chǔ)架構(gòu),充分發(fā)揮了現(xiàn)代計算機的計算特點。但是在此海浪模式中,計算量最大的源函數(shù)與相鄰點無關(guān),無分區(qū)間數(shù)據(jù)交換通訊,并且由于分支判斷少,cache命中率高[5],這與海洋環(huán)流模式NEMO(Nucleus for European Modeling of the Ocean)[6]計算特點有較大的差異。NEMO不但存在大量的數(shù)據(jù)交換,而且分支判斷較多。因此如何使用代碼現(xiàn)代化方法優(yōu)化和提高海洋環(huán)流模式的計算性能,非常具有代表性和必要性。

為了充分利用和發(fā)揮現(xiàn)代高性能計算機的性能,進一步促進海洋環(huán)流數(shù)值模式的應(yīng)用和發(fā)展,本文基于代碼現(xiàn)代化的概念,提出了海洋環(huán)流模式代碼現(xiàn)代化優(yōu)化的方案,并在此基礎(chǔ)上對海洋環(huán)流模式NEMO進行了優(yōu)化,分析和探討了代碼現(xiàn)代化優(yōu)化方案在海洋環(huán)流模式應(yīng)用中的前景和局限性,為今后海洋環(huán)流數(shù)值模擬的研究和發(fā)展提供參考和建議。

1 優(yōu)化方案設(shè)計

1.1 優(yōu)化策略

現(xiàn)代高性能計算機是由多核和眾核處理器、高速緩存和內(nèi)存、高帶寬處理器以及它們之間的通信結(jié)構(gòu)和高速I/O組合而成的[7]。為了充分適應(yīng)平臺的現(xiàn)代化結(jié)構(gòu),充分發(fā)揮平臺的性能,需要發(fā)展高性能和現(xiàn)代化的軟件,包括對原始代碼進行改進,以及為現(xiàn)代計算機重新設(shè)計應(yīng)用程序以獲得最大性能,但這些均需要充分利用高性能計算機硬件資源,這是代碼現(xiàn)代化的基礎(chǔ)。

本文基于代碼現(xiàn)代化優(yōu)化策略,并結(jié)合NEMO海洋環(huán)流模式的特點以及現(xiàn)有平臺架構(gòu),制定了6個階段的優(yōu)化策略,如圖1所示。為了保證結(jié)果的準確性,排除平臺因素引起的誤差,每次均取2次試驗結(jié)果的平均值作為最終結(jié)果。

1.2 試驗算例

試驗算例為全球海洋環(huán)流真實算例ORCA1_LIM,該試驗算例只包含海洋動力學、熱力學(OPA模塊)以及海冰動力學(LIM模塊)。模式共積分120個時間步長,水平網(wǎng)格數(shù)目為362×292,垂向網(wǎng)格數(shù)目為75。為了減小I/O對擴展性的影響,本研究中的測試是針對模式計算性能的無I/O測試。

圖1 代碼現(xiàn)代化優(yōu)化策略Fig.1 Code modernization optimization strategy

1.3 測試環(huán)境

試驗測試環(huán)境如表1所示,單節(jié)點內(nèi)共有兩顆Intel Xeon Gold 6252處理器,每顆處理器有24核心,因此單節(jié)點共有48物理核心,96邏輯核心,內(nèi)存為384 GB DDR4。

表1 測試環(huán)境Table 1 The testing environment

2 NEMO優(yōu)化

2.1 定位熱點函數(shù)

為了更準確地獲取影響NEMO模式計算效率的瓶頸,從而更有針對性地對NEMO模式進行優(yōu)化,首先使用Intel Vtune Amplifier工具定位出模式的熱點函數(shù),即耗時較多的函數(shù),如圖2所示。結(jié)果顯示NEMO的熱點函數(shù)極為分散,即熱點函數(shù)較多,耗時分布較為均衡,本研究選取前五個熱點函數(shù)(tra_ldf_iso,lim_rhg,tra_adv_tvd,ldf_slp、nonosc)進行重點優(yōu)化。

圖2 串行模式各函數(shù)計算時間Fig.2 Serial mode function calculates time

2.2 編譯選項優(yōu)化

在性能分析的基礎(chǔ)上,本研究首先使用編譯選項對模式進行優(yōu)化。使用x HOST選項使編譯器生成處理器支持的最高指令集,該選項在本文測試環(huán)境下等效于x MIC-AVX512;使用ipo啟用文件之間的過程間優(yōu)化,使編譯器對單獨文件中定義的函數(shù)執(zhí)行內(nèi)聯(lián)函數(shù)擴展;使用no-prec-div選項對模式的浮點數(shù)除法進行優(yōu)化,將除法改為乘倒數(shù)的形式進行計算,從而提高運算效率;此外,還使用了fp-model fast=2,qopt-dynamic-align,qopt-prefetch等優(yōu)化選項,具體如表2所示。在經(jīng)過編譯選項優(yōu)化后,NEMO的整體加速比可以達到1.21倍。通過對輸出結(jié)果進行分析,以海表面溫度為例,得出O1和O3兩種編譯選項所引起的誤差絕對值不超過9.3×10-6℃,占僅占平均海表面溫度的0.00005%,且誤差絕對值在1.0×10-8℃以下的網(wǎng)格點占比超過99.8%,滿足精度要求。同樣,對其它變量的檢查也滿足精度要求。因此編譯選項優(yōu)化是有一種有效的優(yōu)化手段。

表2 編譯選項優(yōu)化列表Table 2 The list of compilation option optimization

2.3 標量串行優(yōu)化

標量串行優(yōu)化的目的是通過對熱點函數(shù)代碼進行調(diào)整修改,如去除重復計算、減少條件分支、降低循環(huán)嵌套等,確保代碼使用最少的計算量和合適的精度獲得正確的結(jié)果。下面以traldf_iso熱點函數(shù)為例介紹此優(yōu)化過程。

在如下例子中(圖3a),zdk1t和zdkt計算公式是相同的,即zdk1t(ji,jj,jk-1)=zdkt(ji,jj,jk),因此通過將zdk1t從循環(huán)中移除,在計算完zdkt后再將結(jié)果賦值給zdk1t,可以減少重復性計算,如圖3b所示:

圖3 重復計算優(yōu)化前和優(yōu)化后Fig.3 Repeat calculation before optimization after optimization

當從內(nèi)存中訪問數(shù)組元素時,兩次訪問的尋址間隔會影響Cache命中率,降低計算效率。由于Fortran中的數(shù)組是按照列優(yōu)先的方式存儲,因此調(diào)整數(shù)組元素訪問順序,使得相鄰兩次訪問的尋址間隔變小,可以有效提高Cache命中率。調(diào)整數(shù)組元素訪問順序前后分別如圖4a和4b所示。

另外,減少循環(huán)嵌套也能有效提高模式計算效率。在圖5a的循環(huán)中,z2d是臨時數(shù)組,目的是存儲zftu在三維空間沿K軸方向的積分。由于三重循環(huán)計算效率較低,因此利用Fortran提供的forall命令將該循環(huán)減少為一重循環(huán),如圖5b所示。經(jīng)過上述串行與標量優(yōu)化后,NEMO模式整體計算效率提升1.22倍。

圖4 調(diào)整數(shù)組元素訪問順序優(yōu)化前和優(yōu)化后Fig.4 Adjust access order of array element before optimization after optimization

圖5 循環(huán)嵌套優(yōu)化前和優(yōu)化后Fig.5 Loop nesting before optimization after optimization

2.4 SIMD優(yōu)化

SIMD(Single Instruction Multiple Data)即單指令多數(shù)據(jù)流,以同步方式,在同一時間內(nèi)對多個數(shù)據(jù)執(zhí)行同一條指令。矢量化是現(xiàn)代計算機的一個標志,利用自動矢量化,可以顯著提高軟件的計算性能。在本測試環(huán)境下,利用x HOST選項編譯時,已包含了自動矢量化功能。然而自動矢量化需要遵循矢量化規(guī)則,如循環(huán)中不能存在數(shù)據(jù)前后依賴以及跳轉(zhuǎn)語句等。以下是自動矢量化失敗的例子及采用的優(yōu)化策略。

在圖6a所示的代碼中,由于zdit和zdjt都是由指針方式實現(xiàn)的動態(tài)數(shù)組,Fortran指針可以作為別名指向一塊內(nèi)存地址,因此編譯器不能判斷zdit和zdjt指針之間是否存在數(shù)據(jù)依賴,從而無法實現(xiàn)自動矢量化。通過將假定依賴的語句拆分成2個循環(huán),使每個循環(huán)滿足自動矢量化的條件,可以實現(xiàn)自動矢量化(圖6b)。

圖6 SIMD數(shù)據(jù)依賴優(yōu)化前和優(yōu)化后Fig.6 SIMD data dependency before optimization after optimization

另外,由于代碼中存在大量的指針作為臨時數(shù)組,因此在不影響物理過程的前提下,通過將指針修改為動態(tài)數(shù)組,可以使代碼安全地進行自動矢量化,如圖7所示。經(jīng)過SIMD優(yōu)化后,NEMO模式整體計算效率提升至1.23倍。

圖7 SIMD指針優(yōu)化前優(yōu)化后Fig.7 SIMD pointer before optimization after optimization

2.5 內(nèi)存帶寬優(yōu)化

當CPU需要讀取內(nèi)存中的數(shù)據(jù)時,CPU將首先發(fā)出一個由內(nèi)存控制器執(zhí)行的請求,然后再將數(shù)據(jù)由內(nèi)存返回中央處理器,此過程的時間稱為延時周期。為了縮短延時周期,中央處理器和內(nèi)存之間設(shè)置了L1和L2緩存。由于L1和L2均采用靜態(tài)隨機訪問的原理,可以將其理解為內(nèi)存帶寬[8]。對內(nèi)存帶寬進行優(yōu)化,提高緩存命中率,也是提升軟件計算效率的有效方式。

由于海洋環(huán)流模式復雜的物理過程,模式中使用大量的臨時數(shù)組保存計算過程的中間結(jié)果,這就造成了運算時參與計算的數(shù)組過多的問題。當CPU計算時,如果一條語句需要同時讀取大量數(shù)組,由于內(nèi)存帶寬限制,則會減少緩存命中率,降低計算效率。因此減少或避免使用臨時數(shù)組,是緩解內(nèi)存帶寬占用過高的一種優(yōu)化手段。

如圖8a所示,數(shù)組A,B為臨時數(shù)組,運算時首先通過計算為A,B數(shù)組賦值,并在之后的循環(huán)中訪問該數(shù)組。優(yōu)化策略為將臨時數(shù)組A、B去掉,而把計算賦值的功能重構(gòu)為子程序cal_a,cal_b,并在循環(huán)中需要此中間結(jié)果時調(diào)用。

圖8 內(nèi)存帶寬優(yōu)化優(yōu)化前和優(yōu)化后Fig.8 Memory bandwidth before optimization after optimization

對熱點函數(shù)進行分析后,發(fā)現(xiàn)traldf_iso、nonosc和lim_rhg中均存在大量臨時數(shù)組,可以使用上述策略進行優(yōu)化。圖9顯示了優(yōu)化前后各個函數(shù)的計算時間以及其加速比,可以看出優(yōu)化后單個熱點函數(shù)的計算效率有效提升10%～40%,效果顯著。另外,由于NEMO的熱點函數(shù)呈現(xiàn)分散的特點,此優(yōu)化使得模式整體計算效率加速至1.31倍。

圖9 內(nèi)存帶寬優(yōu)化前后函數(shù)時間對比Fig.9 Comparison of function time before and after memory bandwidth optimization

2.6 MPI多進程擴展

NEMO模式總體是一個循環(huán)過程,其中時間是循環(huán)的主序列。循環(huán)從數(shù)據(jù)讀取開始,以結(jié)果輸出為結(jié)束。循環(huán)部分包含模式中的所有核心計算,總計算量占總循環(huán)過程的95%及以上[9]。本文使用的NEMO 3.6版本基于MPI消息傳遞接口編寫,使用MPI分布式內(nèi)存并行計算可以有效的進行任務(wù)劃分,從而提高計算效率。在經(jīng)過以上串行優(yōu)化之后,接下來利用Intel Trace Analyzer Collector工具對NEMO進行多進程擴展測試,以觀察模式的負載均衡性。

圖10顯示了模式在48進程時的負載均衡性,從結(jié)果可以看出,計算(藍色)和通信(紅色)時間占比比較合理,且分配給每個進程的計算任務(wù)基本一致,其中進程間的計算時間占比最大相差不超過8%,有40個進程相差在1%以內(nèi),基本實現(xiàn)了負載均衡,不需要進行進一步優(yōu)化。

圖10 48進程下負載均衡分析Fig.10 Load balancing analysis of 48 processes

為了更全面地對模式的并行能力進行檢驗,我們進一步測試了模式的擴展性。圖11為海洋環(huán)流模式在不同進程數(shù)下的加速比。在單節(jié)點48進程內(nèi),隨著進程數(shù)的增加,加速效率逐漸降低,這是由于增加進程引起的代價不可消除,如進程間通信量增加,帶寬限制等。另外,程序中必定有不能并行的串行部分,即使進程無限增多,通過并行所產(chǎn)生的加速比都是受限的,即阿姆達爾定律[10]。因此模式無法隨著進程的增加獲得線性加速比。此外,當進程數(shù)超過48時,由于測試環(huán)境沒有開啟超線程,進程可能會被掛起等待,同時進程間通信增加,加速效率會繼續(xù)降低。

在經(jīng)過編譯選項優(yōu)化、標量串行優(yōu)化、SIMD優(yōu)化、內(nèi)存帶寬優(yōu)化幾個步驟后,模式整體性能得到了有效提升,如圖12所示。選取合適的編譯選項可以充分發(fā)揮計算機的性能,使得計算效率提升了21%;通過規(guī)范代碼書寫,去除重復計算,減少循環(huán)嵌套等方法,實現(xiàn)串行代碼優(yōu)化,以及通過修改代碼,使之符合編譯器對矢量化的要求,實現(xiàn)編譯器自動矢量化,這兩步優(yōu)化共實現(xiàn)了2%的性能提升;另外,針對模式占用內(nèi)存帶寬過高的問題,提出了去掉臨時數(shù)組的解決方案,使得優(yōu)化后的單個熱點函數(shù)加速比最高可以達到1.4倍,模式整體加速了8%。經(jīng)過上述優(yōu)化,最終模式的整體計算效率共提升了31%,節(jié)省了接近1/3的計算資源。同時,48進程下的MPI并行模式計算效率相較于串行模式提升了15.8倍,雖然單節(jié)點內(nèi)擴展性較弱,但負載均衡性較好。

圖11 不同進程下的加速比Fig.11 Acceleration ratio of different processes

圖12 各優(yōu)化步驟加速比Fig.12 Acceleration ratio of each optimization step

3 試驗結(jié)果分析

3.1 大量指針被假定數(shù)據(jù)依賴

矢量化有許多限制,例如存在數(shù)據(jù)依賴、函數(shù)調(diào)用(數(shù)學庫調(diào)用除外)、在同一循環(huán)中混合可矢量化類型以及存在與數(shù)據(jù)相關(guān)的循環(huán)退出條件等,則不能進行矢量化。其中數(shù)據(jù)依賴分為寫后讀依賴和讀后寫依賴。寫后讀依賴即計算一個數(shù)組元素時需要依賴該數(shù)組已經(jīng)計算的元素,如果強制將其矢量化,可能后面需要讀取的依賴元素并沒有被計算完成,因此造成讀取數(shù)據(jù)錯誤。除了存在不可避免的計算中的前后依賴之外,NEMO的矢量化報告顯示模式源代碼中存在大量的指針,而由于Fortran中的指針可以作為別名指向內(nèi)存地址,因此編譯器通常無法判斷包含指針的代碼是否存在數(shù)據(jù)依賴,從而無法自動矢量化。本文采取了兩種方法來避免因指針造成的無法矢量化問題,包括拆分循環(huán)和將指針修改為動態(tài)數(shù)組。但需指出的是,這只是針對熱點函數(shù)的局部優(yōu)化,因此效果有限。若要全局優(yōu)化此問題,則需要將代碼重構(gòu)。

3.2 循環(huán)嵌套過多

NEMO的核心計算都包含在循環(huán)嵌套中,源代碼中的函數(shù)多是四重循環(huán),包括一維時間循環(huán)以及三維空間循環(huán)。當算例規(guī)模達到一定程度時,會不可避免地出現(xiàn)效率低下。另外,由于編譯器僅會對最內(nèi)層循環(huán)進行自動矢量化,因此循環(huán)嵌套過多也會造成SIMD優(yōu)化效率低下。針對此問題,一種優(yōu)化思路是遵循內(nèi)大外小的原則,即將迭代次數(shù)多的循環(huán)放在最內(nèi)層,提高SIMD優(yōu)化效率。但更有效的解決方案則是減少循環(huán)嵌套層數(shù),但是需要從根本上改變算法問題,這將是非常復雜的工程。

3.3 內(nèi)存帶寬過高

本研究中的串行與標量優(yōu)化以及SIMD優(yōu)化效果有限,而在經(jīng)過內(nèi)存帶寬分析及優(yōu)化后,熱點函數(shù)的性能提升十分明顯,這說明NEMO海洋環(huán)流模式內(nèi)存帶寬占用過高是個比較顯著的問題。這是由于NEMO計算時涉及大量數(shù)組運算,因此會出現(xiàn)大量緩存未命中的情況,需要頻繁從內(nèi)存中讀取數(shù)據(jù)。而當內(nèi)存帶寬占用過高時,處理器從內(nèi)存獲取數(shù)據(jù)的速度會出現(xiàn)瓶頸。去除臨時數(shù)組是一種解決方案,但卻會增加大量的重復計算,增大計算壓力,這與串行代碼優(yōu)化中的去除重復計算優(yōu)化策略恰恰相反。但針對具體的NEMO內(nèi)存帶寬瓶頸問題,去除臨時數(shù)組后減少的時間遠遠大于增加重復計算帶來的計算時間,因此去除臨時數(shù)組是非常有效的。

4 結(jié) 語

在本文中,首先介紹了NEMO模式特征,并根據(jù)代碼現(xiàn)代化指導步驟制定了NEMO的優(yōu)化策略。隨后給出了試驗算例以及測試環(huán)境,并詳細介紹了NEMO代碼現(xiàn)代化的優(yōu)化過程。結(jié)果表明,通過對選取的熱點函數(shù)進行編譯選項優(yōu)化、標量串行優(yōu)化、SIMD優(yōu)化以及內(nèi)存帶寬優(yōu)化,串行的海洋環(huán)流模式的整體加速比可以達到1.31倍,其中熱點函數(shù)加速效果非常明顯,例如traldf_iso相較于原始性能提升了2.3倍。但因模式熱點函數(shù)分散,選取的熱點函數(shù)并不是模式的全部計算瓶頸,因此整體加速比有限。最后對NEMO在MPI多進程下進行了擴展性測試,測試結(jié)果表明在48進程(測試環(huán)境最大核心數(shù))的加速比為15.8倍,計算和通信時間占比合理,負載均衡。但可能是由于內(nèi)存帶寬占用過高,導致單節(jié)點內(nèi)擴展性較差,這需要進行進一步分析驗證。另外,我們還對優(yōu)化后的模式結(jié)果進行了驗證,以確保優(yōu)化后結(jié)果的準確性。需要注意的是,雖然對于海洋模式NEMO來說,我們的優(yōu)化對模擬結(jié)果影響不大,在誤差允許范圍內(nèi),但是在海氣耦合模式中,由于存在強的非線性海氣相互作用,誤差可能會被放大,因此在優(yōu)化時需要注意并檢驗。總的來說,代碼現(xiàn)代化在不增加任何硬件成本的前提下是一種行之有效的優(yōu)化方案。

此外,在對NEMO進行代碼現(xiàn)代化優(yōu)化的過程中,分析出其主要存在以下幾點影響計算效率的問題,如模式大量使用指針、循環(huán)嵌套過多、內(nèi)存帶寬占用過高等。針對這些問題,本文給出了相應(yīng)的解決方案,如將指針修改為動態(tài)數(shù)組,減少循環(huán)嵌套次數(shù)以及去除臨時數(shù)組等,可以為今后海洋模式的設(shè)計和改進提供參考。

另外,MPI并行化僅在二維方向上進行分解(沿著jpi和jpj)[11]。當MPI進程數(shù)增多時,每個進程間的通信量會成倍增加,通信開銷會越來越大,導致模式性能下降。為了減少通信時間,通過引入MPI/Open MP混合并行將會是一個切實可行的方法。Open MP可以實現(xiàn)共享內(nèi)存層面的任務(wù)并行,它是一種簡單有效的并行方式[11-12]。這種混合并行方法通過減少MPI進程數(shù),可以減小進程間通信開銷,同時利用Open MP打開超線程又能充分利用計算資源,使得并行結(jié)構(gòu)更加合理[13]。由于部分NEMO算例的特征是在三維數(shù)組上沿著jpi、jpj和jpk執(zhí)行的操作,因此可以將Open MP共享式內(nèi)存并行化應(yīng)用在jpk(此試驗中為75)方向上,提高進程內(nèi)的計算能力。在熱點函數(shù)中引入MPI/Open MP并行化后,由于進程間通信數(shù)量減少,將有利于提高模式的擴展性。這是接下來研究的方向之一。