張 琨，賈金芳，黃建強(qiáng)，2，王曉英，嚴(yán)文昕

1(青海大學(xué) 計(jì)算機(jī)技術(shù)與應(yīng)用系，西寧 810016)

2(清華大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)系，北京 100084)

E-mail：qhu_jf@163.com

1 引 言

線性方程組求解是許多科學(xué)計(jì)算和工程應(yīng)用中的核心計(jì)算內(nèi)容，其目前主要有直接法和迭代法兩種求解方式[1-3].直接法由于計(jì)算機(jī)資源的限制難以有效求解大規(guī)模線性方程組，而迭代法具有存儲(chǔ)需求低，系數(shù)矩陣在迭代時(shí)保持穩(wěn)定等優(yōu)點(diǎn)，所以其常用于求解大規(guī)模線性方程組.目前迭代法已經(jīng)由Jacobi方法、SOR方法等經(jīng)典迭代法過(guò)渡到了共軛梯度(CG)方法和廣義共軛殘差(GCR)方法等效率更高的Krylov子空間迭代法[4，5].同時(shí)，為了改善迭代法的收斂性，也相繼產(chǎn)生了多種多樣的預(yù)處理技術(shù)[6-8].

由于求解問(wèn)題規(guī)模不斷增大，單處理器已經(jīng)不能滿(mǎn)足大規(guī)模計(jì)算的需求，因此研究大規(guī)模稀疏線性方程組的快速求解方案具有重要意義.并行計(jì)算是高效求解大規(guī)模稀疏線性方程組的重要途經(jīng)，其主要有多核CPU、GPGPU、集群3種類(lèi)型[9，10].許多學(xué)者為探究不同平臺(tái)的高性能迭代法實(shí)現(xiàn)方案展開(kāi)了相關(guān)研究.Bohacek J等[11]針對(duì)熱擴(kuò)散問(wèn)題在GPU上利用多項(xiàng)式預(yù)處理對(duì)CG法進(jìn)行并行優(yōu)化，與傳統(tǒng)CFD代碼相比最高能夠獲得16000倍的加速，同時(shí)還指出了高效的并行算法不僅要考慮硬件設(shè)備的更新，還要考慮算法本身的改進(jìn).Anzt H等[12]提出一種應(yīng)用于線性系統(tǒng)求解的不完全稀疏近似逆(ISAI)預(yù)處理方法，該方法可以并行生成預(yù)條件子，而且與Block-Jacobi預(yù)處理方法相比可以獲得更好的性能.Jiao Y Y等[13]為滿(mǎn)足地質(zhì)領(lǐng)域求解線性方程組的需要，在MPI、OpenMP、MPI+OpenMP 3種方案中測(cè)試可行且高效的CG法并行方案來(lái)解決球面非連續(xù)變形分析問(wèn)題，在模擬隧道坍塌的實(shí)驗(yàn)中加速比達(dá)到7x左右.Suryanarayana P等[14]基于PETSc并行庫(kù)實(shí)現(xiàn)了一種新的AAR迭代法，并與GMRES、Bi-CGSTAB等多種迭代法進(jìn)行了對(duì)比，在相同的預(yù)處理?xiàng)l件下，新方法能夠更快地完成求解.Bernaschi M等[15]實(shí)現(xiàn)了基于GPU的包含AMG預(yù)處理部分和求解部分的線性求解器，其通過(guò)一種兼容加權(quán)匹配聚合算法實(shí)現(xiàn)最優(yōu)的核函數(shù)配置，從而有效利用了GPU的性能.Aliaga J I等[16]對(duì)ILUPACK中的GMRES算法進(jìn)行了改進(jìn)，優(yōu)化了其中的預(yù)處理子應(yīng)用部分及MGSO部分，使這兩部分的加速比分別提升至3x和10x.Sanjuan G等[17]將Schwarz交替域分解技術(shù)與MPI+OpenMP混合并行技術(shù)結(jié)合，使風(fēng)場(chǎng)計(jì)算中的CG迭代法求解速度提升至近9x.

