王如云，于音弦，臧振濤，曹 迪，陳 林

(1.河海大學(xué)港口海岸與近海工程學(xué)院，江蘇 南京210098;2.河海大學(xué) 力學(xué)與材料學(xué)院，江蘇南京210098)

把地球流體、固體等大型科學(xué)計算問題的計算域劃分成較小的子區(qū)域，并映射到各計算節(jié)點，是并行計算中的一個重要過程，稱為區(qū)域剖分。尋求各計算節(jié)點上負(fù)載平衡的計算網(wǎng)格優(yōu)化劃分方法，稱為區(qū)域剖分法(Domain Decomposition Method)，其結(jié)果將直接影響著并行計算的并行效率。因此，區(qū)域剖分法一直是國內(nèi)外許多高性能計算科研工作者的研究重點之一。

由于能夠?qū)?fù)雜外形進(jìn)行準(zhǔn)確、有效地擬合，無結(jié)構(gòu)網(wǎng)格在有限單元法、有限差分法、有限體積法等數(shù)值方法中都有廣泛應(yīng)用。目前，基于無結(jié)構(gòu)網(wǎng)格進(jìn)行并行計算時，常用的區(qū)域剖分方法有:陣面推進(jìn)法(Advancing Front Technique)［1-2］、圖劃分法(Graph Partitioning Methods)［3］等，其中圖劃分法又包含遞歸譜對剖分法(Recursive Spectral Bisection，RSB)［4-6］、 多 層 次 算 法 (Multilevel Algorithms)［7-8］、Multilevel Recursive Spectral Bisection(MRSB)方法［9］、超圖剖分法(Hypergraph Partitioning)［9-10］等。

陣面推進(jìn)法的優(yōu)勢在于其操作簡單，且易于實現(xiàn)，但其存在計算量大、實施過程較復(fù)雜等問題;遞歸譜對剖分法具有分割量近似最小，使得處理器間信息交換量接近最小的優(yōu)點，但由于該算法所形成的Lagrangian矩陣大小，與總的單元數(shù)量成平方關(guān)系，這就消耗了大量內(nèi)存，同時該算法需要求解Fiedler特征向量也需要大量的計算時間;多層次算法大致可分為粗化、分區(qū)與優(yōu)化3個階段，優(yōu)化階段一般采用 Kernighan-Lin/Fiduccia-Mattheyses(KL/FM)［11］類方法，但該算法實現(xiàn)步驟比較復(fù)雜，KL/FM的使用會加大整體耗時量;MRSB方法結(jié)合了多層次算法與RSB方法，是目前國際上比較通用的方法，該方法具有處理器的開銷較小、計算速度較快等優(yōu)點，但由于有限元模型的幾何拓?fù)湫畔⒃诖只^程中會出現(xiàn)部分丟失現(xiàn)象，造成各個分區(qū)的元素不均衡和邊界過長等不良結(jié)果;超圖剖分方法能夠非常精確地表示通信量、具有外部通信量較小等優(yōu)點，但它存在費用開銷大于其它圖剖分方法的缺點。

根據(jù)計算節(jié)點個數(shù)情況，把所有單元按編號順序，均勻分配到各個計算節(jié)點是一種較為簡單易于實現(xiàn)的區(qū)域剖分方案。但是直接對由自動網(wǎng)格剖分程序生成的單元，按單元編號進(jìn)行區(qū)域剖分其并行效率一般很低。通過研究分析，文中給出單元鄰接矩陣帶寬與外部通信量間的關(guān)系，發(fā)現(xiàn)通過減小帶寬可以達(dá)到減少外部通信量，提高并行效率的目的。減小矩陣帶寬可以通過給網(wǎng)格單元進(jìn)行編號優(yōu)化來實現(xiàn)，其方法較多。如Akhras和DhattG于1976年提出一種通過優(yōu)化網(wǎng)格節(jié)點編號來減小高階稀疏矩陣帶寬的算法［12-13］，簡稱AD(Akhras和Dhatt G)算法就是其中一種。AD算法簡便且利于并行處理。

