基于Akhras-Dhatt算法的區域剖分法

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

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

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

由于能夠對復雜外形進行準確、有效地擬合，無結構網格在有限單元法、有限差分法、有限體積法等數值方法中都有廣泛應用。目前，基于無結構網格進行并行計算時，常用的區域剖分方法有:陣面推進法(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］等。

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

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

基于AD算法思想，文中提出一種區域剖分方法。首先利用AD算法給無結構網格優化編號以減少鄰接矩陣的帶寬，然后根據優化后的單元編號順序和單元屬性，將計算單元均勻地映射到各個計算節點上。通過算例與其它剖分方法比較，結果表明該方法能夠實現各處理器上負載的更為優化的平衡，對處理大型數據具有可行性和有效性。

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

在通常的流體、固體等數值模擬計算中，當前單元物理量的計算需要調用其周邊單元的數據。當將網格進行區域剖分后，會造成部分需要數據傳輸的單元不在同一計算節點上，這就需要進行節點間的外部通信。適當減少外部通信量，可以提高并行計算效率。

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

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

引進

一般數值計算時會調用當前單元上的數據，因此取aii=1。若記，則該值與單元的具體編號無關，是鄰接矩陣的一個不變特征量。

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

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

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

其中

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

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

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

的計算公式，分兩種情況進行研究。

1)當K≤L+1≤N時

2)當1≤L+1≤K時

綜合以上兩種情況可以進一步得到結論:

2)從鄰接矩陣帶狀結構可見，第1個和最后一個計算節點需要進行外部通信的單元總數為，其它計算節點需要進行外部通信的單元總數為

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

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

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

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

2.1 AD算法簡介

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

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

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

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

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

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

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

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

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

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

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

2.3 基于AD算法的區域剖分法

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

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

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

3 算例

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

對于陸地邊界單元需要利用鏡像法構造虛擬單元并計算虛擬單元的物理量的值，故其計算時間會遠大于內單元的計算時間，區域分裂時必須區別對待內單元和邊界單元。

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

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

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

圖2 經典計算單元模板Fig.2 Classic cellcom puting template

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

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

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

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

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

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

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

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

4 結語

在AD算法優化編號的基礎上對網格單元利用序號屬性法進行區域剖分，能得到較高的加速比與并行效率。通過對高達59萬個無結構網格單元進行區域剖分試驗，結果表明利用AD分裂法進行區域剖分，確實能使各計算節點的負載較快速趨于平衡，得到高加速比與并行效率的區域剖分結果，因此該方法是一種可用于大數據量的網格單元進行區域剖分的可行且有效的方法。

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

［2］司海青，王同光，成娟.非結構網格上Eluer方程的區域分裂算法及其并行計算［J］.空氣動力學學報，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］周春華.并行計算中一種非結構網格分割方法［J］.航空學報，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］趙玉靜，吳建軍，符文貞.板料成形數值模擬中網格節點編號優化［J］.鍛壓裝備與制造技術，2009，44(6):73-75.ZHAO Yujing，WU Jianjun，FU 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］胡志宇，鄒琦，陳濤.改進的有限元網格編碼AD算法［J］.三峽大學學報:自然科學版，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)