何康馨,席國江,陳 穎

(西南電子技術(shù)研究所 敏捷智能計算四川省重點實驗室,四川 成都 610065)

0 引言

通信技術(shù)和人工智能的飛速發(fā)展正在加速促進經(jīng)濟社會向數(shù)字化、網(wǎng)絡(luò)化、智能化轉(zhuǎn)型,在此過程中,通感算一體化已經(jīng)成為一種新的趨勢。為適應(yīng)這種新趨勢,敏捷智能計算四川省重點實驗室開發(fā)了“雨燕”敏捷智能集群計算系統(tǒng)[1-2]。該系統(tǒng)可靈活部署于飛機、船舶、車輛等載體上,提供接近大型超算中心的并行集群計算能力?；谠撓到y(tǒng),提出了一種并行廣義最小殘差(Generalized Minimal Residual,GMRES)方法來實現(xiàn)大規(guī)模線性方程組的快速求解。

大規(guī)模線性方程組的求解是計算數(shù)學(xué)和工程計算中的核心問題[3-4],也是很多人工智能應(yīng)用場景的計算基礎(chǔ)。在工程實踐中,由于存儲空間和計算能力的限制,無法采用矩陣分解進行消元計算,因此通常采用迭代法進行求解。常用的迭代方法包括雅可比迭代、高斯-塞德爾迭代、超松弛迭代、共軛梯度法等[5]。在工程應(yīng)用中最成功的方法是Krylov子空間方法,該方法在n維解空間的m(m

為了充分利用集群計算的優(yōu)勢,出現(xiàn)了大量針對GMRES方法并行實現(xiàn)的工作。Chronopoulos和Kim[8]提出了s-step GMRES方法,該方法通過生成s個同時搜索的方向向量來提高算法的并行性,并在Cray-2分級存儲超算上實現(xiàn)。文獻[9]探究了一種 預(yù)處理GMRES方法在多種并行架構(gòu)上的實現(xiàn),包括從分布式和共享內(nèi)存式的中央處理器(Central Processing Unit,CPU)架構(gòu)到圖形處理器(Graphics Processing Unit,GPU)架構(gòu)。還有一些工作針對經(jīng)典GMRES方法進行了并行上的優(yōu)化,文獻[10]通過構(gòu)造牛頓基來增強GMRES方法的數(shù)值穩(wěn)定性并借助QR分解提高其并行性。文獻[11]在分布式架構(gòu)上,通過通信和計算的覆蓋來提高并行效率。

本文基于敏捷集群計算系統(tǒng)提出了一種并行GMRES方法。首先采用矩陣向量乘法構(gòu)造一組Krylov子空間的非正交基,然后對由這些基為列向量構(gòu)成的矩陣做并行高瘦矩陣QR(Tall and Skinny QR,TSQR)分解[12],最后由此構(gòu)造出低維子空間中的一個最小二乘問題并求解,從而得到解在子空間中的最佳逼近。該方法的實現(xiàn)基于多核CPU集群,采用單指令多數(shù)據(jù)(Single Instruction Multiple Data,SIMD)編程,即同一個程序在不同CPU核針對不同數(shù)據(jù)執(zhí)行,數(shù)據(jù)傳遞基于消息傳遞接口(Message Passing Interface,MPI)實現(xiàn)。

1 串行GMRES方法

本節(jié)介紹經(jīng)典GMRES方法[13]?？紤]線性方程組Ax=b,其中A∈Rn×n,x,b∈Rn。Krylov子空間定義如下：

K(m;A;r0)=span(r0,Ar0,…,Am-1r0),

(1)

式中：m為K的維度,且m?n;r0=b-Ax0為初始殘差。給定方程組解的初值x0∈Rn。GMRES方法尋找近似解x∈x0+K,滿足：

b-Ax⊥AK。

(2)

該問題與求解如下優(yōu)化問題等價：

(3)

串行GMRES方法主要包括兩部分：第一部分是通過Arnoldi過程構(gòu)建Krylov子空間的一組標(biāo)準(zhǔn)正交基;第二部分是通過求解一個最小二乘問題在子空間中尋找解的最佳逼近。具體過程如算法1所示。