綜上所述，目前的性能優(yōu)化技術(shù)大多集中在算法改進(jìn)或者單一并行方式領(lǐng)域，而基于異構(gòu)平臺(tái)的多機(jī)并行優(yōu)化相對(duì)較少.本文在CPU+GPU異構(gòu)平臺(tái)上實(shí)現(xiàn)了基于MPI+CUDA混合并行的預(yù)處理共軛梯度算法，并對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比分析，發(fā)現(xiàn)MPI+CUDA混合并行的預(yù)處理共軛梯度算法能夠更高效地求解大規(guī)模對(duì)稱(chēng)正定線性系統(tǒng).

2 線性方程組與預(yù)處理共軛梯度算法

2.1 線性方程組

線性方程組的求解廣泛存在于工程問(wèn)題中.就大規(guī)模數(shù)值計(jì)算而言，使用迭代法求解線性方程組時(shí)涉及到的矩陣及向量運(yùn)算會(huì)占據(jù)大量的計(jì)算時(shí)間和硬件資源.隨著求解問(wèn)題規(guī)模的不斷增大，如何加快大規(guī)模稀疏線性方程組的求解過(guò)程顯得愈加重要.一個(gè)n元線性方程組如公式(1)所示，其中m為方程組個(gè)數(shù)，am，n為方程系數(shù)、x為方程解、b為常數(shù)項(xiàng).

(1)

將線性方程組轉(zhuǎn)換為線性代數(shù)中的概念則可以表示為Ax=b的形式，如公式(2)所示，其中A是m×n的系數(shù)矩陣，x是解向量，b是右端向量.

(2)

2.2 預(yù)處理共軛梯度算法及性能分析

2.2.1 預(yù)處理共軛梯度算法

共軛梯度算法是解決對(duì)稱(chēng)正定線性系統(tǒng)問(wèn)題最有效的迭代方法之一，具有存儲(chǔ)量要求少、不必事先估計(jì)某些參數(shù)、編程簡(jiǎn)單等優(yōu)點(diǎn)[18].為改善迭代法的收斂性，目前主要利用預(yù)處理技術(shù)將原方程Ax=b替換為M(-1)Ax=M(-1)b等形式進(jìn)行求解，其中較為常用的預(yù)條件子有不完全分解預(yù)條件子、稀疏近似逆預(yù)條件子等[19-21].預(yù)處理共軛梯度算法的具體內(nèi)容如算法1所示.

算法1. 預(yù)處理共軛梯度算法

輸入：系數(shù)矩陣A，預(yù)條件子M，初始解x0，右端向量b.

輸出：近似解x.

1.r0=b-Ax0，z0=M-1r0，p1=z0

2.fori= 1： Maximum number of iterationsdo

3.αi= (ri-1，zi-1)/(Api，pi)

4.xi=xi-1+αipi

5.ri=ri-1-αiApi

6.ifconvergedthenexit

7.zi=M-1ri

8.βi= (ri，zi)/(ri-1，zi-1)

9.pi+1=zi+βipi

10.endfor

2.2.2 預(yù)處理共軛梯度算法性能分析

每次迭代過(guò)程中，預(yù)處理共軛梯度算法需要進(jìn)行1次稀疏矩陣向量乘，5次向量?jī)?nèi)積(包括1次殘差計(jì)算)，1次預(yù)條件子應(yīng)用，3次向量數(shù)乘.在預(yù)處理共軛梯度算法的主要計(jì)算量分布測(cè)試中，預(yù)處理部分使用ILU預(yù)條件子進(jìn)行測(cè)試，占比55.51%，稀疏矩陣向量乘占比40.01%，向量?jī)?nèi)積占比2.36%，向量數(shù)乘占比2.12%.由于使用不同的預(yù)條件子的計(jì)算開(kāi)銷(xiāo)會(huì)有較大差異，本文主要針對(duì)其余計(jì)算部分進(jìn)行性能優(yōu)化.單機(jī)設(shè)備的計(jì)算性能有限，而且還有可能無(wú)法滿(mǎn)足大規(guī)模數(shù)據(jù)的存儲(chǔ)需求，所以目前通常采用并行的方式對(duì)預(yù)處理共軛梯度算法進(jìn)行求解及優(yōu)化.

3 預(yù)處理共軛梯度算法并行優(yōu)化

