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

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

0 引 言

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

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

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

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

Ax=b

(1)

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

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

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

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

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

2 三種基于坐標的高效圖分割算法

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

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

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

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

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

3 數(shù)值實驗

本文實驗平臺中的處理器為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。由于所求解的線性方程組均對稱正定，因此全部采用PCG迭代進行求解。初始迭代向量選為全0向量，并在當前殘向量的2范數(shù)與初始殘向量的2范數(shù)之比小于1E-10時，迭代終止。在多重網(wǎng)格中，最粗層線性方程組的規(guī)模控制為100階，且對該線性方程組，采用無填充的不完全LU分解預條件PCG進行求解，終止條件與外迭代相同。

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

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

(1)

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

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

(2)

線性方程組Lin22對應于a1=1、a2=100，且采用下述離散形式:

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

(3)

線性方程組Lin23對應于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)

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

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

(8)

本文所有實驗結(jié)果的時間單位均為秒，在所有表中，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所用的迭代時間如表1所示，在此測試用例中，系數(shù)矩陣每行中的非零元個數(shù)為5，且各向同性。從表1可以看到，當采用Jacobi光滑時，絕大部分情況下，ST優(yōu)勢明顯。當采用Gauss-Seidel光滑時，MI更具有優(yōu)勢。

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

表2 求解線性方程組Lin22的時間 s

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

表3 求解線性方程組Lin23的時間 s

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

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

表4 求解線性方程組Lin31的時間 s

表5 求解線性方程組Lin32的時間 s

表6 求解線性方程組Lin33的時間 s

4 結(jié) 語

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