算法1 串行GMRES方法輸入：A,b,x0,m輸出：xr0=b-Ax0;β=r0;?Arnoldi過程v1=r0/β;forj=1,2,…,m doq=Avj;for i=1,2,…,j do hij=vTiq; q=q-hijvi;end forhj+1,j=q;vj+1=q/hj+1,j;end for?最小化 b-Ax,x∈x0+K求解J(y)=miny∈Rmβe1-Hmy得到y(tǒng)mx=x0+Vmym

在Arnoldi過程中構(gòu)造標(biāo)準(zhǔn)正交基Vm=[v1,v2,…,vm]的過程中有如下關(guān)系：

(4)

(5)

式中：e1=(1,0,…,0)T∈Rm+1。

(6)

在Arnoldi過程中,采用MGS正交化方法來構(gòu)造Krylov子空間的一組標(biāo)準(zhǔn)正交基。每當(dāng)生成 一個基向量,就要與前面已經(jīng)生成的所有向量進行正交化,此過程中存在嚴重的數(shù)據(jù)依賴性,影響算法的并行性。進一步地,需要存儲的n階向量v的數(shù)量隨m線性增長[14],所需計算量隨m平方增長[15],在求解大規(guī)模問題時會帶來巨大災(zāi)難。因此,在下節(jié)的并行GMRES方法中對Arnoldi過程做出改進。

2 并行GMRES方法

在并行GMRES方法中,通過兩步構(gòu)造Krylov子空間的一組基,通過并行矩陣向量乘法構(gòu)造一組非正交基Vm+1=[r0,Ar0,…,Amr0],對Vm+1做TSQR分解來代替Arnoldi過程中的正交化過程。并行GMRES方法的具體流程如算法2所示。算法2的并行過程主要有兩部分：一是并行構(gòu)造Krylov子空間;另一個是并行TSQR分解[16]。

算法2 并行GMRES方法輸入：A,b,x0,m輸出：xr0=b-Ax0; 并行計算Vm+1=[r0,Ar0,A2r0,…,Amr0];并行計算Vm+1=Qm+1Rm+1;Hm=R(,2：m+1);g=R(,1);求解J(y)=miny∈Rmg-Hmy獲得ymx=x0+Vmym

2.1 并行構(gòu)造Krylov子空間

構(gòu)建Krylov子空間的基向量Vm+1=[r0,Ar0,…,Amr0]需要大量的n階矩陣向量乘法,會消耗大量的計算和內(nèi)存資源,該部分的并行效果會直接影響問題的計算效率和可計算規(guī)模。Vm+1的并行構(gòu)造過程如圖1所示。

圖1 Krylov子空間的并行構(gòu)造Fig.1 Parallel construction of Krylov subspace

(7)

(8)

2.2 并行TSQR分解

為了降低數(shù)據(jù)依賴性,將Arnoldi過程中的正交化過程用Vm+1的QR分解來代替。由于m?n,可以采用并行TSQR算法。Vm+1分布在4個計算核中的并行TSQR方法如圖2所示。

圖2 并行TSQR方法Fig.2 Procees of TSQR method

Rij和Qij的下標(biāo)i表示計算核的編號,j表示 第j步TSQR分解過程。圖2中,Vm+1也按行分布式存儲在各個計算核中,在每個計算核中獨立地做局部QR分解,例如在第一步TSQR分解過程中,計算核1,得到：

V1=Q11R11。

(9)

根據(jù)V1的規(guī)模和稀疏性,選擇Gram-Schmidt正交化過程、Householder變換或Givens變換來實現(xiàn)該QR分解。將得到的上三角矩陣兩兩結(jié)合,并做 第二步QR分解,例如在第二步TSQR過程中,計算核1得到：

(10)

重復(fù)上述過程,直至得到最終的上三角矩陣Rm+1。

在并行TSQR分解過程中,采用MPI函數(shù)MPI_recv()和MPI_send()來實現(xiàn)Rij的數(shù)據(jù)分發(fā)。

