吳建平 銀?？?彭 軍 楊錦輝

(國防科學(xué)技術(shù)大學(xué)海洋科學(xué)與工程研究院 湖南 長沙 410073)

0 引 言

稀疏線性方程組是混凝土細(xì)觀數(shù)值模擬[1]、數(shù)值天氣預(yù)報(bào)格點(diǎn)模式[2]、線性彈性問題[3]、油藏模擬[4]等許多科學(xué)與工程計(jì)算領(lǐng)域的核心計(jì)算模塊之一，求解時(shí)間在整個(gè)應(yīng)用問題中占很大比重，其求解效率是影響整個(gè)應(yīng)用問題數(shù)值模擬效率的關(guān)鍵環(huán)節(jié)之一。正因?yàn)榫哂袠O端重要性，因此，直到最近，還廣受關(guān)注并有大量研究[5-6]。Krylov子空間迭代[7-8]是求解這種稀疏線性方程組最有效而健壯的一般性方法之一，當(dāng)系數(shù)矩陣對(duì)稱正定時(shí)，其中共軛斜量(CG)法以計(jì)算經(jīng)濟(jì)、理論上較快的收斂速度應(yīng)用最為廣泛。但這種迭代法的實(shí)際收斂速度與系數(shù)矩陣的特征值分布有關(guān)，分布越集中，收斂速度越快。

多重網(wǎng)格算法[9-10]是求解稀疏線性方程組的另一類有效方法，其顯著特點(diǎn)是通過將光滑與校正有機(jī)結(jié)合，具有潛在的快速收斂性。但這類方法健壯性一般不理想，因此，將其用作預(yù)條件，來加速Krylov子空間迭代法，是將兩者優(yōu)勢有機(jī)結(jié)合的最佳選擇。在針對(duì)一般稀疏線性方程組的代數(shù)多重網(wǎng)格算法中，聚集型代數(shù)多重網(wǎng)格算法以其計(jì)算與存儲(chǔ)上的經(jīng)濟(jì)性、可控性廣受關(guān)注。對(duì)這種多重網(wǎng)格算法，其核心是網(wǎng)格層次結(jié)構(gòu)的選取與確定方法[11]，經(jīng)典的聚集算法中，最有效的是基于強(qiáng)連接的算法及其變種[12-13]，但這些算法主要基于系數(shù)矩陣鄰接圖的局部特性，因而具有容易陷入局部最優(yōu)的潛在缺陷，且不易于進(jìn)行并行實(shí)現(xiàn)。文獻(xiàn)[14]對(duì)多種聚集算法進(jìn)行了比較，發(fā)現(xiàn)基于強(qiáng)連接的算法及其變種一般具有優(yōu)勢。文獻(xiàn)[15]中提出了一種以圖分割為基礎(chǔ)的自頂向下型聚集構(gòu)造算法，有效克服了這些問題。當(dāng)已知線性方程組每個(gè)未知量對(duì)應(yīng)的坐標(biāo)時(shí)，利用坐標(biāo)分割可以更快速地進(jìn)行鄰接圖的分割，從而減少設(shè)置階段所需要的時(shí)間。已有實(shí)驗(yàn)結(jié)果表明，這種算法相對(duì)于基于強(qiáng)連接的經(jīng)典聚集算法，無論在計(jì)算效率，還是參數(shù)敏感性與健壯性方面，都具有明顯優(yōu)勢。

1 基于坐標(biāo)分割的多重網(wǎng)格預(yù)條件

考慮下述稀疏線性方程組：

Ax=b

(1)

式中：A是已知的對(duì)稱正定稀疏矩陣，b為已知右端向量，x為待求的未知解向量。對(duì)式(1)，可以采用預(yù)條件CG法進(jìn)行求解。采用預(yù)條件的目的是希望將式(1)化為另一個(gè)具有同解、且系數(shù)矩陣特征值分布更集中的線性方程組，同時(shí)應(yīng)使得其中牽涉到的輔助計(jì)算量比較少。從這個(gè)意義上說，對(duì)A左乘或右乘一個(gè)近似逆矩陣是可行之道，這在輔助計(jì)算時(shí)又相當(dāng)于近似求解式(1)，因此，任何近似求解式(1)的方法都可以作為預(yù)條件。多重網(wǎng)格型算法具有潛在最優(yōu)收斂性，自然是用于式(1)的有效預(yù)條件選擇之一。