基于AD算法思想，文中提出一種區(qū)域剖分方法。首先利用AD算法給無結(jié)構(gòu)網(wǎng)格優(yōu)化編號以減少鄰接矩陣的帶寬，然后根據(jù)優(yōu)化后的單元編號順序和單元屬性，將計算單元均勻地映射到各個計算節(jié)點上。通過算例與其它剖分方法比較，結(jié)果表明該方法能夠?qū)崿F(xiàn)各處理器上負(fù)載的更為優(yōu)化的平衡，對處理大型數(shù)據(jù)具有可行性和有效性。

1 單元鄰接矩陣帶寬與并行計算效率間的關(guān)系

在通常的流體、固體等數(shù)值模擬計算中，當(dāng)前單元物理量的計算需要調(diào)用其周邊單元的數(shù)據(jù)。當(dāng)將網(wǎng)格進(jìn)行區(qū)域剖分后，會造成部分需要數(shù)據(jù)傳輸?shù)膯卧辉谕挥嬎愎?jié)點上，這就需要進(jìn)行節(jié)點間的外部通信。適當(dāng)減少外部通信量，可以提高并行計算效率。

為研究方便，將整個計算域的網(wǎng)格單元進(jìn)行編號:1，2，…，N。先給出如下 3 個定義:

1)相關(guān)單元:在當(dāng)前單元上進(jìn)行數(shù)值計算時，需要調(diào)用到其它單元上的有關(guān)數(shù)據(jù)，這些被調(diào)用數(shù)據(jù)的單元稱為當(dāng)前單元的相關(guān)單元。若第i個單元是第j個單元的相關(guān)單元，就一般的計算格式而言，第j個單元也是第i個單元的相關(guān)單元，稱單元i和單元j相關(guān)。

引進(jìn)

一般數(shù)值計算時會調(diào)用當(dāng)前單元上的數(shù)據(jù)，因此取aii=1。若記，則該值與單元的具體編號無關(guān)，是鄰接矩陣的一個不變特征量。

2)單元鄰接矩陣:矩陣A=(aij)N×N稱為單元鄰接矩陣，滿足A=AT。

3)按單元序號區(qū)域剖分法(簡稱序號法):假設(shè)計算節(jié)點總數(shù)為r，現(xiàn)把N(N=Kr)個單元上的計算任務(wù)按單元編號順序，均勻分配到各個計算節(jié)點，其中第I(I=1，2，…，r)個計算節(jié)點承擔(dān)單元編號(I－1)K+1，…，IK總共K個單元的計算任務(wù)。

在進(jìn)行按單元編號順序區(qū)域剖分法后，此時單元鄰接矩陣分解成r階分塊矩陣

其中

AI，J=(a(I－1)K+i，(J－1)K+j)K×K，I，J=1，2，…，r顯然，第I個計算節(jié)點與第J個計算節(jié)點間需要交換信息的單元總數(shù) EI，J，就是分塊矩陣 AI，J中的元素和，即

第I個計算節(jié)點與第J個計算節(jié)點間有關(guān)通信量與 EI，J成正比，比例系數(shù)為常數(shù)。當(dāng) I=J 時，EI，J表示第I個計算節(jié)點內(nèi)部通信的單元總數(shù);當(dāng)I≠J時，EI，J表示第I個計算節(jié)點與第J個計算節(jié)點間外部通信的單元總數(shù)。因此，第I個計算節(jié)點與其它所有計算節(jié)點(J=1，…，I－ 1，I+1，…，r)間的外部通信的單元總數(shù)為

假設(shè)鄰接矩陣的半帶寬為L(見圖1)，則位于帶狀區(qū)域外的元素，即當(dāng)|i－j|＞L時有ai，j=0;而在位于帶狀區(qū)域內(nèi)的元素，即當(dāng)|i－j|≤L時ai，j取值可能等于1。假設(shè)帶狀區(qū)域內(nèi)，除主對角線上元素取值為1外，其它位置的元素取值為1的機(jī)會相等，利用特征量m0可導(dǎo)出外部通信的單元總數(shù)

的計算公式，分兩種情況進(jìn)行研究。

1)當(dāng)K≤L+1≤N時

2)當(dāng)1≤L+1≤K時

綜合以上兩種情況可以進(jìn)一步得到結(jié)論:

2)從鄰接矩陣帶狀結(jié)構(gòu)可見，第1個和最后一個計算節(jié)點需要進(jìn)行外部通信的單元總數(shù)為，其它計算節(jié)點需要進(jìn)行外部通信的單元總數(shù)為