3 數(shù)值實驗

本節(jié)基于“雨燕”敏捷集群計算平臺進行數(shù)值實驗,該平臺包含兩顆申威3231 CPU,共64個計算核。測試問題來源于對如下泊松方程的有限元離散[17]：

(11)

式中：f=2π2sin(πx)sin(πy)。

式(11)的變分形式為：

??u·?vdxdy=?f·vdxdy,

(12)

對計算區(qū)域Ω進行三角剖分,具體剖分形式如圖3所示。

圖3 Ω的三角剖分Fig.3 Triangulation of Ω

根據(jù)三角剖分,建立有限維子空間Vh= span{φ1,φ2,…,φN},其中φi是分段線性函數(shù),滿足：

(13)

(14)

將積分記作內(nèi)積,則式(14)變?yōu)椋?/p>

KU=F,

(15)

式中：Kij=(?φi,?φj),Fi=(f,φi)。

圖4 誤差與Krylov子空間維度m的關(guān)系

針對不同規(guī)模的K進行數(shù)值實驗。圖5和圖6展示了一步并行GMRES方法中,CPU時間t和并行效率E分別與計算核數(shù)目p的關(guān)系。并行效率E表示單個計算核的平均加速效果,定義如下：

(a) n=1 600

(c) n=10 000

(a) n=1 600

(c) n=10 000

(16)

式中：t1表示串行GMRES方法消耗的CPU時間,tp表示采用p個計算核的并行GMRES方法所消耗的CPU時間。

算法消耗的CPU時間包括計算時間和通信時間。隨著計算核數(shù)目的增加,一方面,每個計算核處理的數(shù)據(jù)量下降,因此計算時間減少;另一方面,計算核之間的通信量增加,因此通信時間增加。因此圖5中的CPU時間隨著計算核數(shù)目p的增長,減小得越來越慢。當(dāng)p足夠大時,CPU時間甚至可能增加。算法的并行效率與問題規(guī)模相關(guān),當(dāng)問題規(guī)模給定時,如圖6所示,算法的并行效率會隨著計算核數(shù)目的增加而下降。

綜上所述,應(yīng)該根據(jù)問題規(guī)模選擇合適的計算核數(shù)目。不僅需要考慮消耗的CPU時間,還需要考慮并行效率。因為雖然增大p可以減少CPU時間消耗,但若并行效率過低,則是對計算資源的浪費。因此需要選擇合適的p,使得CPU時間消耗較低,同時并行效率較高。據(jù)此,當(dāng)K的規(guī)模為n=10 000、 2 500、1 600時,選擇的計算核數(shù)目分別為16、8、4。

基于該敏捷集群計算平臺,并行GMRES方法實現(xiàn)了10萬規(guī)模線性方程組的快速求解,具體參數(shù)及結(jié)果如表1所示。

表1 10萬規(guī)模線性方程組求解參數(shù)及結(jié)果Tab.1 Computational setting and results of large linear system of order 105

4 結(jié)束語

為了能夠充分發(fā)揮通感算一體化的優(yōu)勢,基于多核敏捷集群計算系統(tǒng),研究了大規(guī)模線性方程組Ax=b的并行解法,提出了并行GMRES方法。并行了串行GMRES方法中計算量最大的部分：一個是構(gòu)建Krylov子空間過程中矩陣向量乘法的并行;另一個是采用并行TSQR分解代替Arnoldi過程中的MGS正交化過程。通過二維泊松方程的有限元離散,得到不同規(guī)模的線性方程組,來驗證算法的有效性。

并行GMRES方法本質(zhì)上減小了數(shù)據(jù)依賴性,提高了并行效率。因此,在求解大規(guī)模線性方程組時,可以充分發(fā)揮多核敏捷集群計算系統(tǒng)的計算和存儲優(yōu)勢。在未來的工作中,為進一步提升大規(guī)模線性方程組的求解效率,一方面可以考慮混合精度算法,另一方面可以開發(fā)基于異構(gòu)計算系統(tǒng)的并行計算方法。