聚集型代數(shù)多重網(wǎng)格具有存儲(chǔ)與計(jì)算上的經(jīng)濟(jì)性，且易于控制，因而廣受關(guān)注。聚集構(gòu)造是這種算法的核心問題之一。在文獻(xiàn)[8]中提出了一種自頂向下構(gòu)建聚集與網(wǎng)格層次結(jié)構(gòu)的方法，該方法的核心思想是從對(duì)應(yīng)于A的鄰接圖G出發(fā)，先利用圖分割算法，在要求子圖之間連接邊的邊權(quán)之和最小的約束條件下，將G分為多個(gè)子圖，使得各個(gè)子圖中節(jié)點(diǎn)個(gè)數(shù)基本相當(dāng)。之后，可以直接將每個(gè)子圖中的節(jié)點(diǎn)聚集為一個(gè)點(diǎn)，從而形成最粗的網(wǎng)格層。

為構(gòu)建其他網(wǎng)格層，可以先從最粗層中每個(gè)點(diǎn)對(duì)應(yīng)的子圖出發(fā)，將每個(gè)子圖又按類似方式分為多個(gè)子圖。這樣，如果假設(shè)每次分為p個(gè)子圖，則從最粗層開始，每細(xì)化一層，該層中對(duì)應(yīng)的子圖個(gè)數(shù)都增加到上一層的p倍。之后，每一層中的每個(gè)子圖就自然形成一個(gè)聚集。按照所要求的條件，在進(jìn)行圖分割時(shí)，將使得子圖之間連接邊的邊權(quán)之和最小，因而可以一定程度上確保子圖之間的依賴關(guān)系較小。這與強(qiáng)連接的意義相近，將連接關(guān)系較強(qiáng)的節(jié)點(diǎn)盡量聚集在一起。但在此時(shí)以自頂向下的方式進(jìn)行的，這與基于強(qiáng)連接的經(jīng)典聚集算法具有本質(zhì)上的不同，經(jīng)典聚集算法是從局部搜索出發(fā)，對(duì)強(qiáng)連接的節(jié)點(diǎn)進(jìn)行聚集。

由于基于圖分割的聚集算法需要反復(fù)用到圖分割算法，如果完全按照使得子圖之間的連接邊權(quán)之和最小的要求來進(jìn)行圖分割，將會(huì)導(dǎo)致設(shè)置階段時(shí)間相當(dāng)慢長。因此，只能采用啟發(fā)式算法。對(duì)一般的圖，可以采用MeTis 等軟件包進(jìn)行圖的分割。但已有實(shí)踐表明，采用這種圖分割軟件進(jìn)行分割時(shí)，可能導(dǎo)致某些子圖不連通。而對(duì)不連通結(jié)點(diǎn)進(jìn)行聚集，從根本上與強(qiáng)連接的思想相違背，本質(zhì)上將會(huì)導(dǎo)致多重網(wǎng)格算法性能的下降。因此，如何進(jìn)行既能確保高效實(shí)現(xiàn)，又能確保子圖連通性的圖分割，是確保這種聚集型代數(shù)多重網(wǎng)格算法有效性的核心問題。

幸運(yùn)的是，目前所求解的很多稀疏線性方程組，其來源多是大規(guī)?？茖W(xué)與工程計(jì)算。有的來源于偏微分方程數(shù)值求解，例如數(shù)值天氣預(yù)報(bào)的格點(diǎn)模式等應(yīng)用問題，是基于空間網(wǎng)格點(diǎn)進(jìn)行離散所得到的。有的雖然不是來源于微分方程，但也基于空間網(wǎng)格點(diǎn)進(jìn)行離散計(jì)算，例如混凝土細(xì)觀數(shù)值模擬等結(jié)構(gòu)力學(xué)中所得到的稀疏線性方程組就是典型例子。對(duì)這些稀疏線性方程組，每個(gè)未知量所對(duì)應(yīng)的空間坐標(biāo)已知，如果充分利用這種已知信息，可以進(jìn)行既高效又確保連通性的圖分割算法設(shè)計(jì)。