3.1 優(yōu)化框架

并行算法在優(yōu)化過(guò)程中為保證稀疏矩陣存儲(chǔ)性能一致，統(tǒng)一使用CSR(Compressed Sparse Row)格式保存稀疏矩陣[22，23].實(shí)驗(yàn)在異構(gòu)平臺(tái)上進(jìn)行測(cè)試，CPU端并行技術(shù)使用MPI，GPU端并行技術(shù)使用CUDA.異構(gòu)并行的預(yù)處理共軛梯度算法執(zhí)行步驟為：

1)MPI進(jìn)程讀取數(shù)據(jù)并初始化CPU和GPU內(nèi)存空間；

2)CUDA核函數(shù)計(jì)算r，然后MPI進(jìn)程進(jìn)行數(shù)據(jù)通信；

3)MPI進(jìn)程進(jìn)行預(yù)條件子應(yīng)用，計(jì)算z和p；

4)CUDA核函數(shù)計(jì)算稀疏矩陣向量乘和向量?jī)?nèi)積，然后MPI進(jìn)程進(jìn)行通信得到α；

5)CUDA核函數(shù)計(jì)算向量數(shù)乘，更新x和r；

6)CUDA核函數(shù)計(jì)算殘差Res，然后MPI進(jìn)程進(jìn)行數(shù)據(jù)通信及迭代退出判斷；

7)MPI進(jìn)程進(jìn)行預(yù)條件子應(yīng)用，計(jì)算z；

8)CUDA核函數(shù)計(jì)算向量?jī)?nèi)積，然后MPI進(jìn)程進(jìn)行數(shù)據(jù)通信得到β；

9)CUDA核函數(shù)計(jì)算向量數(shù)乘，更新p，然后MPI進(jìn)程進(jìn)行數(shù)據(jù)通信并返回步驟4.

實(shí)驗(yàn)方案分為3個(gè)部分，首先在預(yù)處理優(yōu)化部分，方案對(duì)ILU預(yù)條件子和Jacobi預(yù)條件子的加速性能進(jìn)行對(duì)比分析；然后在MPI并行部分，方案將可以重疊的部分計(jì)算操作和通信操作進(jìn)行了整合；最后在MPI+CUDA混合并行部分，方案針對(duì)異構(gòu)并行中存在的數(shù)據(jù)傳輸問(wèn)題、訪存問(wèn)題等，分別通過(guò)使用頁(yè)鎖定內(nèi)存、改變線程任務(wù)分配方式及利用共享存儲(chǔ)器的措施進(jìn)行優(yōu)化.

3.2 預(yù)處理優(yōu)化

本文所測(cè)試的預(yù)條件子為ILU預(yù)條件子和Jacobi預(yù)條件子.Jacobi預(yù)條件子的構(gòu)造過(guò)程相對(duì)簡(jiǎn)單，最后生成的預(yù)條件子M即為DIA(A).ILU預(yù)條件子將系數(shù)矩陣A分解為L(zhǎng)U-R的形式，分解之后的預(yù)條件子可以滿(mǎn)足不同的條件P.ILU預(yù)條件子的構(gòu)造過(guò)程如算法2所示.

算法2. ILU分解

輸入：系數(shù)矩陣A，條件P.

輸出：預(yù)條件子M.

1.fori= 2，3，…，ndo

2.fork= 1，2，…，i-1and(i，k)?Pdo

3.Aik=Aik/Akk

4.forj=k+1，…，nand(i，j)?Pdo

5.Aij=Aij-Aik×Akj

6.endfor

7.endfor

8.endfor

3.3 MPI并行通信優(yōu)化

MPI并行方式中各個(gè)進(jìn)程間的通信開(kāi)銷(xiāo)會(huì)隨著進(jìn)程規(guī)模的增大而增大[24].在預(yù)處理共軛梯度算法中，并行任務(wù)主要在每一個(gè)迭代步內(nèi)執(zhí)行，其中MPI主要負(fù)責(zé)任務(wù)的粗粒度劃分、進(jìn)程間數(shù)據(jù)通信、進(jìn)程內(nèi)部數(shù)據(jù)拷貝及非并行部分的計(jì)算.