3)由于計算節(jié)點間的外部通信比計算節(jié)點內(nèi)部通信耗時多，因此各個計算節(jié)點上需要進(jìn)行外部通信單元總數(shù)減少(轉(zhuǎn)成內(nèi)部通信)，亦即減少了外部通信量，從而可以提高各個計算節(jié)點的計算效率。

因此，減小單元鄰接矩陣的帶寬，可以達(dá)到減小外部通信量，從而達(dá)到提高并行計算效率的目的。

圖1 鄰接矩陣示意Fig.1 Connectivity matrix sketch map

2 減小鄰接矩陣帶寬的AD算法

2.1 AD算法簡介

AD算法［12］根據(jù)如下3個指標(biāo)的大小決定節(jié)點新的編號:(1)節(jié)點和:若以當(dāng)前節(jié)點為頂點的單元有m個，對各單元所有頂點的編號求和，然后求出m個單元頂點編號總和，即節(jié)點和;(2)節(jié)點商:節(jié)點和除以m得到的商稱為節(jié)點商;(3)最大最小編號和:若以當(dāng)前節(jié)點為頂點的單元中，最大的頂點編號與最小的頂點編號之和。

“三農(nóng)”主持人首先應(yīng)該是對農(nóng)村、農(nóng)民、農(nóng)業(yè)有感情的人。熟悉廣播節(jié)目的聽眾都知道，并不是聲音甜美就是親切，并不是語言樸實就是貼心，親切和貼心都源自節(jié)目的編播人員對聽眾的情感。

AD算法分為2個階段:第1階段使用節(jié)點商和最大最小編號和作為排序指標(biāo)，按升序排列，直到整體的迭代步數(shù)大于最大迭代次數(shù)，或者帶寬序列中相鄰3次數(shù)值相同時轉(zhuǎn)入第2階段;第2階段使用節(jié)點和及節(jié)點商作為排序指標(biāo)，按升序排列，直到滿足類似于第1階段的終止條件。

2.2 減小單元鄰接矩陣帶寬的AD算法

最初設(shè)計AD算法主要是為了減少矩陣儲存量，其手段是通過對節(jié)點編號進(jìn)行優(yōu)化，使得矩陣帶寬減小而達(dá)成目的。通過前面對單元鄰接矩陣帶寬與并行計算效率間的關(guān)系分析得到，減小單元鄰接矩陣的帶寬，可以達(dá)到減小外部通信量，從而提高并行效率的目的。因此，為了提高并行計算效率，基于AD算法思想，通過引進(jìn)類似節(jié)點商的單元商對單元編號進(jìn)行優(yōu)化，使得矩陣帶寬減小，從而設(shè)計了一種減小單元鄰接矩陣帶寬的AD算法。

仿照優(yōu)化節(jié)點編號的AD算法，引出如下變量:單元和:若與當(dāng)前單元有通信關(guān)系的單元有m個，對這m個單元的編號求和，該和稱為單元和;單元商:當(dāng)前的單元和除以與當(dāng)前單元有通信關(guān)系的單元總數(shù)(m)得到的商稱為單元商。減小鄰接矩陣帶寬的AD算法的排序指標(biāo):AD算法的每階段均由兩個指標(biāo)進(jìn)行升序排列，存在序號振蕩現(xiàn)象［14］，在一定程度上降低了算法的效率。此外，2個階段之間也存在類似問題，即在第1階段排好的節(jié)點編號到第2階段時，由于指標(biāo)的變更，有可能局部次序被打亂，從而在一定程度上降低了算法的效率。為了避免AD算法存在的相互矛盾和震蕩問題，文中在此只采用單元商作為減小鄰接矩陣帶寬的AD算法的排序指標(biāo)。

減小鄰接矩陣帶寬的AD算法終止條件:首先引入如下3個變量:

1)逆序比值Wn:按單元序號檢查前后相鄰兩單元的單元商，如果后一單元商比前一單元商小，則稱為一個逆序，所有逆序的總數(shù)與單元總數(shù)的比值稱為逆序比值;

2)半帶寬L:令當(dāng)前單元與離它最遠(yuǎn)的相關(guān)單元的編號差為L0，所有單元中最大的L0則為半帶寬;