2 三種基于坐標(biāo)的高效圖分割算法

考慮無向圖G=(V,E)，其中V={1,2,…,n}，E為連接邊集合，假設(shè)V中結(jié)點(diǎn)i對(duì)應(yīng)的坐標(biāo)為(xi,yi,zi)，現(xiàn)在需要將其分割為p個(gè)子圖。本節(jié)對(duì)這一問題，給出正方分割、最小界面分割與逐步單向分割等三種分割算法及其高效實(shí)現(xiàn)。

正方分割算法以每個(gè)子圖接近于正方體或正方形的方式進(jìn)行分割。在該算法中，先找出不同x、y、z坐標(biāo)的個(gè)數(shù)，分別記為nxs、nys、nzs，并記下每個(gè)x、y、z坐標(biāo)出現(xiàn)的次數(shù)，第i個(gè)x坐標(biāo)出現(xiàn)的次數(shù)記為ixs(i)，第j個(gè)y坐標(biāo)出現(xiàn)的次數(shù)記為iys(j)，第k個(gè)z坐標(biāo)出現(xiàn)的次數(shù)記為izs(k)。為進(jìn)行高效實(shí)現(xiàn)，先采用歸并排序?qū)⑺薪Y(jié)點(diǎn)按照x坐標(biāo)優(yōu)先、再y坐標(biāo)優(yōu)先、后z坐標(biāo)優(yōu)先的順序，與坐標(biāo)從小到大的順序進(jìn)行排序，每個(gè)結(jié)點(diǎn)對(duì)應(yīng)的坐標(biāo)相應(yīng)進(jìn)行調(diào)整。之后，分別從x、y、z三個(gè)方向，通過二分查找的方式，從已排序序列中找出每個(gè)方向上不同坐標(biāo)的個(gè)數(shù)，并記錄下該方向每個(gè)坐標(biāo)出現(xiàn)的次數(shù)。以y方向?yàn)槔?，?duì)排序后的結(jié)點(diǎn)依次進(jìn)行掃描，先置nys=1，并將iys(1)置為第一個(gè)結(jié)點(diǎn)的y坐標(biāo)，之后，對(duì)每個(gè)新結(jié)點(diǎn)，查找其對(duì)應(yīng)y坐標(biāo)在iys(1)到iys(nys)是否存在，如果不存在，則nys加1，并置iys(nys)等于該y坐標(biāo)。

在進(jìn)行圖分割時(shí)，按盡量使得每個(gè)子圖中的結(jié)點(diǎn)形成近似正方體或正方形的方式來進(jìn)行組織。對(duì)三維問題，即當(dāng)nzs=1時(shí)，按子圖接近于正方形的方式進(jìn)行分割，否則按子圖接近于正方體的方式進(jìn)行分割。這樣可以最大程度上減少各子圖之間的連接邊數(shù)量。設(shè)將G分為p=px×py×pz個(gè)子圖，在這種方式下，即盡量使得nxs:px=nys:py=nzs:pz，并盡量使得各方向上的坐標(biāo)在該方向的分割中進(jìn)行平均分布。

最小界面分割算法遍歷所有可能的分割，并以每個(gè)子圖表面積或周長之和最小的方式進(jìn)行實(shí)際分割。在該算法中，先采用與正方分割類似的方式，找出不同坐標(biāo)的個(gè)數(shù)，以及每個(gè)坐標(biāo)出現(xiàn)的次數(shù)。再確定在分割數(shù)各種形式的分解p=px×py×pz中，哪一種最優(yōu)。這里確定分割組織的方式與正方分割不同，是通過因數(shù)分解，使得每個(gè)子圖表面積(三維情形)或周長(二維情形)總和最小。具體地，先將p進(jìn)行因素分解，由此可以確定px、py、pz的各種不同組合形式，之后，對(duì)每種組合形式，在二維情況下可以計(jì)算每個(gè)子圖的周長之和，在三維情況下，可以計(jì)算每個(gè)子圖的表面積之和，最后采用使得總和最小的組合方式作為最終的分割組織形式。

