MPI+CUDA聯合加速重力場反演的并行算法

趙鍇坤 朱炬波 谷德峰 韋春博

1 中山大學物理與天文學院天琴中心,廣東省珠海市大學路2號,519082

高精度、高分辨率的重力場模型在氣候變化、資源勘探、災害預防、國防安全等領域有著非常重要的作用。重力場解算的難度主要來源于龐大的計算量與巨大的內存占用。法方程的構建是重力場反演過程中的關鍵步驟,而法矩陣的規模與重力場階數相關[1],是直接影響反演速度與內存占用的重要因素。當重力場階數過大時,法方程的求解速度會大幅下降,故高階重力場的解算尤為困難[2]。隨著計算機技術的發展,已有一些新的技術與高性能算法可用于重力場的快速反演。鄒賢才等[3]使用OpenMP對重力場解算過程進行并行化處理,計算速度提升近6倍;陳秋杰等[4]使用MKL(math kernel library)科學計算庫與OpenMP并行計算庫,顯著提高了重力場解算的效率;周浩等[5-6]基于MPI并行技術對法方程構建、大規模矩陣相乘、求逆等關鍵節點引入并行計算,解算120階重力場模型僅耗時40 min。

除CPU平臺的并行計算方法與高性能算法外,近年來計算機圖形學也得到快速發展,圖形處理器(GPU)的計算能力甚至已經超過CPU。GPU編程工具的進步也使得各種算法成功移植到GPU上,混合計算平臺已經實現了針對不同學科應用的大幅提速[7]。Hou等[8]提出一種結合CUDA和OpenMP的混合并行方法,能處理大規模的數據;劉世芳等[9]實現了基于MPI和CUDA的稠密對稱矩陣三角化算法,取得了較好的加速效果。

上述研究均顯示出并行計算在重力場相關研究中能起到顯著作用。為提高使用最小二乘法反演重力場模型的解算速度,本文針對復雜的重力場反演流程,聯合MPI與CUDA設計合理的并行方案,同時實現內存耗用與計算速度等多方面的優化。

1 基于動力積分法的重力場模型解算方法

動力積分法是較早應用于重力場反演的基礎方法之一,具有很高的求解精度,其計算流程可以概括為:

1)使用先驗力學模型計算各種攝動加速度對待估參數的偏導數,并通過積分器獲得位置對各待估參數的偏導,構建設計矩陣;

2)通過積分器獲得參考軌道序列,并與實測數據作差,獲得O-C殘差序列b;

3)使用最小二乘法求解方程Ax=b,獲得衛星初始狀態改進量ΔX和重力場位系數的改進量Δp;

4)對初始狀態和先驗重力場模型進行改進,重復上述流程,直到殘差滿足收斂條件。

下面將詳細闡述法方程的構建與求解過程,以便對線性的重力場反演任務進行分配與規劃。

首先,將設計矩陣A按照參數類型進行分塊,表示為:

A=[LG]

(1)

式中,L為低軌衛星位置對局部參數的偏導數,局部參數是指如衛星初始狀態向量、加速度計測定的非保守力等與時間有關的參數;G為低軌衛星位置對全局參數的偏導數,全局參數是指如重力場位系數等與時間無關的參數。通過該方式對A進行分塊,有利于多弧段法方程的合并。

將分塊后的單歷元觀測方程構建為法方程:

ATAx=ATy

(2)

(3)

式中,ΔX為局部參數的改正量,Δp為全局參數的改正量。上式簡記為:

(4)

由于單個弧段的所有歷元的待估參數完全一致,故單弧段的法方程形式與式(4)一致。但多弧段間的局部參數存在區別,需采用式(5)實現多弧段法方程合并:

(5)

使用最小二乘法將式(5)整理成法方程形式并化簡為:

(6)

消去局部參數即可解得重力場位系數改進量:

(7)

將式(7)代入式(6)即可求得衛星初始狀態向量的改進量:

(8)

動力法反演重力場模型的詳細流程見圖1。

max(ΔX)表示衛星狀態修正量在各個方向上的最大值,下同

經測試,解算過程中耗時最多的環節為數值積分器計算參考軌道、大規模矩陣相乘以及大規模矩陣求逆。除計算速度外,由于設計矩陣大小近似為3倍歷元個數乘以n位系數,而n位系數又與重力場模型階數n相關(式(9)),故當重力場模型階數n過高時,設計矩陣的大小會劇增,需要兼顧計算過程中的內存及顯存耗用,故有必要根據不同需求分別制定不同的并行加速策略:

n位系數=(n+3)(n-1)

(9)

2 針對復雜計算邏輯的MPI并行策略

MPI提供的各接口可以使不同進程間的信息相互傳遞[5],進而實現復雜任務的并行化處理。使用MPI實現按弧段并行構建設計矩陣及求解法方程,具體任務分配及運行流程見圖2。