在MPI并行算法中，每次迭代都會(huì)進(jìn)行數(shù)據(jù)通信，但是除計(jì)算Ap所進(jìn)行的稀疏矩陣向量乘(SpMV)通信之外，其他向量?jī)?nèi)積(Dot)的部分通信結(jié)果可以延遲參與計(jì)算.為減少算法的通信開(kāi)銷(xiāo)，每次迭代步內(nèi)將部分計(jì)算操作與通信操作進(jìn)行重疊.由于計(jì)算Ap的SpMV操作需要保證p向量在各個(gè)進(jìn)程中是完整的，而且也不存在可以與通信進(jìn)行重疊的計(jì)算任務(wù)，所以在SpMV操作之前使用阻塞的Allgather函數(shù)進(jìn)行數(shù)據(jù)通信.在后續(xù)3組需要通信的Dot操作中，通信結(jié)果存在可以延遲計(jì)算的部分，同時(shí)有相應(yīng)的計(jì)算任務(wù)可以與通信進(jìn)行重疊，所以使用非阻塞的Iallreduce函數(shù)進(jìn)行數(shù)據(jù)通信.在將計(jì)算操作和通信操作進(jìn)行重疊之后，通信優(yōu)化算法可以減少使用阻塞通信函數(shù)帶來(lái)的額外開(kāi)銷(xiāo).

3.4 GPU并行優(yōu)化

3.4.1 數(shù)據(jù)傳輸優(yōu)化

CPU與GPU協(xié)同計(jì)算時(shí)最主要的開(kāi)銷(xiāo)在于兩者之間的數(shù)據(jù)傳輸，并行程序只對(duì)部分必須進(jìn)行全局更新的向量進(jìn)行數(shù)據(jù)傳輸.在每次迭代過(guò)程中無(wú)法避免的數(shù)據(jù)傳輸內(nèi)容包括向量p、向量r、向量z、殘差Res以及向量?jī)?nèi)積的中間結(jié)果subDot.對(duì)于這些需要在CPU與GPU間頻繁傳輸?shù)臄?shù)據(jù)在開(kāi)辟內(nèi)存空間時(shí)采用頁(yè)鎖定內(nèi)存(pinned memory)，因?yàn)椴僮飨到y(tǒng)不會(huì)對(duì)頁(yè)鎖定內(nèi)存執(zhí)行分頁(yè)和交換操作，使得該內(nèi)存能夠始終駐留在物理內(nèi)存中，而GPU可以通過(guò)直接內(nèi)存訪問(wèn)方式直接在CPU和GPU間拷貝數(shù)據(jù)，從而提高數(shù)據(jù)傳輸效率.

3.4.2 訪存優(yōu)化

將算法主要計(jì)算操作移植到GPU上之后，全局存儲(chǔ)器的合并訪存(coalesced access)特性也會(huì)對(duì)程序性能造成一定影響.CUDA模型的最小執(zhí)行單位為線程束(warp)，當(dāng)同一線程束內(nèi)的線程訪問(wèn)連續(xù)的存儲(chǔ)空間時(shí)，能夠減少訪存操作的次數(shù).向量計(jì)算中的一維數(shù)據(jù)能夠保證線程束中的線程訪存是連續(xù)的，但在矩陣計(jì)算中則需要根據(jù)使用的稀疏矩陣存儲(chǔ)格式調(diào)整線程計(jì)算任務(wù)的分配策略.實(shí)驗(yàn)采用的稀疏矩陣存儲(chǔ)格式為CSR，該格式使用3個(gè)數(shù)組分別保存稀疏矩陣每行的偏移量、列索引及非零元素值.計(jì)算任務(wù)以線程為最小單位分配時(shí)，會(huì)出現(xiàn)同一線程束內(nèi)的線程訪問(wèn)不連續(xù)內(nèi)存空間的情況，從而造成訪存性能下降.為保證矩陣計(jì)算過(guò)程中的訪存連續(xù)性，計(jì)算任務(wù)以線程束為最小單位進(jìn)行分配.以線程束為最小單位調(diào)度任務(wù)之后，同一線程束內(nèi)的線程訪問(wèn)連續(xù)的內(nèi)存空間，從而可以利用全局存儲(chǔ)器的合并訪存特性提高訪存性能.