逐步單向分割算法基于分割數(shù)的素因子分解，按素因子從大到小的順序，每次沿不同坐標(biāo)數(shù)最大的方向進(jìn)行分割，直到所有素因子遍歷完為止。在該算法中，先對(duì)p進(jìn)行因數(shù)分解，并對(duì)各素因子按從大到小的順序進(jìn)行排序，之后，對(duì)G先按第一個(gè)素因子進(jìn)行分割，并對(duì)分割所得子圖，再按第二個(gè)素因子進(jìn)行第二層分割，如此繼續(xù)，直到所有素因子訪問完畢為止。在對(duì)某個(gè)子圖按照某個(gè)素因子進(jìn)行分割時(shí)，按照與正方分割中相同的方式找出每個(gè)方向上的不同坐標(biāo)個(gè)數(shù)，之后，分割沿先坐標(biāo)個(gè)數(shù)最大、再坐標(biāo)個(gè)數(shù)其次、最后坐標(biāo)個(gè)數(shù)最小的優(yōu)先順序進(jìn)行。因此，在進(jìn)行每次分割之前，也需要進(jìn)行排序，這依然采用歸并排序來實(shí)現(xiàn)。這里，對(duì)p的所有素因子按從大到小的順序進(jìn)行排序，并依次進(jìn)行分割，是為了確保最后分割所得每個(gè)子圖中擁有盡可能相等的結(jié)點(diǎn)個(gè)數(shù)，這在多重網(wǎng)格算法構(gòu)造中，有利于每層中每個(gè)子圖中的結(jié)點(diǎn)數(shù)最大程度地相等，從而最終有利于確保多重網(wǎng)格算法的有效性。

3 數(shù)值實(shí)驗(yàn)

本文實(shí)驗(yàn)平臺(tái)中的處理器為Intel(R) Xeon(R) CPU E5-4640 0 @ 2.40 GHz(cache 20 480 KB)，操作系統(tǒng)為Red Hat Linux 2.6.32-431.el6.x86_64，編譯器為Intel FORTRAN Version 15.0.0.090。由于所求解的線性方程組均對(duì)稱正定，因此全部采用PCG迭代進(jìn)行求解。初始迭代向量選為全0向量，并在當(dāng)前殘向量的2范數(shù)與初始?xì)埾蛄康?范數(shù)之比小于1E-10時(shí)，迭代終止。在多重網(wǎng)格中，最粗層線性方程組的規(guī)?？刂茷?00階，且對(duì)該線性方程組，采用無填充的不完全LU分解預(yù)條件PCG進(jìn)行求解，終止條件與外迭代相同。

實(shí)驗(yàn)針對(duì)6個(gè)稀疏線性方程組進(jìn)行，包括Lin21、Lin22、Lin23、Lin31、Lin32與Lin32。所有線性方程組都采用有限差分離散得到，其中Lin21、Lin22與Lin23從離散二維偏微分方程的Dirichlet邊值問題得到：

-a1?2u/?x2-a2?2u/?y2=f

(1)

定義域?yàn)?0,1)×(0,1)，函數(shù)f與邊界條件從真解u=1得到。在每個(gè)維度上有n+2個(gè)點(diǎn)，其中n=1 024，且對(duì)任意連續(xù)函數(shù)u，將值u(xi,yj)定義為ui,j，其中xi與yj分別沿x與y方向上一致分布。線性方程組Lin21對(duì)應(yīng)于a1=a2=1，且采用下述離散形式：

-ui,j-1-ui-1,j+4ui,j-ui+1,j-ui,j+1=h2fi,j

(2)

線性方程組Lin22對(duì)應(yīng)于a1=1、a2=100，且采用下述離散形式:

-100ui,j-1-ui-1,j+202ui,j-ui+1,j-100ui,j+1=h2fi,j

(3)

線性方程組Lin23對(duì)應(yīng)于a1=a2=1，且采用的離散形式為:

-ui-1,j-1-4ui,j-1-ui+1,j-1-4ui-1,j+20ui,j-4ui+1,j-

ui-1,j+1-4ui,j+1-ui+1,j+1=h2fi,j