3)最大單元商差 ΔA:ΔA=max(|A(i)－A0(i)|)，i=1，2，…，n，A(i)，A0(i)分別表示本次迭代、前次迭代時單元i的單元商。

減小單元鄰接矩陣帶寬的AD算法如下:(1)根據(jù)數(shù)值計算方法需要調(diào)用其它單元數(shù)據(jù)的情況，建立單元鄰接矩陣;(2)根據(jù)單元鄰接矩陣，計算各單元的單元商;(3)根據(jù)單元商的大小對單元進(jìn)行重新排序，記錄新舊單元編號互相對應(yīng)的關(guān)系;(4)根據(jù)新舊單元編號的對應(yīng)關(guān)系，以及原單元鄰接矩陣的信息，得到新單元編號下的鄰接矩陣;(5)計算逆序比值、半帶寬以及ΔA，如果滿足算法的終止條件轉(zhuǎn)(6)，否則轉(zhuǎn)(2);(6)程序結(jié)束，輸出結(jié)果。

2.3 基于AD算法的區(qū)域剖分法

2.3.1 按單元序號和屬性的區(qū)域剖分法 按序號法(見1中定義)進(jìn)行區(qū)域剖分，雖然能控制各子區(qū)域單元數(shù)量近似相等，但由于不同屬性單元(如邊界單元、內(nèi)單元)需要的計算時間是不同的，為了平衡各計算節(jié)點上的負(fù)載，在單元編號順序基礎(chǔ)上，需進(jìn)一步考慮單元屬性的影響。假設(shè)網(wǎng)格總數(shù)為N，計算節(jié)點個數(shù)為r，網(wǎng)格單元屬性有m種，不同屬性的單元數(shù)量分別為 n1，n2，…，nm，在序號法的基礎(chǔ)上，給出一種改進(jìn)算法。該算法具體思想是，對于任意的第j種屬性單元(j=1，2，…，m)，按數(shù)量 kj=［nj/r］將該類型單元按順次分配到各計算節(jié)點，盡量使各計算節(jié)點上每種類型單元的數(shù)量接近相等，且各節(jié)點上的單元總數(shù)接近N/r，各節(jié)點上的單元數(shù)總和仍為N。該算法稱為按單元序號和屬性的區(qū)域剖分法，簡稱序號屬性法。

2.3.2 基于AD算法的區(qū)域剖分法 序號屬性法一般能得到較好的分裂結(jié)果，各計算節(jié)點上的執(zhí)行時間達(dá)到了大體均衡的效果。但該方法沒有考慮外部通信量的問題，因此使用該區(qū)域剖分方法產(chǎn)生的外部通信量可能較大，在一定程度上影響了并行效率。因此在進(jìn)行序號屬性法之前，引入減少單元鄰接矩陣帶寬的AD算法對單元編號進(jìn)行優(yōu)化，以期盡可能減小鄰接矩陣的帶寬，進(jìn)而減少外部通信量，達(dá)到并行效率的進(jìn)一步提高，這就是基于AD算法的區(qū)域剖分法(簡稱AD分裂法)思想。

具體算法由兩部分組成:(1)根據(jù)鄰接矩陣情況，利用減小單元鄰接矩陣帶寬的AD算法對單元編號進(jìn)行優(yōu)化;(2)利用序號屬性法進(jìn)行區(qū)域剖分。

3 算例

測試平臺為win 732位操作系統(tǒng)，2G內(nèi)存的奔騰E5400臺式機(jī)。算例以研究全球潮波運動為背景，基于球面二維淺水波方程的離散求解WENO方法，計算單元采用的是只具有陸地邊界，而無水邊界的球面二維無結(jié)構(gòu)網(wǎng)格單元。由于離散格式的需要，計算當(dāng)前單元信息時不僅要用到第1層相鄰單元的信息，同時還要用到第2層相鄰單元總共9個相鄰單元的信息(稱為大模版，見圖2)。

對于陸地邊界單元需要利用鏡像法構(gòu)造虛擬單元并計算虛擬單元的物理量的值，故其計算時間會遠(yuǎn)大于內(nèi)單元的計算時間，區(qū)域分裂時必須區(qū)別對待內(nèi)單元和邊界單元。