3.4.3 多級(jí)存儲(chǔ)器優(yōu)化

圖1 內(nèi)積并行Fig.1 Dot parallel mode

4 實(shí)驗(yàn)與分析

4.1 實(shí)驗(yàn)環(huán)境和數(shù)據(jù)

實(shí)驗(yàn)在4個(gè)高性能節(jié)點(diǎn)上進(jìn)行測(cè)試，每個(gè)節(jié)點(diǎn)兩塊CPU，型號(hào)為Intel(R) Xeon(R) CPU E5-2692 v2 @ 2.20GHz，內(nèi)存為64GB；節(jié)點(diǎn)配備GPU為T(mén)esla T4，顯存為16GB；CUDA模型版本為v10.1.實(shí)驗(yàn)所測(cè)試的數(shù)據(jù)來(lái)自SuiteSparse Matrix Collection數(shù)據(jù)集[25]，該數(shù)據(jù)集包含大量來(lái)自實(shí)際工程中的稀疏矩陣數(shù)據(jù).如表1所示，實(shí)驗(yàn)選用了計(jì)算流體力學(xué)問(wèn)題、電路仿真問(wèn)題、結(jié)構(gòu)性問(wèn)題等8個(gè)領(lǐng)域共11個(gè)測(cè)試數(shù)據(jù)，同時(shí)測(cè)試數(shù)據(jù)的元素量從3087898到117406044不等，具有一定數(shù)據(jù)規(guī)?？缍?

表1 矩陣數(shù)據(jù)Table 1 Matrix data

4.2 實(shí)驗(yàn)結(jié)果與分析

4.2.1 預(yù)處理結(jié)果

對(duì)于條件數(shù)過(guò)大的系數(shù)矩陣，使用迭代法進(jìn)行求解時(shí)通常會(huì)出現(xiàn)收斂速度慢甚至不收斂的情況.雖然使用預(yù)處理技術(shù)會(huì)增加部分計(jì)算開(kāi)銷(xiāo)，但為了提高迭代法的求解效率，利用預(yù)處理技術(shù)來(lái)改善迭代法的收斂性是有必要的.

實(shí)驗(yàn)測(cè)試了ILU預(yù)條件子和Jacobi預(yù)條件子兩種方式對(duì)CG算法收斂性的影響.由于不同領(lǐng)域問(wèn)題對(duì)誤差精度要求不一致，將問(wèn)題精度設(shè)置為統(tǒng)一標(biāo)準(zhǔn)會(huì)出現(xiàn)不同測(cè)試數(shù)據(jù)差異過(guò)大的情況，因此實(shí)驗(yàn)對(duì)相同迭代次數(shù)下的誤差變化進(jìn)行對(duì)比.在迭代200次的條件下，無(wú)預(yù)處理算法和兩種預(yù)處理算法的最終殘差如圖2所示.從圖中可以看出，在加速迭代法收斂的方面ILU預(yù)條件子比Jacobi預(yù)條件子效果更好.ILU預(yù)條件子構(gòu)造和使用的成本都比Jacobi預(yù)條件子更高，而且ILU預(yù)條件子在應(yīng)用時(shí)也沒(méi)有與Jacobi預(yù)條件子相當(dāng)?shù)牟⑿行?因此對(duì)于預(yù)處理子的選擇需要綜合考慮應(yīng)用成本和誤差精度要求兩個(gè)方面，當(dāng)誤差精度要求較低且需要控制應(yīng)用成本時(shí)可以選擇使用Jacobi預(yù)條件子，當(dāng)對(duì)誤差精度要求較高時(shí)則需要使用ILU預(yù)條件子.

圖2 預(yù)條件子對(duì)比Fig.2 Comparison of preconditioner

4.2.2 無(wú)預(yù)處理算法并行結(jié)果

