基于GPU集群加速的電阻率三維數值模擬

杜 偉， 吳小平

(1.中國科學技術大學 地球和空間科學學院 地震與地球內部物理實驗室，合肥 230026;2.蒙城地球物理國家野外科學觀測研究站，蒙城 233500)

0 引言

直流電法作為一種重要的勘探手段，廣泛應用于礦產資源勘察、水文地質等領域[1]。三維直流電法正演模擬作為反演的基礎，計算速度直接影響反演解釋的效率。正演主要的耗時步驟是大規模稀疏線性方程組的求解過程，在三維條件下共軛梯度迭代法較直接法更節省內存，計算時間更短，故被廣泛采用[2-4]。

迭代法最后將稀疏線性方程組的求解歸結為矩陣和向量的操作[14]。GPU作為一種向量處理機，適合大量并行線程同時進行，對于矩陣和向量的操作有天然的優勢。CUDA(Compute Unified Device Architecture，統一計算架構)是NVIDIA公司推出的一種編程模型，將GPU的通用計算普遍化。CUDA在2006年剛推出時綁定了C語言，之后經過多年發展，現已支持MATLAB、Python、Fortran等語言。2009年PGI公司聯合NVIDIA公司共同開發了支持FORTRAN的CUDA FORTRAN語言，保留了FORTRAN動態內存分配、模塊、向量操作等特性；添加了可操作GPU的語句；提供了cuBLAS和cuSPARSE等函數庫用于數值計算[5-6]。

國內、外學者在基于GPU的稀疏線性方程組求解方面做了很多研究，其工作主要圍繞如何加快稀疏矩陣向量(Sparse matrix-vector multiplication, SpMV)，如何有效存儲稀疏矩陣與如何使用GPU的不同層次的內存等問題展開[8,11]。目前，對于基于GPU集群的大規模線性方程組求解研究較少，主要有兩個問題需要解決：①矩陣預處理，傳統的IC/ILU/SSOR等預處理方法要求解三角線性方程組，其求解具有固有串行性[8]，不能有效地在多GPU上實現；②通信問題，NVIDIA推出的GPU-direct技術，實現了GPU在同一個節點內的GPU P2P直接通信，但跨節點的GPU數據傳輸需要將數據從GPU端拷貝到主機端，通過主機端的數據通信，然后拷貝到GPU端，數據通信延遲仍然是問題。

有學者使用了最簡單的Jacobi預處理，在GPU集群上實現了共軛梯度法[9-10]，但Jacobi預處理在大多稀疏線性方程中效果不理想，有效的并行預處理方法有待研究。基于節點內GPU-GPU直接訪問的技術，有學者實現了基于一個節點內的多GPU共軛梯度算法[9,11]。Cevahir 等[12]使用超圖劃分以及主機MPI通信實現了不同節點GPU的共軛梯度法；趙蓮等[13]實現了基于多級k路圖劃分的GPU集群的近似逆預處理共軛梯度法，其他的基于GPU集群的共軛梯度法研究較少。以上的節點間GPU通信均通過CPU內存中轉，效率較低。

MVAPICH2[19]是一個支持CUDA-Aware[7]的MPI通信實現，配合IB網絡(InfiniBand network)可實現支持跨節點GPU到GPU直接通信。筆者使用可以并行的對稱逐次超松弛迭代近似逆預處理技術，以及支持GPU-GPU通信的MVAPICH2 MPI，采用RCM算法將矩陣重新排序，減小矩陣帶寬，提高數據局部性，成功解決了以上問題，并在GPU集群上獲得了很好的擴展性。

1 三維電阻率數值模擬

當在地表使用雙電極供電時，三維直流電阻率控制方程為：

▽·(σ▽φ)=-I[δ(r-r+)-δ(r-r-)]

(1)

其中：φ為電勢；σ為電導率；r+和r-分別供電的正負電極位置；I為供電電流。采用非結構化網格離散研究區域，并在地表和其他邊界面添加邊界條件，采用有限元法求解上述控制方程的泛函極值最終得到大型稀疏線性方程組[3,23]：

Au=b

(2)