在該數(shù)值模擬方法中，內(nèi)單元物理量計算時間t1=3.187 5×10－6s;邊界單元物理量計算時間t2=3.203 125×10－4s;一個單元上物理量的內(nèi)部和外部通信時間分別為 t3=3.125 × 10－8s，t4=3.906 25 ×10－6s。子計算節(jié)點 i上的執(zhí)行時間為:，其中分別表示子計算節(jié)點i上的內(nèi)單元、外單元、內(nèi)部通信以及外部通信的單元數(shù)量。

以下2個算例涉及到的計算網(wǎng)格單元，都是借助于地表水模擬軟件(Surface Modeling System，SMS)生成的全球無開邊界的二維無結(jié)構(gòu)網(wǎng)格單元中抽取的部分區(qū)域上的網(wǎng)格單元。

算例1 抽取的部分區(qū)域上的網(wǎng)格單元數(shù)量為54 311個，區(qū)域示意如圖3所示。采用不同數(shù)量的計算節(jié)點，分別運用序號法、序號屬性法以及目前使用較廣的Metis軟件多層次K路圖分區(qū)法，與AD分裂法進(jìn)行比較，所得結(jié)果如表1所示。在計算節(jié)點數(shù)為4時，應(yīng)用不同的區(qū)域剖分方法，所得區(qū)域剖分后的效果如圖3～6所示。圖上4種灰度分別代表不同計算節(jié)點上的計算單元。

圖2 經(jīng)典計算單元模板Fig.2 Classic cellcom puting template

表1 算例1中不同數(shù)量計算節(jié)點下的區(qū)域剖分效果Tab.1 Dom ain decom position result of different processors in example1

圖3 海域島嶼分布Fig.3 Seaisland distribution map

圖4 4個計算節(jié)點Metis圖分區(qū)法區(qū)域剖分結(jié)果Fig.4 Result of Graph Partitioning method by metis in 4 processors

算例2 抽取的部分區(qū)域上的網(wǎng)格單元數(shù)量為588 755個，其網(wǎng)格單元數(shù)量約為算例1的10倍。采用不同數(shù)量的計算節(jié)點，分別運用序號法、序號屬性法以及目前使用較廣的Metis軟件多層次K路圖分區(qū)法，與AD分裂法進(jìn)行比較，所得加速比、并行效率和外通倍數(shù)的結(jié)果如表2所示。

圖5 4個計算節(jié)點序號屬性法區(qū)域剖分結(jié)果Fig.5 Result of the serial number and property method in 4 processors

圖6 4個計算節(jié)點AD分裂法區(qū)域剖分結(jié)果Fig.6 Result of AD domain decom position method in 4 processors

表2 算例2中不同數(shù)量計算節(jié)點下區(qū)域剖分結(jié)果Tab.2 Dom ain decom position result of different processors in exam p le 2

以上結(jié)果表明，將網(wǎng)格單元進(jìn)行序號法分裂時，雖然簡單迅速，但并行效率不高。序號屬性法得到的并行效率有提高，但外通信單元數(shù)量較多，不利于并行效率的提高。Metis軟件圖分區(qū)法的主要目的是使各個計算節(jié)點上的單元數(shù)量相等，且外通信數(shù)量盡量少，但是對于采用模擬間斷水流的高精度WENO數(shù)值模型，邊界單元由于需要構(gòu)建虛擬單元以及賦值，其計算時間遠(yuǎn)大于內(nèi)單元;Metis軟件沒有區(qū)別對待內(nèi)單元以及邊界單元，分裂結(jié)果中各個計算節(jié)點上內(nèi)單元與邊界單元個數(shù)不均衡，使得各個子計算節(jié)點執(zhí)行時間不均衡，導(dǎo)致并行效率不高。AD分裂法能保證各計算節(jié)點在負(fù)載平衡的基礎(chǔ)上，有效減少外部通信單元數(shù)量，該方法考慮了內(nèi)單元與邊界單元的差異，最后各個計算節(jié)點上的執(zhí)行時間大體均勻，從而得到更高的并行效率。從圖6可以看出，AD分裂法得到的圖像中各子區(qū)域內(nèi)單元相對集中，各計算節(jié)點邊界單元個數(shù)均衡。當(dāng)對大數(shù)據(jù)進(jìn)行分裂時，AD分裂法也可以得到較高的并行效率，這一定程度上說明該方法在處理大數(shù)據(jù)時具有可行性與有效性。