(4)

線性方程組Lin31、Lin32與Lin33從離散三維偏微分方程的Dirichlet邊值問題得到：

-a1?2u/?x2-a2?2u/?y2-a3?2u/?y2=f

(5)

定義域?yàn)?0,1)×(0,1)×(0,1)，函數(shù)f與邊界條件從真解u=1得到。在每個(gè)維度上有n+2個(gè)點(diǎn)，其中n=128，且對(duì)任意連續(xù)函數(shù)u，將值u(xi,yj,zk)定義為ui,j,k，其中xi、yj與zk分別沿x、y與z方向上一致分布。線性方程組Lin31對(duì)應(yīng)于a1=a2=a3=1，且采用下述離散形式：

-ui-1,j,k-ui,j-1,k-ui,j,k-1+6ui,j,k-ui+1,j,k-

ui,j+1,k-ui,j,k+1=h2fi,j,k

(6)

線性方程組Lin32對(duì)應(yīng)于a1=1、a2=10與a3=100，且采用下述離散形式：

-ui-1,j,k-10ui,j-1,k-100ui,j,k-1+222ui,j,k-ui+1,j,k

-10ui,j+1,k-100ui,j,k+1=h2fi,j,k

(7)

線性方程組Lin33對(duì)應(yīng)于a1=a2=a3=1，且采用的離散形式為：

27ui,j,k=h2fi,j,k

(8)

本文所有實(shí)驗(yàn)結(jié)果的時(shí)間單位均為秒，在所有表中，p的含義如第1節(jié)中所述，V、W分別表示采用V-與W-循環(huán)的多重網(wǎng)格，K25與K35分別表示采用參數(shù)t=0.25與t=0.35的K-循環(huán)網(wǎng)格。所有表中SQ表示正方分割，MI表示最小界面分割，ST表示逐步單向分割。

求解線性方程組Lin21所用的迭代時(shí)間如表1所示，在此測試用例中，系數(shù)矩陣每行中的非零元個(gè)數(shù)為5，且各向同性。從表1可以看到，當(dāng)采用Jacobi光滑時(shí)，絕大部分情況下，ST優(yōu)勢明顯。當(dāng)采用Gauss-Seidel光滑時(shí)，MI更具有優(yōu)勢。

求解線性方程組Lin22的時(shí)間結(jié)果如表2所示，在此測試用例中，系數(shù)矩陣每行中的非零元個(gè)數(shù)依然為5，但該問題具有很強(qiáng)的各向異性。在表2中，“-”表示在1 000次迭代內(nèi)未能達(dá)到迭代收斂條件。從表中可見，采用Jacobi光滑時(shí)，依然是絕大多數(shù)情況下，算法ST更具有優(yōu)勢。在采用Gauss-Seidel光滑時(shí)，對(duì)V-與W-循環(huán)，在p較小時(shí)，一般MI最優(yōu)，當(dāng)p取8時(shí)，SQ具有明顯優(yōu)勢；對(duì)K-循環(huán)，一般ST更有優(yōu)勢。特別值得注意的是，在p較小時(shí)，對(duì)K-循環(huán)，SQ與MI都存在1 000次迭代內(nèi)未達(dá)到收斂標(biāo)準(zhǔn)的情況，但ST不存在這個(gè)問題，這說明ST具有更好的健壯性。

表2 求解線性方程組Lin22的時(shí)間 s

求解線性方程組Lin23的時(shí)間結(jié)果如表3所示，在此測試中，系數(shù)矩陣每行中的非零元個(gè)數(shù)為9，該問題具有一定的各向異性。對(duì)迭代階段所用時(shí)間而言，對(duì)W-與V-循環(huán)，大部分情況下，MI更好。對(duì)K-循環(huán)而言，特別是在p較大時(shí)，ST更有優(yōu)勢。

表3 求解線性方程組Lin23的時(shí)間 s

求解線性方程組Lin31、Lin32與Lin33的結(jié)果分別如表4、5與6所示。從表4可見，對(duì)系數(shù)矩陣每行中只有7個(gè)非零元的三維同構(gòu)問題Lin31，當(dāng)采用Jacobi光滑時(shí)，ST一致性地具有優(yōu)勢。當(dāng)采用Gauss-Seidel光滑時(shí)，MI一致性地具有優(yōu)勢。