其中：u為三維離散網格各節點的電勢；A為最終的組裝矩陣；b為描述源分布的向量。A矩陣具有對稱稀疏正定的特性(SPD)，且其維度較大，方程組一般采用稀疏存儲并用迭代法求解。

2 SSORAI預處理

預處理可以降低原矩陣的條件數，加快迭代求解收斂。傳統的SSOR預處理[14]具有構造簡單、預處理效果較好的特點，但迭代過程中需要求解兩個三角矩陣線性方程組，并行度不高。并行的預處理方法有矩陣稀疏近似逆(sparse approximate inverse，SAI)、多項式預處理(Polynomial Preconditioners)等[14]方法。對稱逐次超松弛迭代近似逆[15](SSORAI)是一種近似逆方法，SSORAI預處理的思想即對SSOR預處理矩陣直接求逆，采用Neumann級數展開并截斷的方法直接構造出近似逆矩陣。

對于對稱矩陣A=E+D+F，考慮其SSOR預處理：

M=(D+ωE)D-1(D+ωF)=CCT

(3)

其中：D為對角線矩陣；E和F為上三角和下三角矩陣；0<ω<2為松弛因子，

C=(D+ωE)D-1/2

=D(I+ωD-1E)D-1/2

(4)

對式(4)求逆：

C-1=D1/2(I+ωD-1E)-1D-1

(5)

假設矩陣譜半徑ρ(D-1E)<1，由Neumann級數展開：

C-1=D1/2(I+ωD-1E)-1D-1

=D1/2(I-ωD-1E+(ωD-1E)2-…)D-1

(6)

截取低階項作為預處理矩陣的逆：

C-1?D1/2(I-ωD-1E)D-1=K

(7)

所以SSOR預處理矩陣的近似逆為：

(8)

圖1 SSOR預處理與SSORAI預處理效果對比Fig.1 Comparison of SSOR preconditioner and SSORAI preconditioner

圖2 RCM排序效果對比Fig.2 Result of RCM sorting(a)未排序前；(b) RCM算法重排效果

3 基于GPU集群的實現

3.1 RCM排序

使用SSORAI預處理后矩陣可以采用矩陣行塊劃分的方式直接并行[17]，但是GPU間的通信仍是瓶頸。分析GPU間通信的原因，發現由于矩陣元素分布過于分散(圖2(a))，迭代過程中SpMV需要使用到所有的向量，每次計算GPU都需要將所有的向量收集起來，導致通信量較大。為提高數據的局部性，減少GPU間數據通信，采用帶寬縮減的方式減少矩陣帶寬。常用的帶寬縮減算法有CM算法[21](Cuthill-McKee algorithm)、RCM算法[16](reverse Cuthill-McKee algorithm)、GPS算法[22]等算法，考慮到RCM算法具有簡單高效的特點，故采用之。

RCM算法將稀疏矩陣看作圖的鄰接矩陣，采用廣度優先遍歷的方式重新獲取圖的排列[16]，具體算法為：

1)確定一個起始節點r并令x1=r。

2)(主循環)從i=1、…、n，尋找所有與xi節點相鄰并且未被編號的節點，并依次按照節點度增序編號。

3)(反序)將序列反序得到新的排序y1、y2、…、yn，其中yi=xn-i+1。

RCM算法時間復雜度與鄰接矩陣的非零元個數成正比，計算時間較短,重排后的矩陣帶寬極大地減少(表1)。

表1 RCM排序時間和帶寬縮減對比Tab.1 Time of RCM sorting and bandwidth after the sorting

假設P為一個置換矩陣，則cond(PA)=cond(A),即置換矩陣作用于矩陣，變換后的矩陣條件數不變。RCM實質就是對矩陣進行置換操作，所以經過RCM排序后的矩陣條件數不變，即RCM排序不改變矩陣條件數，但可以極大地縮減矩陣帶寬。

3.2 GPU集群計算和通信