圖2 MPI任務分配及計算流程

在MPI中,各進程編號通常為從0開始的整數。為便于表述,稱0號進程為主進程,其他進程為分進程。由于軌道數據、加速度計數據、姿態數據等均按日期存儲,故需要使用主進程讀取各軌道文件,按弧段整理完畢后使用MPI_Send命令將數據發送到各分進程,再啟動積分器執行后續任務。針對內存有限的單機計算,則可以按天數分配任務,為配合部分對數據連貫性有要求的函數(如批量讀取加速度計數據的函數),直接由各分進程讀取連續數日的數據文件,并按弧段進行存儲。為使各分進程的任務均勻分配,分進程總數可定為任務總數的因數。

在計算法方程前,式(1)中低軌衛星位置對重力場位系數的偏導數矩陣G的內存耗用最大。以單弧段3 h共2 160個歷元為例,不同重力場模型階數n對應的矩陣G的規模見表1。

表1 不同階數n對應的矩陣G的大小

一般而言,集群的單個節點內存在16 GB以上,故在解算法方程前可以保證內存充足,無需通過分塊處理、寫文件存儲等方式防止內存溢出。

在使用最小二乘法求解時,需要通過式(2)構建法方程,單弧段法方程的構建結果見式(4)。當位系數個數充足時,式(4)中的Ngg大小將超過矩陣G。不同重力場模型階數n對應的Ngg大小見表2。

表2 不同階數n對應的矩陣Ngg的大小

顯然,當重力場模型階數達到一定程度時,Ngg占用的內存會過大,需要使用合適的方法來控制其大小。

3 矩陣相乘的并行加速方法

矩陣相乘是重力場解算流程中耗時最多、內存占用最大的環節。針對大規模矩陣相乘,可利用CUDA將程序編譯為核函數(kernel)并發送至設備(device),進而在GPU上實現矩陣相乘的并行運算。CUDA是由NVIDIA推出的運算平臺和編程模型,包含有CUDA指令集架構(ISA)與GPU內部的并行計算引擎,使開發者編寫的程序能在GPU上高速運行。另外,當重力場階數過大時,反演流程中的各矩陣,尤其是Ngg會占用大量顯存,而顯存容量一般較小,會嚴重限制可解算重力場階數的上限。本節將利用CUDA和MPI實現大規模矩陣相乘的并行加速與內存壓縮,同時解決耗時過長、內存與顯存占用過大的問題。

法方程構建過程中涉及多個矩陣,其中以Ngg規模最大、占用內存最多,因此本文主要針對Ngg矩陣進行研究。

3.1 使用CUDA實現并行加速

矩陣G的大小為3倍歷元個數乘以n位系數,是設計矩陣A的主要成分。由式(9)可以推算出串行計算Ngg=GTG的乘法次數約為3×歷元個數×n4。而使用CUDA并行計算時,能使用大量線程(thread)同時計算,每個線程負責矩陣中單個元素的計算,使乘法次數下降到3×歷元個數,大幅減少了計算次數。

整個過程在設備上進行,計算完畢后再將結果轉移到主機。在GPU上計算的函數稱為核函數,其邏輯見圖3。設線程的id為idx,則idx號線程負責的元素位置為(x/列數,y%列數),%為取模運算符;再將GT的第x行與G的第y列對應元素相乘后求和,即可求得Ngg的(x,y)號元素的值。

圖3 CUDA計算矩陣相乘的簡單邏輯

針對線性代數運算,CUDA提供了專門的子程序庫cuBLAS(CUDA basic linear algebra subroutine library),其level-3函數cublasDgemm適用于double型常規矩陣相乘運算,比未經處理的簡單矩陣并行相乘的核函數計算更快,后文的實驗也將使用cuBLAS中的函數進行計算。

針對前述內存占用過大的問題,可根據Ngg的特殊性進行優化,主要有以下2點:

1)Ngg具有對稱性,對于C=AAT型的矩陣相乘,cuBLAS提供了cublasDsyrk函數進行運算,并將結果矩陣以上三角或下三角的形式存儲。但由于cublasDsyrk本身沒有將結果再對稱地運算,故可在CPU上完成對稱操作。經過該處理后,顯存耗用將減少一個矩陣G的大小,且計算量將減半。

2)當Ngg占用內存接近或超過顯存上限時,需要對矩陣G進行合適的分塊處理。分塊數量以及分塊方式的差異都會對矩陣相乘的計算速度產生影響,可根據實際情況靈活調整。

3.2 使用MPI減少內存占用

在矩陣相乘環節,采用矩陣分塊計算的方法僅能控制顯存,而無法壓縮矩陣占用的內存。為進一步提高相同硬件配置下MPI并行進程的總數或能解算的重力場階數上限,在上節的基礎上,使用MPI進一步對重力場反演任務進行劃分,通過去除直接計算Ngg的步驟實現分進程的內存壓縮。主要流程(圖4)如下:

圖4 使用MPI減少分進程內存峰值的流程

4)代入式(7)計算最終結果。

該計算流程省去了分進程中Ngg占用的大部分內存,當重力場階數較高時能將單進程內存峰值減少一半以上。同時,由于該過程不涉及文件讀寫,計算速度不會有明顯的下降。

4 算例分析

本節主要測試算法的加速效果。使用30 d仿真軌道數據模擬重力場月解的計算場景,解算不同階數的時變重力場模型,“真實”重力場模型選用GGM05S,并不加入任何誤差到仿真軌道數據中,先驗重力場模型選用EGM96,弧長統一為3 h,共240個弧段。

使用的計算機配置如下:處理器為Intel(R) Xeon(R) W-2235 CPU @ 3.80 GHz,內存為32 GB,六核心,顯卡為quadro P1000,顯存為4 GB,CUDA核心共640個。

4.1 MPI并行加速效果

本文程序基于Windows 10系統開發,使用微軟的Microsoft MPI進行并行計算,為適配Visual Studio 2010,選用的MPI版本為8.1。

以30階重力場解算為例,使用已經搭建完畢的重力場解算程序計算,僅修改并行的進程數量進行效率比對,統計結果見表3與圖5。

表3 啟動不同數量分進程時的重力場解算耗時比較

圖5 啟動不同數量分進程時的重力場解算耗時比較

結果顯示,隨著并行進程數量的增多,計算速度逐漸加快,但由于MPI并行本質上是邏輯的并行而非硬件層面的并行,計算速度仍然受到CPU性能的約束,所以進程數量越多,單進程的計算速度越慢。當進程數量超過CPU核心數目時,單進程的效率將大幅降低。在本實驗中,分進程數量從5增加至10時,單進程計算耗時明顯增加,而整體計算耗時并沒有顯著減少,內存峰值卻加倍,造成了不必要的資源占用。故使用MPI進行并行計算時,應綜合考慮實際任務需求與硬件性能進行適當的選擇,建議總進程數量不超過CPU核數。

4.2 CUDA對矩陣相乘的加速效果

本文使用的CUDA版本為8.0,主要用于大規模矩陣相乘環節。測試所使用矩陣G的大小均為2 160×N,N表示重力場位系數個數。分別計算30階、60階和120階的重力場模型對應的Ngg,使用的方法包括串行計算和CUDA并行計算,后者又包括簡單核函數并行以及cuBLAS庫函數并行,計算結果見表4。

表4 不同方法計算矩陣相乘耗時比較

結果顯示,使用CUDA并行計算矩陣相乘的速度遠快于串行計算,且矩陣越大,加速效果越明顯。但顯存容量通常要比內存更小,在使用CUDA時應綜合考慮顯存限制與計算性能,對大矩陣進行合適的分塊處理。

4.3 30階和120階重力場解算實例

使用上述方法搭建動力學法反演重力場模型的并行計算平臺。使用時長為30 d的仿真軌道數據進行測試,完成SST-HL模式下30階重力場模型的解算,共迭代3次。使用GGM05S重力場模型進行外符合精度評估,計算出的位系數階誤差RMS見圖6。

圖6 30階重力場解算結果階誤差RMS

結果顯示,30階重力場模型的解算精度達到10-14量級,且從第2次迭代開始收斂。可認為本文搭建的重力場并行解算平臺有足夠的低階重力場解算能力。本次共啟用10個MPI分進程,單個弧段耗時3～6 s,迭代3次耗時約430 s。串行模式下單弧段耗時105～110 s,迭代3次耗時約81 600 s。對比可知,本文的并行計算程序在計算30階重力場模型時可加速約180倍。

將重力場模型階數提升到120階,此時串行模式已經無法完成解算任務。使用與上一節仿真實驗相同的條件,計算出的位系數階誤差RMS見圖7。

圖7 120階重力場解算結果階誤差RMS

結果顯示,經過3次迭代得到的重力場模型的階誤差RMS可達10-13量級,證明了計算方法的正確性。該算例啟用5個分進程進行計算,每個弧段耗時100～200 s不等,全過程耗時約6.3 h。同時,該結果僅由單機計算,對于高性能的集群同樣適用。

5 結 語

本文針對基于最小二乘法的重力場模型解算過程,聯合MPI與CUDA實現了并行加速。主要得出以下幾點結論:

1)使用MPI和CUDA可分別在全局和局部實現重力場解算過程的并行加速。前者負責重力場反演任務拆分,加速比與進程數正相關,但建議進程數量不超過CPU核數;后者負責法矩陣相關計算,加速比可達380以上。

2)使用單機解算30階重力場時加速比可達180;解算120階重力場并迭代3次僅耗時約6.3 h,且該方法在高性能集群上同樣適用。