從表5可見，對(duì)每行中非零元個(gè)數(shù)較少的各向異性問題Lin32，采用Jacobi光滑時(shí)，一般ST具有優(yōu)勢。當(dāng)采用Gauss-Seidel光滑時(shí)，對(duì)V-與W-循環(huán)，一般SQ更優(yōu)，而對(duì)K循環(huán)，一般ST更優(yōu)。從表6可見，對(duì)系數(shù)矩陣每行中具有很多非零元的各項(xiàng)同性問題Lin33，一般采用MI更有優(yōu)勢，采用SQ時(shí)，效果也相差不大。

表4 求解線性方程組Lin31的時(shí)間 s

表5 求解線性方程組Lin32的時(shí)間 s

表6 求解線性方程組Lin33的時(shí)間 s

4 結(jié) 語

本文針對(duì)基于坐標(biāo)分割的聚集型代數(shù)多重網(wǎng)格預(yù)條件，提出了正方分割、最小界面分割與逐步單向分割等三種坐標(biāo)分割的算法，對(duì)其進(jìn)行了高效實(shí)現(xiàn)，并通過從二維與三維模型偏微分方程離散所得同構(gòu)與異構(gòu)、每行中非零元個(gè)數(shù)較少與較多等多種類型的問題，對(duì)其進(jìn)行了對(duì)比實(shí)驗(yàn)與分析研究。實(shí)驗(yàn)結(jié)果表明，當(dāng)采用Jacobi光滑時(shí)，除了系數(shù)矩陣每行中非零元個(gè)數(shù)較多的情況之外，一般采用逐步單向分割更有優(yōu)勢。當(dāng)采用Gauss-Seidel光滑時(shí)，一般采用最小界面分割算法更優(yōu)。此外，當(dāng)系數(shù)矩陣每行中的非零元個(gè)數(shù)較多時(shí)，一般采用最小界面分割算法更有優(yōu)勢。對(duì)K-循環(huán)多重網(wǎng)格，一般采用逐步單向分割算法也更具有優(yōu)勢，而且，這種分割的健壯性更好，即使對(duì)強(qiáng)各向異性的問題，最終所得多重網(wǎng)格預(yù)條件迭代的效果也相當(dāng)穩(wěn)健。需要說明的是，本文主要針對(duì)模型偏微分方程離散所得稀疏線性方程組進(jìn)行了算法的測試與驗(yàn)證，雖然既針對(duì)各項(xiàng)同性與各項(xiàng)異性的問題進(jìn)行了實(shí)驗(yàn)驗(yàn)證，又針對(duì)每行中非零元個(gè)數(shù)較少與較多的情況進(jìn)行了測試分析，但對(duì)實(shí)際大規(guī)?？茖W(xué)與工程計(jì)算領(lǐng)域中的需要求解的稀疏線性方程組，一方面可能所得系數(shù)矩陣每行中的非零元個(gè)數(shù)更多，且每個(gè)坐標(biāo)可能對(duì)應(yīng)于多個(gè)自由度，如對(duì)混凝土細(xì)觀數(shù)值模擬，采用六面體單元離散時(shí)，對(duì)應(yīng)系數(shù)矩陣中每行含有約81個(gè)非零元素，且每個(gè)坐標(biāo)通常對(duì)應(yīng)于3個(gè)應(yīng)變分量。同時(shí)，這里模型問題對(duì)應(yīng)的系數(shù)矩陣都是M矩陣，而實(shí)際應(yīng)用問題中對(duì)應(yīng)的系數(shù)矩陣有時(shí)并不具有這種性質(zhì)。雖然本文給出的算法對(duì)待求解的稀疏線性方程組并沒有上述限制條件，但實(shí)際應(yīng)用效果有待實(shí)驗(yàn)驗(yàn)證。因此，下一步將就本文給出的算法，對(duì)混凝土細(xì)觀數(shù)值模擬等實(shí)際應(yīng)用領(lǐng)域中的稀疏線性方程組，進(jìn)行分析驗(yàn)證。