經過RCM算法排序后的矩陣具有帶寬有限的特征(圖2(b))，在GPU集群上采用矩陣行塊劃分的方式，同時對相應的向量也進行連續劃分，讓各自的矩陣和向量在GPU內部使用，具體劃分方式如圖3所示。在共軛梯度算法中，向量的加減、內積等操作可以完全并行。但進行矩陣向量乘積時，每個GPU不僅需要自身的向量，還需要秩相鄰的兩個GPU的半帶寬(half-band)的向量，例如秩為Rank 2的GPU進行SpMV時需要的向量為如圖3中A+B+C的組合，需要通過MPI通信獲取。MVAPICH2配合IB網絡支持GPU到GPU的直接通信，但是通信時GPU設備空閑，造成資源浪費，通信延遲仍然是一個瓶頸。利用GPU通信和計算重疊，可以進一步提高程序效率。

為了利用通信計算重疊，將每個GPU負責的矩陣塊劃分為3部分(圖3)矩陣塊A、B、C，同時將向量也劃分為三塊，并在秩向量前后各擴展half-band的空間，以接收存儲前后秩傳遞的數據。其中，各部分進行SPMV運算時，第一部分矩陣塊A需要前一個秩的半個帶寬向量，第二部分矩陣塊B只需要自身的向量，第三部分C需要后一個秩的半個帶寬的向量。為了編程的方便，在第一個秩前面和最后一個秩后面添加MPI的虛擬進程MPI_PROC_NULL[18](圖3)。這樣劃分以后，計算和通信順序為：計算矩陣塊B的SpMV，同時通信向量第一塊數據；等待通信結束，計算矩陣塊C的SpMV，同時通信向量第三塊數據；等待通信結束，計算矩陣塊A的SpMV，從而完成整塊矩陣的SpMV計算。

圖3 矩陣行塊劃分及向量連續劃分以及MPI通信示意圖Fig.3 Partition of matrix and vector (left) anf data transmission in each node with MPI (right)

為了統一編程，在RANK0前面和RANKN后面分別添加一個MPI的虛擬進程NULL[18]。第一次通信：將自身前半個帶寬向量通信到上一個秩向量列最后(矩陣第三塊進行SpMV計算需要)；第二次通信：將自身后半個帶寬向量通信到下一個秩列最前(矩陣第一塊進行SpMV計算需要)。

4 實驗結果

對三維直流電法模型研究區域使用TetGen[20]軟件進行非結構化網格剖分，使用不同參數生成不同的網格數目進行測試。程序完全由FORTRAN和CUDA FORTRAN編寫，實現了通用接口，可以方便地供現有FORTRAN程序調用。實驗使用CentOS 6.5操作系統(CPU：Intel Xeon E5-2670 @ 2.60GHz)，GPU為NVIDIA Tesla K40c(2880 core@745MHz，Peak single precision floating point performance: 4.29 Tflops，Peak double precision floating point performance: 1.43 Tflops)，驅動版本為352.39 ，開發工具為CUDA Toolkit 7.5。實驗中向量和稀疏矩陣運算調用NVIDIA函數庫cuBLAS和cuSPARSE。GPU集群通信使用MVAPICH2-2.1配合Mellanox InfiniBand。編譯使用PGI發布的支持CUDA FORTRAN的pgf90編譯器。為保證計算精度，計算一律采用雙精度浮點數和長整型整數。

圖4 兩層模型視電阻率計算結果Fig.4 Apparent resistivity of the two-layer model

圖5 低阻模型與高阻模型計算結果Fig.5 Results of the low resistivity model and the high resistivity model(a)低阻模型；(b)高阻模型

4.1 GPU計算結果及精度檢驗

為檢驗程序計算正確性和精度，本次實驗設計了具有解析解的兩層模型。取模型第一層電阻率為100 Ω·m ，層厚為 20 m ，第二層電阻率為10 Ω·m。點電源位于原點，電流強度為1 A ，限制迭代結束標志為 ‖Au-b‖/‖b‖ <10-10(所有計算均使用此收斂條件)。本程序計算的單極-單極視電阻率結果見圖4。與解析解對比，靠近源數值計算誤差較大，最大為3.39%，遠離源誤差不超過1.0%，計算視電阻率曲線與理論曲線趨勢一致。