4 結(jié)語

在AD算法優(yōu)化編號的基礎(chǔ)上對網(wǎng)格單元利用序號屬性法進(jìn)行區(qū)域剖分，能得到較高的加速比與并行效率。通過對高達(dá)59萬個無結(jié)構(gòu)網(wǎng)格單元進(jìn)行區(qū)域剖分試驗，結(jié)果表明利用AD分裂法進(jìn)行區(qū)域剖分，確實能使各計算節(jié)點的負(fù)載較快速趨于平衡，得到高加速比與并行效率的區(qū)域剖分結(jié)果，因此該方法是一種可用于大數(shù)據(jù)量的網(wǎng)格單元進(jìn)行區(qū)域剖分的可行且有效的方法。

［1］Lohner Rainald.A parallel advancing front grid generation scheme［J］.International Journal for Numerical Methods in Engineering，2001，51(6):663-678.

［2］司海青，王同光，成娟.非結(jié)構(gòu)網(wǎng)格上Eluer方程的區(qū)域分裂算法及其并行計算［J］.空氣動力學(xué)學(xué)報，2006，24(1):102-108.SIHaiqing，WANG Tongguang，CHEN Juan.Domain decompositions and parallel algorithms to solve Euler equations on the unstructured grid［J］.Acta Aerodynamica Sinica，2006，24(1):102-108.(in Chinese)

［3］ Henrickson B，Kolda TG.Graph partitioningmodels for parallel computing［J］.Parallel Computing，2000，26(12):1519-1534.

［4］ Pothen A，Simon H，Liou K P.Partitioning sparsematrices with eigenvectors of graphs［J］.SIAM JMath Anal Appl，1990，11(3):430-452.

［5］ Simon H D.Partitioning of unstructured problems for parallel processing［J］.Comput Eng，1991，2(2/3):35-148.

［6］周春華.并行計算中一種非結(jié)構(gòu)網(wǎng)格分割方法［J］.航空學(xué)報，2004，25(3):229-232.ZHOU Chunhua.A method of non-structured mesh partition for parallel computation［J］.Acta Aeronautica et Astronautica Sinica，2004，25(3):229-232.(in Chinese)

［7］ Walshaw C，Cross M.Mesh partitioning:a multilevel balancingand refinement algorithm［J］.SIAM Journal on Scientific Computing，2000，22(1):63-80.

［8］ Karypis G，Kumar V.A fast and high qualitymultilevel scheme for partitioning irregular graphs［J］.SIAM Journal on Scientific Computing，1998，20(1):359-392.

［9］Catalyurek U V，Boman E G，Heaphy R，et al.Hypergraph based dynamic load balancing for adaptive scientific computations［C］//Proc of IEEE International Conference on High Performance Computing.Long Beach:IEEE Press，2007:367-380.

［10］ Trifunovi A，KnottenbeltW J.Parallelmultilevel algorithms for hypergraph partitioning［J］.Journal of Parallel and Distributed Computing，2008，68(4):563-581.

［11］ Kernighan BW，LIN S.An efficient heuristic procedure for partitioning graphs［J］.The Bell System Technical Journal，1970，49(2):291-307.

［12］ Akhras G，DhattG.An automatic node relabeling scheme forminimizing amatrix or network bandwidth［J］.International Journal for Numerical Methods in Engineering，1976，10(4):787-797.

［13］趙玉靜，吳建軍，符文貞.板料成形數(shù)值模擬中網(wǎng)格節(jié)點編號優(yōu)化［J］.鍛壓裝備與制造技術(shù)，2009，44(6):73-75.ZHAO Yujing，WU Jianjun，F(xiàn)U Wenzhen.An algorithm for optimizing the node number of the numerical simulation in sheet forming［J］.China Metalforming Equipment and Manufacturing Technology，2009，44(6):73-75.(in Chinese)

［14］胡志宇，鄒琦，陳濤.改進(jìn)的有限元網(wǎng)格編碼AD算法［J］.三峽大學(xué)學(xué)報:自然科學(xué)版，2009，31(2):63-65.HU Zhiyu，ZOUQi，CHEN Tao.An improved AD algorithm for finite elementsmesh code optimization［J］.Journalof China Three Gorges University:Natural Sciences，2009，31(2):63-65.(in Chinese)