無(wú)預(yù)處理算法中的主要計(jì)算內(nèi)容為向量數(shù)乘、向量?jī)?nèi)積、稀疏矩陣向量乘.串行算法和3種并行算法的運(yùn)行性能對(duì)比如圖3所示.為與MPI+CUDA混合并行進(jìn)程數(shù)量保持一致，MPI并行一共使用4個(gè)節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)啟動(dòng)一個(gè)進(jìn)程.由于使用MPI并行會(huì)存在部分額外的開(kāi)銷(xiāo)，實(shí)驗(yàn)所測(cè)試的矩陣數(shù)據(jù)平均加速比為2.99x(最高3.41x).MPI并行的開(kāi)銷(xiāo)主要來(lái)自于各進(jìn)程間的數(shù)據(jù)通信，不同算例的通信次數(shù)是一致的，通信開(kāi)銷(xiāo)有差異主要是因?yàn)闇y(cè)試算例的數(shù)據(jù)規(guī)模不同.當(dāng)測(cè)試矩陣非零元素?cái)?shù)量不多而階數(shù)較大時(shí)，會(huì)出現(xiàn)矩陣每行計(jì)算量(Density)過(guò)小的情況，最后導(dǎo)致MPI并行性能提升不明顯.MPI并行中矩陣每行計(jì)算量較小的算例如G3_circuit(4.83)和thermal2(6.99)加速比只有2.5x左右，而矩陣每行計(jì)算量較大的算例如Emilia_023(44.42)、BenElechi1(53.48)加速比能達(dá)到3x以上.

圖3 無(wú)預(yù)處理算法性能對(duì)比Fig.3 Performance comparison of original algorithm

MPI+CUDA混合并行在MPI并行的基礎(chǔ)上將計(jì)算任務(wù)移植到GPU上執(zhí)行，該方式并沒(méi)有更改MPI并行各進(jìn)程間的通信模式，但是可以提高各進(jìn)程計(jì)算部分的性能.MPI+CUDA混合并行將計(jì)算任務(wù)進(jìn)行GPU移植之后充分發(fā)揮了GPU的多線程能力，在測(cè)試數(shù)據(jù)上的平均加速比為17.25x(最高27.12x).MPI+CUDA混合并行的性能變化趨勢(shì)與MPI并行基本一致，但有一點(diǎn)區(qū)別是總的計(jì)算量不能過(guò)小，否則不能發(fā)揮GPU的多線程優(yōu)勢(shì)，甚至可能由于CPU與GPU之間數(shù)據(jù)傳輸開(kāi)銷(xiāo)造成性能下降的情況.如測(cè)試數(shù)據(jù)cfd2(25.02)，雖然矩陣每行都有一定的數(shù)據(jù)量，但由于總的計(jì)算量過(guò)小，最后MPI+CUDA混合并行的加速比不足10x.CUDALib并行算法使用Nvidia并行計(jì)算庫(kù)中的cuSPARSE和cuBLAS提供的代數(shù)計(jì)算接口函數(shù)進(jìn)行實(shí)現(xiàn)，在各個(gè)測(cè)試數(shù)據(jù)上的平均加速比為13.13x(最高19.42x).

4.2.3 ILU預(yù)處理算法并行結(jié)果

ILU預(yù)處理算法在無(wú)預(yù)處理算法的基礎(chǔ)上增加了預(yù)條件子應(yīng)用的計(jì)算內(nèi)容.ILU預(yù)條件子應(yīng)用需要在每次迭代過(guò)程中增加兩次三角方程求解.因?yàn)槿欠匠糖蠼庥?jì)算具有明顯的串行性質(zhì)，實(shí)驗(yàn)將這一部分計(jì)算內(nèi)容按照串行的方式執(zhí)行.由于存在一部分串行計(jì)算內(nèi)容，ILU預(yù)處理算法與無(wú)預(yù)處理算法相比整體性能出現(xiàn)下降的現(xiàn)象.串行算法、MPI并行算法、CUDALib并行算法及MPI+CUDA混合并行算法的性能對(duì)比如圖4所示.在實(shí)驗(yàn)測(cè)試的矩陣數(shù)據(jù)上MPI并行算法性能平均提升47%(最高54%)，CUDALIB并行算法性能平均提升71%(最高92%)，MPI+CUDA混合并行算法性能平均提升84%(最高97%).ILU預(yù)處理算法在應(yīng)用預(yù)條件子時(shí)的計(jì)算量與SpMV相當(dāng)，而且由于預(yù)條件子應(yīng)用沒(méi)有進(jìn)行并行處理，所以當(dāng)算法中其他部分如SpMV、Dot、AXPY等通過(guò)并行方式提高計(jì)算性能時(shí)，預(yù)條件子應(yīng)用會(huì)占據(jù)程序執(zhí)行的主要時(shí)間.