另外本程序計算兩個典型的低阻和高阻模型(圖5)，取供電電極AB=4 000 m，即A(0,2000,0)，B(0,2000,0)，埋深h=200 m。圍巖電阻率為 100 Ω·m，低阻和高阻模型分別設置塊體電阻率為 10 Ω·m 和1 000 Ω·m。地表四極觀測的視電阻率如圖5所示，可明顯看出地下低阻體和高阻體的存在，視電阻率圖對稱性較好。從以上解析解對比及模型計算，可以看出本程序計算是正確可靠的。

4.2 多GPU擴展性研究

為研究GPU集群的擴展性，本次實驗選用了不同大小的網格做測試(表2)。實驗結果表明，在GPU上采用RCM排序后的SSORAI算法明顯快于未排序的計算。原因是由于排序后，矩陣帶寬明顯減小，調用cuSPARSE庫函數cusparseDcsrmv進行矩陣向量乘時，可以合并訪存以及減少線程的分歧，進而減少計算時間。建議在非結構化網格剖分中，盡量使用帶寬縮減算法減少帶寬，以提高矩陣向量乘的性能。

表2 RCM排序后SSORAI-CG計算時間對比Tab.2 Time of SSORAI-CG algorithm with RCM sorting and without sorting

圖6展示了筆者提出的方法對于大規模稀疏線性方程組的求解具有很好的擴展性。在較小矩陣規模時，由于通信占用時間較長，加速比較低；當矩陣規模增大時，加速比亦增大。在矩陣規模超過300×104時，2塊GPU卡相對1塊GPU卡有1.86倍的加速比，4塊GPU卡相對于1塊GPU卡有3.35倍加速比。本次試驗中，最大矩陣維度約為900×104，此時2塊GPU卡能能達到1.95倍加速比，4塊GPU卡能達到3.76倍加速比，效果較為理想。

圖6 GPU集群擴展性Fig.6 Expansibility of GPU clusters with our algorithm

5 總結

筆者采用SSORAI預處理共軛梯度法，求解直流電法正演生成的線性方程組，提高了預處理的并行性；該方法迭代次數較SSOR預處理略多，單次計算時間在小規模矩陣時較SSOR少；但隨著矩陣規模增大，SpMV耗時增加，所以有使用多GPU加速的必要。

針對多GPU上負載均衡的問題，使用RCM算法進行矩陣重新排序，減少矩陣帶寬，提高數據局部性。實驗結果表明，RCM算法的時間與矩陣規模呈線性關系，耗時較少；同時，由于排序后矩陣數據局部性較好，計算的時間較未排序有明顯地下降，隨著矩陣規模增大能有1.6倍～3.9倍的加速。建議對于大帶寬矩陣迭代求解時，先使用帶寬縮減算法對矩陣進行帶寬縮減，可以大幅提高計算性能。對排序后的矩陣采用矩陣行塊劃分，采用通信計算重疊，基于多GPU實現的預處理共軛梯度法獲得了很好的加速比和可擴展性。基于2塊GPU卡實現時，隨著矩陣規模不同，能獲得1.47～1.95的加速比；基于4塊GPU卡實現時，在矩陣規模較小計算量較小時，通信所占比重較大，只能獲得小于2倍的加速比，但隨著矩陣規模增大，矩陣規模大于三百萬時能獲得3.35倍以上的加速比，本次實驗中，在889×104規模時獲得了3.76倍加速比。

致謝：

本論文的數值計算得到了中國科學技術大學超級計算中心的計算支持和幫助。

[1] 傅良魁.電法勘探教程[M]. 北京：地質出版社, 1983. FU L K, Electrical prospecting tutorials[M] . Beijing: Geological Publishing House, 1983. (In Chinese)

[2] WU XIAOPING. A 3-D finite-element algorithm for DC resistivity modeling using shifted incomplete cholesky conjugate gradient method[J]. Geophys. J. Int., 2003, 154 (3): 947-956.

[3] 王威, 吳小平. 電阻率任意各向異性三維有限元快速正演[J]. 地球物理學進展, 2010,25 (4): 1365-1371. WANG W, WU X P . Rapid finite element resistivity forward modeling for 3D arbitrary anisotropic structures[J]. Progress in Geophysics, 2010,25(4):1365-1371. (In Chinese)

