李維釗,王偉
易海陸圓(山東)數(shù)字技術(shù)有限公司,山東 青島 266237
數(shù)值預報屬于微分方程的初值問題,初始場的正確性對于預報結(jié)果具有決定性作用。利用數(shù)值模式對大氣、海洋等進行數(shù)值預報,已經(jīng)成為大氣、海洋環(huán)境預測預報的主要方式。數(shù)值預報的資料同化,實質(zhì)上是利用不同來源的觀測資料不斷修正模式的預報,從而迫使數(shù)值預報產(chǎn)生更加接近實況的模式環(huán)境動力狀態(tài)的過程,它為下一時刻的數(shù)值預報提供更加準確的初始場,從而提高數(shù)值預報的準確性。無論在理論研究還是工程實踐中,資料同化已經(jīng)成為數(shù)值預報的關(guān)鍵研究內(nèi)容。
變分法和濾波法是現(xiàn)代資料同化中最有代表性的技術(shù)方法。1959年Cressman在該方法的基礎(chǔ)上采用迭代求解方法形成了實際可以操作的逐步迭代法;1963 年Gardin 提出了最優(yōu)插值法,建立了資料同化的理論基礎(chǔ)——統(tǒng)計估計理論;20世紀60年代后期,隨著觀測資料的時空范圍和種類的增加,出現(xiàn)了變分法的資料同化方法,1969年Thompson等更是提出一套較為完善的四維變分理論和工程方法,但是對于強非線性問題性能較差;1960 年,數(shù)學家Kalman 提出了卡爾曼濾波方法;1994年Evensen等首次提出的集合卡爾曼濾波方法,可以應(yīng)用于強非線性系統(tǒng),經(jīng)過大量的理論研究和應(yīng)用實踐,證明集合卡爾曼濾波以及由此發(fā)展出來的集合最優(yōu)插值具有很好的數(shù)據(jù)資料同化效果和不同場景的適應(yīng)能力,得到廣泛的認可。
集合卡爾曼濾波方法解決了誤差協(xié)方差矩陣預報中遇到的問題,在實際應(yīng)用中,由于大量不同來源、不同應(yīng)用目標、不同空間解析度和對不同物理分量預測預報的需求,導致資料同化的數(shù)據(jù)量和計算量巨大,算法的計算效率成為制約科研和應(yīng)用的一個重要因素?,F(xiàn)代資料同化算法基本上是通過數(shù)值計算的方式迭代求解,特別是矩陣求逆算法,在變分法和集合卡爾曼濾波算法中,成為資料同化的關(guān)鍵技術(shù)之一,也是計算量的主要來源。本文從以下兩個方面提高海量資料同化的效率問題:
(1)面向大規(guī)模計算環(huán)境的并行策略及實現(xiàn)方式;
(2)面向異構(gòu)計算平臺的使用奇異值分解(SVD)的矩陣求逆批量快速算法。
本研究使用C/C++完整地實現(xiàn)了包括SVD在內(nèi)的海洋資料同化軟件。資料同化的實現(xiàn)程序大多采用FORTRAN語言實現(xiàn),而某些國產(chǎn)自主的計算平臺在支持FORTRAN 語言的生態(tài)方面存在不足,與此形成鮮明對照的是工程人員選擇更具備生態(tài)優(yōu)勢的C/C++語言實現(xiàn)程序功能。本研究同時也證明了,只要合理設(shè)計,借助于現(xiàn)代編譯系統(tǒng)的優(yōu)化功能,在同樣的應(yīng)用場景下,C/C++可以達到FORTRAN語言相同的計算效率。
基于文獻[1]的集合卡爾曼濾波及集合最優(yōu)插值,國內(nèi)外出現(xiàn)大量的基于FORTRAN語言和線性代數(shù)庫的算法在高性能計算環(huán)境或計算卡下的實現(xiàn),并進行了大量的優(yōu)化工作,具體見文獻[2-6],盡管如此,資料同化算法仍然是非常耗時的,尤其是對于國產(chǎn)異構(gòu)計算平臺,提升基礎(chǔ)算法的計算效率是應(yīng)用推廣的關(guān)鍵因素之一。
參考文獻[2],集合最優(yōu)插值算法的表述如下:
其中:φa是分析場數(shù)據(jù),φ是背景場數(shù)據(jù),α是區(qū)間( 0,1] 的隨機因子;A是靜態(tài)樣本數(shù)據(jù),A′是擾動預測矩陣,由A中心化而來,A′A′T是背景誤差協(xié)方差的估計;H是測量算子,代表真實模型與預測模型之間的映射;γ是擾動測量誤差向量,γγT是觀測誤差協(xié)方差;d′=d-Hφ是觀測增益,其中d是擾動測量向量。經(jīng)過一系列的推導,參考文獻[2],可以得到如下的公式:
其中U、Λ 滿足Λ=ΣΣT,對進行SVD分解,得到U向量和奇異值Σ,結(jié)合矩陣的性質(zhì)等,得到如下步驟計算:
靜態(tài)樣本的選取參考文獻[7],采用文獻[8]中提到的局部集合最優(yōu)插值方法進行資料同化。本文基于國產(chǎn)異構(gòu)計算平臺,采用CPU和類GPU卡協(xié)同方式實現(xiàn)海量資料同化的并行計算。
以資料同化為例,一般的大規(guī)模計算應(yīng)用需要完成如下幾項工作:
2.1.1 數(shù)據(jù)讀取
從存儲系統(tǒng)將數(shù)據(jù)讀入內(nèi)存。主要的數(shù)據(jù)包括靜態(tài)樣本數(shù)據(jù)、背景場數(shù)據(jù)和觀測場數(shù)據(jù),本研究中靜態(tài)樣本數(shù)據(jù)量在數(shù)百GB,通過數(shù)個IO代理進程從系統(tǒng)的存儲設(shè)備中將數(shù)據(jù)讀入CPU內(nèi)存,再將數(shù)據(jù)分發(fā)給資料同化的計算進程。數(shù)據(jù)在文件中的數(shù)據(jù)結(jié)構(gòu)對IO的性能有顯著的影響,原則上,應(yīng)該優(yōu)化數(shù)據(jù)在文件中的存儲方式,使得打開/關(guān)閉文件的操作數(shù)量盡量少,并且文件操作的并發(fā)數(shù)量要與存儲系統(tǒng)的配置相匹配。
2.1.2 數(shù)據(jù)分發(fā)
采用并行計算的方式,資料同化軟件需要將靜態(tài)樣本、背景場和觀測場的數(shù)據(jù)通過通信網(wǎng)絡(luò)分發(fā)到計算進程。
2.1.3 資料同化計算
資料同化計算部分主要包括三個部分:SVD之前的數(shù)據(jù)插值等前處理;SVD計算任務(wù);SVD 后處理。SVD 前、后的處理主要為矩陣的乘法運算。同時,SVD 的主要計算量來源是矩陣的二對角化。單純用CPU 或計算卡進行SVD,非常耗時,因此,在異構(gòu)計算平臺上,發(fā)揮CPU 處理存在復雜邏輯的計算,結(jié)合計算卡的數(shù)值計算性能,實現(xiàn)高效SVD的工程實現(xiàn)算法,是提高資料同化的計算效率的關(guān)鍵技術(shù)。SVD算法在量子計算、人工智能、機器視覺、信號降噪等領(lǐng)域應(yīng)用廣泛,實現(xiàn)高效SVD 的工程實現(xiàn)算法具有重要的實用價值。
2.1.4 數(shù)據(jù)回寫
計算進程計算完成后,需要將結(jié)果數(shù)據(jù)回寫到分析場數(shù)據(jù)文件中。
已有大量的文獻對資料同化工程實現(xiàn)的并行策略進行了研究,并取得不同程度的性能提升。與文獻[3]和文獻[2]相似,本研究資料同化的并行策略可以總結(jié)為計算任務(wù)和流水線并行策略,其核心為任務(wù)的編排。
計算任務(wù)的并行度,取決于任務(wù)粒度、通信方式和負載平衡三個方面因素。因此,資料同化任務(wù)的分解和調(diào)度是選擇并行策略的基礎(chǔ)。
海洋資料同化的一般流程如圖1所示。
圖1 海洋資料同化流程圖Fig.1 Ocean data assimilation flowchart
按照海洋資料同化流程,結(jié)合國產(chǎn)異構(gòu)計算平臺的環(huán)境配置,實現(xiàn)并行化設(shè)計,海洋資料同化的任務(wù)模塊劃分如下:
2.2.1 數(shù)據(jù)訪問任務(wù)模塊
同化軟件啟動Nio個IO 代理進程實現(xiàn)海洋資料同化的數(shù)據(jù)文件的輸入和輸出。數(shù)據(jù)訪問的效率取決于存儲系統(tǒng)的并行特性和文件的數(shù)據(jù)結(jié)構(gòu)。其中文件的數(shù)據(jù)結(jié)構(gòu)是關(guān)鍵。在HYCOM 數(shù)值模式生成的數(shù)據(jù)中,共有120 個包括垂向36層網(wǎng)格數(shù)據(jù)的文件。資料同化程序的代理進程讀取相應(yīng)的數(shù)據(jù),并將數(shù)據(jù)按照計算進程的數(shù)量及序號,平均分配到相應(yīng)的計算進程執(zhí)行同化計算。原則上,120份垂直36層的靜態(tài)樣本數(shù)據(jù)可以組織到1或多個文件中,工程實現(xiàn)時,文件的數(shù)量NF和每個文件中的數(shù)據(jù)結(jié)構(gòu)對數(shù)據(jù)IO 的效率影響明顯。以溫度樣本為例,原程序文件組織形式為NF=120 個文件,每個文件中的溫度的數(shù)據(jù)結(jié)構(gòu)如圖2所示。
圖2 靜態(tài)樣本數(shù)據(jù)原始組織方式示意圖Fig.2 Schematic diagram of the original organization of static sample data
由于每個計算進程分配的數(shù)據(jù)為36層數(shù)據(jù)中相同區(qū)域的網(wǎng)格點數(shù)據(jù),如圖3所示(Np為計算進程的數(shù)量)。
圖3 靜態(tài)樣本數(shù)據(jù)在計算進程中的分配示意圖Fig.3 Schematic diagram of static sample data allocation in the calculation process
其中,總的文件數(shù)量為NF=120,每個文件中都有圖3 所示的一份數(shù)據(jù)。每層數(shù)據(jù)的讀取均遵循“打開文件—讀取數(shù)據(jù)—關(guān)閉文件”的文件讀取操作,每個進程執(zhí)行文件讀取操作的數(shù)量K=NF×36,文件讀取操作的總數(shù)量為KTotal=Nio×NF×36,當Nio=8,NF=120 時,KTotal=34,560 次文件操作。由于存在多個代理進程讀取數(shù)據(jù),120 個文件的讀寫存在大量的沖突,超出存儲系統(tǒng)的并發(fā)性能,造成文件讀取效率低下。數(shù)據(jù)IO 優(yōu)化的目標有兩個:一是通過優(yōu)化數(shù)據(jù)結(jié)構(gòu),減少文件操作的次數(shù),增加單次文件連續(xù)讀操作的數(shù)據(jù)量;二是將文件合理放置在存儲系統(tǒng)不同物理存儲單元上,實現(xiàn)數(shù)據(jù)訪問的高并發(fā)性,如文獻[2],在基于Lustre 的文件系統(tǒng)上,通過將文件放置在不同OST上,取得明顯的性能優(yōu)化。本文考慮到存儲文件系統(tǒng)的多樣性和復雜性,通過優(yōu)化文件的數(shù)據(jù)結(jié)構(gòu)實現(xiàn)IO性能的提升。根據(jù)局部集合最優(yōu)插值算法,網(wǎng)格點的同化計算順序不影響最終結(jié)果,因此,設(shè)定1個最大數(shù)量的IO進程數(shù)Nio-max=32,在資料同化進行之前,將靜態(tài)樣本數(shù)據(jù)的文件存儲格式按照如圖4組織(文件數(shù)量NF=120)。
圖4 靜態(tài)樣本數(shù)據(jù)優(yōu)化組織方式示意圖Fig.4 Schematic diagram of static sample data optimization organization method
由圖4可知,數(shù)據(jù)優(yōu)化的方式就是將原來按照層序號連續(xù)放置的一維方式,改變?yōu)榘凑誟層、進程]的二維方式,適應(yīng)計算進程對數(shù)據(jù)結(jié)構(gòu)要求,遵循同樣的文件操作方式,文件讀操作的最大數(shù)量變?yōu)椋∟io=8 時)KTotal=Nio×NF=8×120=960,不同的文件代理進程讀取不同的文件,結(jié)合數(shù)據(jù)分發(fā)任務(wù)模塊的具體特點,避免了進程間文件訪問的沖突,實驗證明,優(yōu)化措施帶來顯著的文件性能IO 提升,并且具有一定的普遍適用性。表1 是Nio=4,8,12,Nf=120 時,兩種讀寫方式的運行時間對比。優(yōu)化后運行效率提高到原來的10倍(Nio=8 時性能最優(yōu),不同存儲系統(tǒng)存在差異)。
表1 不同讀寫策略讀寫效率對比Table 1 Comparison of reading and writing efficiency of different reading and writing strategies
2.2.2 數(shù)據(jù)分發(fā)任務(wù)模塊
通過文件讀取操作的優(yōu)化,文件中的數(shù)據(jù)結(jié)構(gòu)得到優(yōu)化,IO代理進程在讀取數(shù)據(jù)后,不需要進一步組織數(shù)據(jù),就可以把讀取的數(shù)據(jù)直接分發(fā)給計算進程,計算進程接收到數(shù)據(jù)后,啟動執(zhí)行同化計算的任務(wù)模塊,按照網(wǎng)格點實施同化計算,文獻[2]中對采用MPI 相關(guān)通信機制進行了詳細說明。研究證明,通過對文件數(shù)據(jù)結(jié)構(gòu)的優(yōu)化,結(jié)合MPI的通信機制,充分利用分布式計算中心通信網(wǎng)絡(luò)的特性,可以實現(xiàn)高效的數(shù)據(jù)讀取和分發(fā)。這不是本研究的重點,請參考文獻[2]。
2.2.3 計算任務(wù)模塊
采用局部集合最優(yōu)插值的資料同化算法,HYCOM 數(shù)值模式中,地球按照經(jīng)緯度劃分為4,500×3,200 的網(wǎng)格,對周圍一定范圍內(nèi)的觀測數(shù)據(jù),按照網(wǎng)格點逐一進行同化計算,即公式(2)和公式(3)-(7)。因此,資料同化具有良好的計算并行性,將每個網(wǎng)格點的計算作為一個基本的計算任務(wù),劃分為M=4,500×3,200=14,400,000 個基本計算單元,連同數(shù)據(jù),將M 個計算單元按照一定的規(guī)則(如文獻[2]的負載均衡)分配到Np個計算進程上,由Np個計算進程并行完成。文獻[2]中對于全球網(wǎng)格的非海洋網(wǎng)格點進行剔除后,重新進行計算負載的分配,計算效率獲得顯著提升。文獻[2]中通過計算負載均衡整個同化效率得到提升,而有效計算核時并沒有改變,同化計算本身的效率并沒有提高,本文第3部分,提出了一種適合在國產(chǎn)異構(gòu)計算平臺顯著提升同化計算效率(核心是SVD 算法的實現(xiàn))的方法。
整個資料同化軟件的操作可以分為:數(shù)據(jù)讀入、數(shù)據(jù)分發(fā)、同化計算和數(shù)據(jù)回寫四個部分。四個部分順序執(zhí)行,而采用流水線的設(shè)計思路,在數(shù)據(jù)讀入運行開始后,當數(shù)據(jù)讀入到一定數(shù)量,啟動數(shù)據(jù)分發(fā),而計算進程在接收到必要的數(shù)據(jù)后,啟動同化計算,當計算進程完成計算后,一次性將結(jié)果數(shù)據(jù)回寫到代理進程,由代理進程將最終的數(shù)據(jù)寫入文件。這樣實現(xiàn)四個部分的操作在時間上的重疊,如圖5所示。
圖5 任務(wù)模塊流水線操作示意圖Fig.5 Schematic diagram of task module pipeline operation
本文的第2部分和文獻[2]中,通過對資料同化的任務(wù)編排、流水線操作和負載均衡進行優(yōu)化,可以有效提高包括資料同化在內(nèi)的并行計算應(yīng)用的計算效率,而資料同化的基本計算單元包括SVD算法在內(nèi)的矩陣運算,計算量大,現(xiàn)有算法并行化程度低,對其進行優(yōu)化可以從基礎(chǔ)算法方面提高計算速度。
SVD 算法通常采用基于QR 迭代和Jacobi旋轉(zhuǎn)的方法實現(xiàn),文獻[4]、[5]分別對基于Jacobi旋轉(zhuǎn)和QR迭代的方法,單獨使用計算卡或通過CPU-GPU 協(xié)同并行計算實現(xiàn)SVD 算法做了研究。本文通過研究QR迭代的SVD求解算法,提出CPU-GPU 協(xié)同并行計算實現(xiàn)SVD 算法的詳細步驟和驗證結(jié)果。
對于批量矩陣的SVD算法,優(yōu)化策略如下:
(1)根據(jù)計算卡雙精度浮點處理器內(nèi)核的數(shù)量、計算卡內(nèi)存數(shù)量和矩陣的大小,確定SVD算法每個批次的矩陣數(shù)量Nm,該數(shù)量通過對計算卡的測試獲得。
(2)GPU承擔矩陣的二對角化計算任務(wù),并負責Givens變換系數(shù)計算奇異向量U、V矩陣的任務(wù)。
(3)CPU 承擔QR 迭代,保存每一步Givens變換的系數(shù),用于計算卡計算U、V矩陣。
(4)合理規(guī)劃和使用計算卡內(nèi)存和主機內(nèi)存,設(shè)計CPU和計算卡的計算任務(wù)調(diào)度機制,實現(xiàn)高效的調(diào)度算法。
基于QR 迭代的方法是求解二對角矩陣奇異值的最快速的算法。矩陣的SVD 求解,可以分為三步:第一步是通過Householder變換,將矩陣轉(zhuǎn)換為二對角矩陣;第二步采用QR 迭代,得到矩陣的奇異值,保存Givens變換的系數(shù);第三步使用Householder 變換時得到U、V,然后再使用QR迭代過程中Givens變換的系數(shù),迭代計算奇異向量U、V矩陣。
3.2.1 SVD算法原理
SVD的定義如下:
一個m×n矩陣,定義A=UΣVT,其中U是一個m×m矩陣,U是一個n×n矩陣,Σ 是一個m×n矩陣,其除主對角線上元素外,其他元素為0,并且主對角線上元素σi>0。
在資料同化應(yīng)用中,由于m=n,即待分解的矩陣為方陣,且為對稱方陣,因此,U=V。在資料同化的SVD操作,見公式(1)-公式(7)。
3.2.2 SVD解算的實現(xiàn)流程和數(shù)據(jù)結(jié)構(gòu)
從QR 迭代法求解矩陣的奇異值和奇異向量的一般算法出發(fā),根據(jù)國產(chǎn)異構(gòu)計算平臺的特點,本文設(shè)計了SVD 解算的實現(xiàn)流程和數(shù)據(jù)結(jié)構(gòu),具體流程如圖6。
圖6 異構(gòu)計算平臺SVD算法流程Fig.6 SVD algorithm process for heterogeneous computing platforms
說明:
(1)CPU 程序負責申請主機內(nèi)存和計算卡的存儲,CPU內(nèi)存使用鎖頁內(nèi)存,從而提高了內(nèi)存訪問的速度。
(2)從流程圖中可以看到,在SVD 解算過程中,CPU 和計算卡均需要等待,設(shè)CPU 執(zhí)行時間為tcpu,計算卡執(zhí)行時間為tgpu,4 次數(shù)據(jù)交換時間為tdata,1 次批量SVD 的矩陣數(shù)量為n,所以執(zhí)行1 個矩陣SVD 總的運行時間為t=
(3)在SVD流程中,CPU和GPU均存在等待時間,且等待時間基本上等于對方的執(zhí)行時間。
為實現(xiàn)高效的SVD 解算,在主機鎖頁內(nèi)存中開辟2份緩沖區(qū),分別存放n個矩陣的數(shù)據(jù)及計算結(jié)果數(shù)據(jù),按照如圖7 的方式執(zhí)行SVD 算法流程。
通過合理配置批量SVD 矩陣的數(shù)量,使得CPU 執(zhí)行時間、GPU 執(zhí)行時間和數(shù)據(jù)交換時間達到平衡,可以顯著提升批量SVD的計算效率。
3.2.3 矩陣二對角化(Householder變換)
Householder 變換最初由A.C Aitken 在1932年提出,是一種線性變換,其變換矩陣被稱作Householder矩陣,其定義具體如下:
設(shè)ω∈Rn,‖ω‖2=1,定義H=I-2ωωT(H∈Rn×n)。
Householder變換可以將矩陣的某些元素置零而矩陣的范數(shù)不變,正是利用此性質(zhì),將待分解的矩陣先變?yōu)槎蔷仃嚕粗鲗蔷€和上對角線之外的元素都為零。如下所示3×3 二對角矩陣:
二對角變換算法如下:
(1)設(shè)矩陣An×n是n 行n 列的輸入矩陣。如下:
(2)x=[a11,a21,…,an1],計算v=x±‖x‖2e1,其中e1為除第一個元素為1,其他元素是0的單位向量。
(3)計算ω=v/‖v‖2。
(4)計算H=I-2ωωT。
(5)計算A′n×n=HAn×n,A′n×n為:
(6)x=[a′12,a′13,…,a′1n],計算v=x±‖x‖2e1。
(7)計算ω=v/‖v‖2。
(8)計算H=I-2ωωT。
(9)計算A″n×n=HA′n×n,A″n×n為:
(10)將A″n×n從a″22開始,組成新的矩陣An-1×n-1,重復(2)-(9),直到矩陣的大小變?yōu)?×2。
由于二對角采用計算卡計算,并且二對角過程中的ω組成U 和V 向量矩陣需要保存,按照如下原則進行計算,提升計算速度:
(1)利用矩陣運算的性質(zhì)減少乘法次數(shù)。具體如下:
并且在計算ω=v/‖v‖2時從而有:
可以明顯減少乘法次數(shù),提升計算速度。
(2)批量計算。以‖v‖2的計算為例,具體如下:
如果按照常規(guī)計算,設(shè)計一個循環(huán),對v的每個元素進行平方后再求和,每次循環(huán)做一次乘法和一次累加,這樣的計算并沒有充分發(fā)揮GPU 的內(nèi)核緩存性能;優(yōu)化的算法是每個循環(huán)進行k次的乘法和k次的累加,k一般取4,8,16。對于某國產(chǎn)異構(gòu)計算平臺的計算卡,當k=8 時,計算速度最快。在二對角化過程中,對類似的運算(主要是矩陣乘法等運算)進行如是處理,可以明顯提升計算速度。
3.2.4 QR迭代
二對角化計算完后,得到二對角矩陣Bn×n和正交向量Un×n、Vn×n,由于二對角化是在異構(gòu)計算平臺的計算卡上運行的,在進行QR 迭代前,要將數(shù)據(jù)搬移到CPU的內(nèi)存中,見圖8,而正交向量Un×n、Vn×n將在QR 迭代后,使用保存的Givens變換的系數(shù)進行計算,得到奇異向量。
圖8 Givens變換矩陣元素的數(shù)據(jù)結(jié)構(gòu)Fig.8 The data structure of givens transformation matrix elements
二對角矩陣Bn×n如下所示(以n=3為例):
對Bn×n沿主對角線和上對角線元素順次進行Givens 變換(具體見相關(guān)文獻),使得次對角線的元素為零,并保存Givens 變換的系數(shù)c和s,不斷迭代,直到將次對角線元素全部變?yōu)榱悖▽嶋H計算時設(shè)置一個閾值,當小于該閾值時就停止迭代)。
由于QR 迭代存在大量的分支,并且QR 迭代是串行進行的,因此,不適合使用計算卡并行實現(xiàn),該部分計算由更適合進行邏輯判斷的CPU 進行。計算完成后,對得到矩陣的主對角線元素進行如下處理:
(1)如果主對角線元素大于零,不變;如果小于零,取絕對值,同時,將對應(yīng)的U和V矩陣對應(yīng)向量的元素符號取反。
(2)對主對角線按照由大到小排序,同時調(diào)整U和V矩陣對應(yīng)向量的順序,以保持一致。
以上兩步操作在完成奇異向量計算后進行。
3.2.5 奇異向量的計算
QR迭代完成后,將保存的Givens變換的系數(shù)搬移到計算卡內(nèi)存,對二對角化后得到的U和V,與Givens變換矩陣相乘,最終得到奇異向量U和V。由于QR 迭代產(chǎn)生的Givens 變換矩陣必須按照相同的順序作用到二對角化后得到的U和V,設(shè)計如圖8數(shù)據(jù)結(jié)構(gòu):
其中,Nj:第j次迭代的矩陣數(shù)量;SR和ER:進行矩陣運算時,Givens變換矩陣非零部分是2×2 的大小,與U和V對應(yīng)的行列的子矩陣進行運算,SR和ER分別是起始行和結(jié)束行的位置。
在計算卡中計算奇異矩陣時,通過遍歷該數(shù)據(jù)結(jié)構(gòu),完成奇異向量的高效計算。
包括資料同化在內(nèi)的大多數(shù)的科學計算算法均由Fortran 語言實現(xiàn),F(xiàn)ortran 的語言特性是接近數(shù)學公式描述,具有很高的執(zhí)行效率,可以直接對矩陣和復數(shù)進行運算。盡管Fortran具備上述優(yōu)點,但是相對于C/C++語言來說,開發(fā)環(huán)境、運行環(huán)境和生態(tài)方面,差距巨大,尤其是對于非科學計算領(lǐng)域的工程人員,不掌握Fortran編程開發(fā)能力。加之某些異構(gòu)計算平臺不支持Fortran開發(fā)環(huán)境,而均支持C/C++開發(fā)和運行環(huán)境。
鑒于上述原因,本文所涉及的資料同化軟件,全部使用C/C++開發(fā),對于SVD解算核心算法,使用C語言開發(fā),并采用模塊化設(shè)計思路,矩陣乘法、矩陣二對角化、QR分解、向量計算等均是單獨的模塊,封裝成動態(tài)或靜態(tài)鏈接庫,可以通過設(shè)置不同的參數(shù)靈活調(diào)用。經(jīng)實際驗證,采用C語言設(shè)計實現(xiàn)的SVD算法在相同的CPU上的運行速度,與LAPACK庫相比,運行速度相當,而對于采用CPU和計算卡協(xié)同計算相比,具有顯著的性能提升,具體見3.4實測結(jié)果。
整個程序采用C/C++設(shè)計,可以運行在Linux操作系統(tǒng)和MPI并行環(huán)境下,也可以運行在Windows 操作系統(tǒng)下,Windows 版本提供方便的圖形化操作界面,便于配置和交互。
實際測試的環(huán)境為某國產(chǎn)異構(gòu)計算平臺,每個計算節(jié)點1 顆CPU,2 張國產(chǎn)加速卡。分為SVD核心算法測試和整個同化軟件的測試。
類GPU環(huán)境與資源配置如下:
(1)CPU 節(jié)點配置1 塊32 核心的CPU,主頻2.0 GHz 處理器雙精度浮點理論值0.256 TFLOPS。
(2)異構(gòu)節(jié)點,在CPU節(jié)點基礎(chǔ)上,配置4塊16G 內(nèi)存的類GPU 加速卡,單卡雙精度浮點理論值6.96TFLOPS。
3.4.1 SVD核心算法測試
測試方法分為兩種,一種是采用隨機數(shù)測試,具體如下:
(1)對比測試的庫為rocm的dgesvd_batched函數(shù)。
(2)采用8,192個100階方陣,矩陣元素為隨機數(shù)。
測試結(jié)果如表2。
表2 SVD核心算法對比測試數(shù)據(jù)(方法1)Table 2 SVD core algorithm comparison test data(method 1)
從實際測試數(shù)據(jù)看,自主開發(fā)的批量SVD算法的計算速度是ROCM 庫的9.51 倍,計算速度顯著提升。
另一種測試方法是在同化程序優(yōu)化前后,對SVD 解算部分的單獨進行計時,進而測試SVD在國產(chǎn)異構(gòu)計算平臺同化程序中的性能提升,由于采用并行計算,為確保測算準確性,抽取優(yōu)化前后運行時解算矩陣數(shù)量最多的進程的運行時間。具體如表3。
表3 SVD核心算法對比測試數(shù)據(jù)(方法2)Table 3 SVD core algorithm comparison test data(method 2)
同化程序中矩陣的大小均為100×100,雖然測試中矩陣元素的數(shù)據(jù)不同,由于矩陣數(shù)量較多,因此,測試數(shù)據(jù)的有效性可以保證。從測試數(shù)據(jù)可知,計算效率提升到23.33 倍,提升效果顯著。
方法1 和方法2 效率的不同主要原因是兩種測試的數(shù)據(jù)不同,方法1 的數(shù)據(jù)為隨機數(shù)據(jù),迭代次數(shù)多,而方法2 的數(shù)據(jù)在進行二對角化前,主對角線的元素絕對值遠遠大于其他元素,導致對角化后QR的迭代次數(shù)較少。
3.4.2 同化軟件測試
整個同化軟件在采用上述的各項優(yōu)化措施后,經(jīng)測試計算性能得到顯著提升。
運行環(huán)境如下:
硬件測試環(huán)境:CPU為20顆21-core Processor,本地硬盤為240GB M.2 Nvme SSD,內(nèi)存為8*16GB DDR4 2666MHz,存儲為Parastor300S。
軟件測試環(huán)境:操作系統(tǒng):CentOS7.6,國產(chǎn)異構(gòu)加速卡環(huán)境:DTK-21.10.1庫,MPI:intelmpi/2020.1.217,編譯器:mpiifort,hipcc。
測試算例如下:
基于國產(chǎn)異構(gòu)計算架構(gòu)的全球高分辨率海洋資料同化高效軟件,共分為三個算例,分別是對三個不同時間的海表溫度的海洋資料同化。算例1、2、3分別是對2020年9月29日、2013年9月22日、2012年9月9日的海表溫度進行資料同化。每個算例包括優(yōu)化前程序與優(yōu)化后程序,優(yōu)化前程序是原始的同化程序,優(yōu)化后程序是使用自研batchedsvd 算法以及多種策略優(yōu)化的同化程序。
具體測試數(shù)據(jù)采用3 個算例平均,具體如表4。
表4 海洋資料同化優(yōu)化前后性能比較Table 4 Ocean data assimilation performance comparison
從測試數(shù)據(jù)看,程序運行速度提高到原來的2.5 倍以上,IO 運行速度提升到20 倍左右,計算部分的加速比顯著提升。
本文所涉及的同化軟件及SVD算法方面在大規(guī)模計算應(yīng)用方面的貢獻如下:
(1)使用C/C++完整實現(xiàn)同化軟件,同化軟件便于可以集成或被集成到其他軟件系統(tǒng),對于豐富大規(guī)模計算尤其是國產(chǎn)異構(gòu)平臺的生態(tài)體系、促進同化軟件的發(fā)展意義重大。
(2)使用C 語言基于國產(chǎn)異構(gòu)計算平臺,完整實現(xiàn)了SVD 算法,包括矩陣乘法、轉(zhuǎn)置、Householder 變換、Givens 變換、矩陣二對角化、QR迭代,在批量矩陣SVD速度方面取得顯著提升,對于科學計算、人工智能、信號降噪、視頻推薦等領(lǐng)域的工程應(yīng)用意義重大,具體見文獻[9-11]。
(3)以資料同化應(yīng)用為例,對于數(shù)據(jù)和計算雙密集型的應(yīng)用在國產(chǎn)異構(gòu)平臺的設(shè)計思路、優(yōu)化方法做了一般性的探索,并形成了較為成熟的應(yīng)用開發(fā)經(jīng)驗和指導原則,對于同類型的大規(guī)模計算應(yīng)用開發(fā)具有一定的借鑒意義。
下一步的工作:
(1)同化軟件中,除SVD算法外,插值、距離計算等計算任務(wù)的計算量比較大,并且適合使用計算卡,配合篩選算法,實現(xiàn)CPU和計算卡的協(xié)同計算方式,進一步提升計算速度。
(2)開展以現(xiàn)有SVD解算為基礎(chǔ),面向國產(chǎn)異構(gòu)計算平臺,使用C/C++語言,開發(fā)和豐富跨平臺的線性代數(shù)的相關(guān)算法庫。
本文主要針對國產(chǎn)異構(gòu)計算平臺環(huán)境下的數(shù)據(jù)和計算雙密集型的同化算法的高效實現(xiàn)給出相關(guān)策略、方法和創(chuàng)新,可對類似的應(yīng)用提供借鑒,特別是采用C/C++實現(xiàn)的以SVD 解算為核心的自主開發(fā)線性代數(shù)基本計算庫,豐富了國產(chǎn)異構(gòu)計算平臺的通用計算卡底層生態(tài),擴展了國產(chǎn)異構(gòu)計算平臺的應(yīng)用領(lǐng)域。
利益沖突聲明
所有作者聲明不存在利益沖突關(guān)系。