圖4 ILU預(yù)處理算法性能對(duì)比Fig.4 Performance comparison of ILU preconditioned algorithm

4.2.4 Jacobi預(yù)處理算法并行結(jié)果

Jacobi預(yù)條件子的應(yīng)用比ILU預(yù)條件子的應(yīng)用計(jì)算量更少，每次迭代過(guò)程中應(yīng)用Jacobi預(yù)條件子的計(jì)算量與一次向量數(shù)乘相當(dāng)，所以在Jacobi預(yù)處理并行算法中主要的計(jì)算開(kāi)銷(xiāo)還是集中在稀疏矩陣向量乘部分.因?yàn)榇腥蝿?wù)所占據(jù)的計(jì)算量較少，Jacobi預(yù)處理算法在各矩陣數(shù)據(jù)上的運(yùn)行時(shí)間與無(wú)預(yù)處理算法整體接近，能獲得較好的加速效果.串行算法、MPI并行算法、CUDALib并行算法及MPI+CUDA混合并行算法的性能對(duì)比如圖5所示.在實(shí)驗(yàn)測(cè)試的矩陣數(shù)據(jù)上MPI并行算法平均加速比為2.92x(最高3.46x)，CUDALib并行算法平均加速比為9.56x(最高14.3x)，MPI+CUDA混合并行算法平均加速比為12.76x(最高19.45x).

圖5 Jacobi預(yù)處理算法性能對(duì)比Fig.5 Performance comparison of Jacobi preconditioned algorithm

4.2.5 可擴(kuò)展性分析

如圖6所示，實(shí)驗(yàn)測(cè)試了不同并行算法加速比隨著節(jié)點(diǎn)數(shù)目增加的變化趨勢(shì).在使用MPI并行方式的算法中，隨著節(jié)點(diǎn)數(shù)目的增加，無(wú)預(yù)處理算法和Jacobi預(yù)處理算法的加速比保持上升趨勢(shì)，但I(xiàn)LU預(yù)處理算法的加速比提升并不明顯，這是由于ILU預(yù)處理算法中的預(yù)條件子應(yīng)用會(huì)增加較多額外的開(kāi)銷(xiāo)，從而限制了整體加速比的提升.在使用MPI+CUDA混合并行方式的算法中，隨著節(jié)點(diǎn)數(shù)目的增加，無(wú)預(yù)處理算法和Jacobi預(yù)處理算法的加速比同樣保持上升趨勢(shì)，并且與MPI并行方式的對(duì)應(yīng)算法加速比相比還有大幅度的提升，這得益于GPU的多線程優(yōu)勢(shì)，而ILU預(yù)處理算法的加速比由于預(yù)條件子應(yīng)用開(kāi)銷(xiāo)的存在依然沒(méi)有明顯提升.

圖6 可擴(kuò)展性對(duì)比Fig.6 Comparison of scalability

5 結(jié) 語(yǔ)

大規(guī)模稀疏線性方程組的求解問(wèn)題一直都是人們研究的重點(diǎn)方向.隨著各種工程問(wèn)題的求解規(guī)模增加，單機(jī)設(shè)備已經(jīng)難以高效完成大規(guī)模計(jì)算任務(wù).本文首先實(shí)現(xiàn)了3種不同的預(yù)處理共軛梯度并行算法，然后對(duì)結(jié)果進(jìn)行了對(duì)比分析.實(shí)驗(yàn)結(jié)果表明，在MPI并行、CUDALib并行和MPI+CUDA混合并行3種方式中，MPI+CUDA混合并行方式能夠充分發(fā)揮多機(jī)異構(gòu)設(shè)備的性能，達(dá)到更高的計(jì)算效率.在下一步工作中，我們將進(jìn)行更多預(yù)條件子異構(gòu)性能優(yōu)化的研究.