[4] 吳小平.利用共扼梯度算法的電阻率三維有限元正演[J]. 地球物理學報, 2003，46(3):428-432. WU X P . A 3D finite element resistivity forward modeling using conjugate gradient algorithm[J]. Chinese Journal Geophysics, 2003,46(3): 428-432. (In Chinese)

[5] RUETSCH G, FATICA M. CUDA Fortran for scientists and engineers: best practices for efficient CUDA Fortran programming[M]. Elsevier, 2013.

[6]CUDA Fortran programming guide and reference,2016.[online]Available:https://www.pgroup.com/doc/pgicudaforug.pdf

[7] WANG H, POTLURI S, BUREDDY D, et al. GPU-aware MPI on RDMA-enabled clusters: Design, implementation and evaluation[J]. Parallel and Distributed Systems, IEEE Transactions on, 2014, 25(10): 2595-2605.

[8] LI R, SAAD Y. GPU-accelerated preconditioned iterative linear solvers[J]. The Journal of Supercomputing, 2013, 63(2): 443-466.

[9] GEORGESCU S, OKUDA H. Conjugate gradients on multiple GPUs[J]. International Journal for Numerical Methods in Fluids, 2010, 64(10‐12): 1254-1273.

[10]GRIEBEL M, ZASPEL P. A multi-GPU accelerated solver for the three-dimensional two-phase incompressible Navier-Stokes equations[J]. Computer Science-Research and Development, 2010, 25(1-2): 65-73.

[11]CEVAHIR A, NUKADA A, MATSUOKA S. Fast conjugate gradients with multiple GPUs[M]. Computational Science-ICCS. Springer Berlin Heidelberg, 2009.

[12]CEVAHIR A, NUKADA A, MATSUOKA S. High performance conjugate gradient solver on multi-GPU clusters using hypergraph partitioning[J]. Computer Science-Research and Development, 2010, 25(1-2): 83-91.

[13]趙蓮, 趙永華, 陳堯, 等. GPU集群加速近似逆預條件CG并行求解器[J]. 計算機科學與探索, 2015, 9(9): 1084-1092. ZHAO L, ZHAO Y H, CHEN Y, et al. Approximate inverse preconditioned CG parallel solver on GPU cluster[J]. Journal of Frontiers of Computer Science and Technology, 2015, 9(9): 1084-1092. (In Chinese)

[14]SAAD Y. Iterative methods for sparse linear systems[M]. Siam, 2003.

[15]HELFENSTEIN R, KOKO J. Parallel preconditioned conjugate gradient algorithm on GPU[J]. Journal of Computational and Applied Mathematics, 2012, 236(15): 3584-3590.

[16]GEORGE A, LIU J W. Computer solution of large sparse positive definite[M]. Prentice Hall Professional Technical Reference. 1981.

[17]陳國良. 并行計算: 結構· 算法· 編程[M]. 北京：高等教育出版社, 2011. CHEN G L. Parallel computing: structure, algorithm and programming [M]. Beijing: Higher Education Press, 2011. (In Chinese)

[18]都志輝, 李三立, 陳渝, 等. 高性能計算之并行編程技術--MPI 并行程序設計[M]. 北京: 清華大學出版社, 2001. DU Z H,LI S L,CHEN Y,et al. Parallel programming technology in highperformance computing- MPI parallel programming[M]. Beijing: Tsinghua University Press, 2001.(In Chinese)

[19]MVAPICH2: http://mvapich.cse.ohio-state.edu/

[20]TetGen: http://wias-berlin.de/software/tetgen/

[21] CUTHILL E, MCKEE J. Reducing the bandwidth of sparse symmetric matrices[C]. Proceedings of the 1969 24th national conference. ACM, 1969: 157-172.

[22] GIBBS N E, POOLE, JR W G, STOCKMEYER P K. An algorithm for reducing the bandwidth and profile of a sparse matrix[J]. SIAM Journal on Numerical Analysis, 1976, 13(2): 236-250.

[23]WANG W.,WU X.,SPITZER K. Three-dimensional DC anisotropic resistivity modelling using finite elements on unstructured grids[J]. Geophys. J. Int,2013,193(2):734-746.