王 姝，王小鴿，楊廣文

清華大學 計算機系，北京100084

1 引言

隨著燃燒數(shù)值模擬計算技術的發(fā)展，模擬的問題越來越復雜，模擬的精度要求越來越高，傳統(tǒng)模擬方法采用結構網(wǎng)格，計算較為簡單，但是受到其結構特性的限制，不能靈活地根據(jù)模擬情況改變網(wǎng)格的結構以適應新的計算需求，因此逐漸被非結構網(wǎng)格所替代。相對于結構網(wǎng)格而言，非結構網(wǎng)格可以對復雜的幾何外形結構進行更為準確的描述；非結構網(wǎng)格中的一個點周圍的點數(shù)和單元數(shù)不是固定的，這一特點使得非結構網(wǎng)格易于調(diào)整網(wǎng)格結構，以適應計算的需要。

在進行燃燒數(shù)值模擬計算的過程中，部分區(qū)域的計算精度要求比其他區(qū)域高，如化學反應劇烈的區(qū)域、流場梯度變化劇烈的區(qū)域。為了提高計算精度，需要對此類區(qū)域的網(wǎng)格進行自適應加密，通過減小計算網(wǎng)格的尺度以達到需要的精度。另外，在某些不需要高計算精度的區(qū)域，為了減少計算量，對網(wǎng)格進行自適應減疏。

由于數(shù)值計算規(guī)模越來越大，傳統(tǒng)的串行程序已經(jīng)不能滿足高效的計算需求，將數(shù)值模擬計算并行化成為迫切需求。同時，在計算并行化以后，由于計算集中在化學反應劇烈的區(qū)域，進程之間的計算負載會嚴重不平衡，將導致并行效率大大降低。非結構網(wǎng)格和網(wǎng)格自適應結合使用可以達到模擬計算的精度要求，但是含自適應功能以后，網(wǎng)格數(shù)目的變動會使負載隨之變化。在這種情況下，需要在進程之間進行負載平衡調(diào)度，使計算負載均勻分布在所有進程上，使得非結構網(wǎng)格自適應計算模擬在提高精度的同時，縮短程序運行的時間，提高數(shù)值模擬計算的效率。

在自適應非結構網(wǎng)格上實現(xiàn)程序并行及負載平衡，最大的挑戰(zhàn)是如何在原有數(shù)據(jù)結構的基礎上正確、高效地執(zhí)行程序。本文是基于文獻[1]中的燃燒數(shù)值模擬和文獻[2]中的非結構網(wǎng)格自適應方式，在充分研究了非結構網(wǎng)格的數(shù)據(jù)結構和非結構網(wǎng)格自適應變化特點的基礎上，提出了一套有效的負載平衡策略，使得進程之間的負載較為平均，減少了程序總體運行時間；為了使負載轉移不受非結構網(wǎng)格存儲方式的限制，本文對非結構網(wǎng)格的存儲方式進行了改進，使得負載轉移更加容易實現(xiàn)。

2 相關工作

在非結構網(wǎng)格上實現(xiàn)并行和負載平衡計算逐漸成為數(shù)值模擬計算的一個重要研究方向，近幾十年有許多這方面的系統(tǒng)和算法。國際上的許多靜態(tài)負載平衡方法[3-6]是在編譯階段預測負載，并將數(shù)據(jù)和計算映射到進程上，這對于常規(guī)并行問題是很有效的，但是當問題變得非常規(guī)時（例如需要引用非本地數(shù)據(jù)或者計算迭代的開銷在運行時無法預測）性能會變得很差。SUPPLE[7]結合了靜態(tài)和動態(tài)負載平衡的方法，首先通過靜態(tài)的方式在編譯階段預測負載和通信，程序運行時采用塊為單位進行負載轉移。還有很多對負載的劃分基于圖劃分的方法[8-11]，它們的網(wǎng)格劃分方法受到非結構網(wǎng)格和自適應機制的約束，很難應用到現(xiàn)有的數(shù)值模擬計算中。

在國內(nèi)的研究中，DSMC[12]并行算法的負載平衡的指標是每個進程上的網(wǎng)格上的分子數(shù)量大致相當，采用交替二分法劃分網(wǎng)格，這種劃分方法的負載平衡效率較好，但會產(chǎn)生許多不規(guī)則的通信，容易造成程序效率低下。文獻[13]采用預處理方法采用METIS 進行負載平衡和任務劃分，并提出一種網(wǎng)格單元重排序算法提高訪存命中率難點以及研究的目標。這種動態(tài)負載平衡效果較好，但是需要改變原有數(shù)據(jù)結構的存儲方式。還有文獻[14]采用子區(qū)域數(shù)遠遠小于處理器數(shù)目的方式實現(xiàn)動態(tài)負載，這樣的處理方式可以得到較為均衡的負載，但是由于問題劃分得過小，需要經(jīng)常通信，會降低程序執(zhí)行的效率。文獻[15]在全局范圍內(nèi)實現(xiàn)網(wǎng)格的重新劃分，這種處理方式會帶來大量、經(jīng)常的通信。

國外的負載平衡系統(tǒng)和算法絕大部分是不開源的，無法應用到現(xiàn)有的數(shù)值模擬計算中。而國內(nèi)的處理方法大都需要靈活網(wǎng)格的存儲方式，例如文獻[13]的網(wǎng)格重排序和劃分、文獻[15]的網(wǎng)格重新劃分等。要在盡量不改變非結構網(wǎng)格自適應存儲方式的前提下達到負載平衡的目的，需要一種對數(shù)據(jù)存儲方式要求不高的負載平衡方法。

在自適應非結構網(wǎng)格上實現(xiàn)動態(tài)負載平衡需要考慮其數(shù)據(jù)結構和自適應方式。每個網(wǎng)格單元上的數(shù)值需要由其相鄰網(wǎng)格單元上的數(shù)值計算得出，如圖1 的偽代碼所示。

圖1 單元格計算的偽代碼

非結構網(wǎng)格中，相鄰網(wǎng)格的存儲地址是通過兩個網(wǎng)格共享邊的單元格屬性數(shù)組獲得的。如圖2 所示。

圖2 邊的編號、方向和單元格的關系

從圖2 中可以看出，編號為I的單元格的四條邊從右邊開始逆時針編號，依次為NEG1(I) ，NEG2(I)，NEG3(I)，NEG4(I)。邊被賦予了方向，其中水平邊以向右為正方向，豎直邊以向上為正方向。沿著邊的正方向，左邊的單元格為該邊的左單元格，右邊的單元格為該邊的右單元格。編號為J的邊的左右單元格分別為NL0(J)和NR0(J)。要想獲得單元格I的鄰居，以右鄰居為例，則應先取它的右邊NEG1(I)，然后取這條右邊的右屬性單元NR0(NEG1(I))，從而得到它的右鄰居。

網(wǎng)格的自適應加密方式如圖3 所示，1 號網(wǎng)格加密變成4 個子單元格，子單元格還可以再加密變成更小的4 個子單元格，以3 號網(wǎng)格為例。

圖3 網(wǎng)格自適應加密方式

偏微分方程求解和網(wǎng)格自適應需要用到鄰居單元格的數(shù)據(jù)，而鄰居單元格通過圖2 所示的關系得到。并行求解時，兩個相鄰網(wǎng)格被分配到不同的進程上，它們共享的邊（即重疊邊）需要分別在這兩個進程上存放。在自適應加密或減疏以后，由于兩個重疊邊的分裂情況不同，需要更新分裂后子邊對應的新重疊邊。如圖4（a）所示，初始時有兩條均未分裂的重疊邊I，J，邊J先分裂成兩個子邊，而I未分裂（如圖4（b）所示），這時邊I的重疊邊仍為J，邊J的子邊SUB1(J)和SUB2(J)的重疊邊為邊I。當邊I也發(fā)生分裂時，邊J子邊的重疊邊才發(fā)生更新。如圖4（c）所示，I的子邊SUB1(I)的重疊邊為SUB1(J)，SUB2(I)的重疊邊為SUB2(J)。

本文使用的實驗是一個化學反應流計算程序。程序執(zhí)行1 280 步，采用8 個進程（便于在圖上顯示）。初始網(wǎng)格為16×4，化學反應穩(wěn)定后的網(wǎng)格規(guī)模在20 000 左右?；瘜W反應使用常微分方程求解，流動使用偏微分方程求解。通過實驗得出，各進程解常微分方程的時間差異很大，而解偏微分方程的時間差異不大，因此計算不平衡部分集中在解常微分方程部分。本文關注的重點在于求解常微分方程計算時間的差異，而非通常意義下的負載，即在求解過程中做了多少、什么種類的運算。為說明簡便，本文后續(xù)內(nèi)容中的負載均指求解常微分方程的計算負載。考慮進程間的負載變動情況，統(tǒng)計程序執(zhí)行計算的時間（s），為使結果更清晰，將時間t取對數(shù)ln t，各進程的執(zhí)行時間情況如圖5所示。

圖4 原網(wǎng)格相鄰網(wǎng)格被分配到不同的進程示意圖

圖5 不同時間步上各進程的負載分布

從圖5 可以看出，在程序執(zhí)行初期，從1 號進程開始，負載逐漸變大，并且隨著時間的推進各進程的負載陸續(xù)增加，這是第一階段；直到第400 步左右，7，8 號進程的負載劇烈增加，這是第二階段；到第600 步以后，各進程的負載趨于穩(wěn)定，這是第三階段。出現(xiàn)這種情況的原因是：第一階段化學反應不明顯，流場梯度變化劇烈的區(qū)域網(wǎng)格發(fā)生加密，隨著該區(qū)域的位置不斷轉移，各進程的負載隨之變化。在第二階段，在7，8 進程的計算區(qū)域突然發(fā)生了劇烈的化學反應，導致這些區(qū)域的負載急劇增加。在第三階段，化學反應趨于穩(wěn)定，負載變動緩慢。

在文獻[1]的工作基礎上，為了提高計算效率，將原有并行程序進行優(yōu)化。本文在原有數(shù)據(jù)結構上，根據(jù)非結構網(wǎng)格和網(wǎng)格自適應方式的特點，對非結構網(wǎng)格的存儲方式進行了優(yōu)化，提出了一套高效的負載轉移方法。

3 需要考慮的問題

要在自適應非結構網(wǎng)格上實現(xiàn)動態(tài)負載平衡，需要考慮以下六個方面：

（1）非結構網(wǎng)格

在傳統(tǒng)的結構網(wǎng)格上，網(wǎng)格單元使用多維數(shù)組存放，并且可以根據(jù)該數(shù)組的下標關系來訪問鄰接單元格，每個單元格的鄰接單元格是確定的，這便于網(wǎng)格的剖分；而非結構網(wǎng)格通過一維數(shù)組存放單元格，鄰居單元格的數(shù)據(jù)調(diào)用方式如圖2 所示。

（2）自適應

當一個網(wǎng)格單元與其相鄰網(wǎng)格的梯度超過閾值時，網(wǎng)格在型心點處分裂成四個，稱為加密。原來的網(wǎng)格稱為父網(wǎng)格，四個新產(chǎn)生的網(wǎng)格稱為子網(wǎng)格。當一個父網(wǎng)格單元與其相鄰網(wǎng)格的梯度低于某一閾值時，其四個子網(wǎng)格合并為一個，稱為減疏。網(wǎng)格的自適應會導致某些網(wǎng)格上的計算急劇增加，加劇了進程之間的負載不平衡性，更加大了負載平衡的難度。非結構網(wǎng)格的這種特點使得動態(tài)負載過程中的負載劃分變得困難。

（3）負載的變動

數(shù)值模擬計算一般包括兩個方面：求解常微分方程和求解偏微分方程，本文使用的程序中，化學反應使用常微分方程求解，流動使用偏微分方程求解。通過分析程序可以得出：不同進程之間的化學反應計算負載差異巨大；流動計算僅與網(wǎng)格數(shù)目有關，而化學反應穩(wěn)定以后，各個進程中的網(wǎng)格數(shù)量基本一致。常微分方程通過牛頓迭代求解，在劇烈化學反應開始之前，常微分方程求解計算負載僅與網(wǎng)格數(shù)有關；在劇烈化學反應開始之后，牛頓迭代次數(shù)急劇上升，同一單元格上的計算負載相差幾百倍，因此可以得出結論：同一個網(wǎng)格的計算負載與當前牛頓迭代次數(shù)相關。

由于負載的變動，需要解決的問題包括負載預測的方法和負載轉移的時機。即負載的變動與哪些因素有關，如何將這些因素與負載聯(lián)系起來，以及在執(zhí)行哪些計算時，需要負載轉移。

（4）負載轉移的粒度

根據(jù)數(shù)值計算系統(tǒng)的計算特點，求解偏微分方程需要由其鄰居單元格上的數(shù)值計算得出，采用原始粗網(wǎng)格中的一列作為負載轉移的單位，可以降低負載轉移后尋找鄰居單元格的難度。

（5）非結構網(wǎng)格自適應存儲方式帶來的困難

圖6 原始網(wǎng)格

通過結合圖6 和圖7 來分析非結構網(wǎng)格自適應時存儲數(shù)組的特點。假設原始網(wǎng)格中有1～9 號非父（沒有子單元格）的單元格，按列先序一起存儲在圖7（a）所示的數(shù)組中。非父單元格占據(jù)數(shù)組的前半部分位置，父單元格占據(jù)數(shù)組的后半部分位置，中間未填滿的部分留空。原始未加密網(wǎng)格的存儲方式，如圖7（a）所示，此時通過計算，2 號單元格需要加密。那么2 號單元格的子單元格編號為10，11，12，13，放入非父單元格末尾，如圖7（b）所示。此時2 號單元格已經(jīng)加密成為父單元格，需要將它挪入父單元格存儲區(qū)域，即從數(shù)組末尾開始往前存放。然后將末尾的13 號單元格填補上2號單元格在非父單元格存儲區(qū)域中留下的空白，如圖7（c）所示。

圖7 網(wǎng)格自適應存儲方式

由于負載轉移的單位為以原始粗網(wǎng)格中的一列，而非結構網(wǎng)格的存儲方式是不以列為單位的，這樣一來，要從如圖7 所示的非結構網(wǎng)格存放的一維數(shù)組中挑選出第一列（1，3，10，11，12，13 號單元格）上的數(shù)據(jù)執(zhí)行轉移，不便于處理。這里需要解決的問題是如何改進網(wǎng)格自適應存儲方式，達到第4 章中提出的以列為單位執(zhí)行負載轉移。

（6）負載轉移策略

負載轉移策略的目標是使通信量盡可能少，每個進程上的負載趨于一致。因此在每次負載轉移之前，需要一個進程統(tǒng)計所有進程的負載情況，各進程負載的最大值與平均負載的比值為Max/Avg 。在負載平衡的情況下，Max/Avg為1。不平衡率越大，負載越不平衡。在進行負載調(diào)度時，判斷全局的Max/Avg 是否大于某個給定的閾值，如果大于，則執(zhí)行負載轉移。這個閾值的大小是根據(jù)負載轉移的開銷和不平衡情況綜合考慮的。遷移的開銷越大，閾值的設置也需要大些，以避免負載平衡工作的“得不償失”。

4 解決方案

本文的實驗使用的計算節(jié)點采用兩個Intel Xeon X5670六核處理器（2.93 GHz，12 MB Cache），160 GB SATA 硬盤。節(jié)點內(nèi)存48 GB，結果取多次實驗的平均值。

結合上一章中提出非結構網(wǎng)格的特點和網(wǎng)格自適應方式，針對需要解決的問題，本文提出的負載平衡方案具體如下。

4.1 負載預測方法

經(jīng)過分析串行程序各部分運行時間，解常微分方程時間占解常微分方程和解偏微分方程總時間的34.4%。偏微分方程計算時間只與網(wǎng)格數(shù)目有關，常微分方程計算約90%的時間花費在牛頓迭代上，且臨近時間步的牛頓迭代次數(shù)相近，因此將本時間步牛頓迭代次數(shù)作為下一時間步負載的估計。

為了更加精確地預測負載，需要選取能夠進一步表征常微分方程計算負載的變量進行預測。在本實驗情況下，選取化學反應放熱量Q。在放熱量達到某個閾值時，開始劇烈化學反應，將此閾值作為切換負載預測策略的分割點。采用這種自動切換負載預測方法的原因是：開始劇烈化學反應之前，放熱量Q 不穩(wěn)定，無法作為負載預測的判據(jù)，此時解常微分方程的負載與牛頓迭代次數(shù)成正比；在劇烈化學反應開始之后，放熱量Q 與負載的相關性更高。通過結合兩種預測方式，可以使負載的預測更接近真實情況。這種負載預測方法的效果在第5 章中分析。

4.2 負載轉移時機

考慮到第3 章中提出的負載變動方式，由于解偏微分方程的計算負載只與網(wǎng)格數(shù)有關，化學反應穩(wěn)定后各進程上的網(wǎng)格數(shù)目差異不大，求解偏微分方程計算負載基本平衡，因此只需在解常微分方程時執(zhí)行動態(tài)負載平衡策略。

由主控進程統(tǒng)計所有計算進程的負載，若Max/Avg 大于某個閾值，則需要進行負載轉移。

另外，為了避免自適應時重疊邊的更新成為難點，采取在自適應之后執(zhí)行負載轉移這種方式。這樣僅需把計算負載轉移出去，計算完成之后更新單元格的狀態(tài)即可。

4.3 數(shù)據(jù)存儲方式的調(diào)整

由第3 章的描述可以得知，非結構網(wǎng)格的存儲結構特點是：單元格在存儲時沒有規(guī)律可循，要從所有網(wǎng)格中挑出某一列單元轉移到其他進程上，需要使用數(shù)組標記每個單元格所在的列，當執(zhí)行負載轉移時，遍歷標記數(shù)組，從中挑出所需的單元格，并放入一塊連續(xù)的數(shù)組中。為了避免這種方式帶來的麻煩，改為修改非結構網(wǎng)格的存儲方式，將整個網(wǎng)格存儲地址以列為單位進行劃分，如圖8 所示，各列分別占據(jù)一塊連續(xù)的空間，與圖7 相同，2 號單元格加密，其所在列的任何單元格的增加和變動都在這一塊空間中，要轉移第一列（1，3，10，11，12，13 號單元格）不需挑選，直接轉移即可。這使得以列為單位的負載轉移成為可能，并且程序在執(zhí)行時有較高的空間局部性。

圖8 以列為單位的存儲方式

4.4 負載轉移策略

本文采用的轉移策略為進程一對一進行負載轉移，這種遷移是一種簡單易行的方法，在其基礎上可以實現(xiàn)更為復雜的遷移策略。

負載平衡分為兩個階段：

負載轉移過程：將進程按照預測的負載大小排序，負載轉移的規(guī)則為最大負載的進程轉移給最小負載的進程，次大負載的進程轉移給次小負載的進程，以此類推。在過載進程（即進程負載大于平均負載）的內(nèi)部，負載以列為單位進行排序，要轉移的負載為從最大負載列開始，隔列轉移。這樣做的原因是由于燃燒的穩(wěn)定性使得相鄰網(wǎng)格列上的負載接近，通過排序并隔行轉移可以使得轉移的源進程與目的進程的負載相近，避免過多或過少轉移負載，引入新的負載不平衡。主控進程標記轉移的方向及數(shù)量，并廣播給所有計算進程。然后根據(jù)上一步收到的標記數(shù)組，負載小于平均負載的進程（欠載進程）準備接收數(shù)組，并擴大計算涉及的輸入及輸出數(shù)組，從過載進程接收計算需要的輸入數(shù)組包，拆包將數(shù)據(jù)放入適當?shù)奈恢茫贿^載進程準備發(fā)送數(shù)組，將計算需要的輸入數(shù)組打包，并發(fā)送給對應的欠載進程。

在計算完成之后，原來的欠載進程將計算后改變的輸出數(shù)組打包，發(fā)送給原來的過載進程。并恢復計算涉及的輸入及輸出數(shù)組為原來的規(guī)模。原來的過載進程準備接受數(shù)組包，并將其放入原先的位置。

5 效果的評價

本文從每一個計算步中，進程的負載統(tǒng)計、有負載平衡策略與無負載平衡策略時程序的運行時間、進程之間的負載不平衡程度等方面對負載平衡的策略進行分析。

（1）負載不平衡程度的改善

由于所有進程中負載最大的進程執(zhí)行時間最長，直接影響到程序執(zhí)行的總時間；而平均負載是理想情況下每個進程的執(zhí)行時間。

有負載轉移和無負載轉移時，最大負載與平均負載的比值如圖9 所示，橫軸為程序執(zhí)行的步驟，灰色實線為無負載轉移時最大負載與平均負載的比值，虛線部分為使用牛頓迭代次數(shù)預測負載策略時最大負載與平均負載的比值，黑色實線部分為結合兩種預測方式后最大負載與平均負載的比值。從圖中可以看出，使用牛頓迭代次數(shù)預測在絕大部分時間可以較好地降低最大負載與平均負載的比值，但是在400 步左右由于化學反應劇烈，這種策略就不夠有效。加入放熱量Q 預測負載后，在化學反應劇烈的階段也能很好地預測負載變動趨勢，最大負載與平均負載的比值降低到無負載轉移的一半左右。

圖9 最大負載與平均負載的比值

有負載轉移和無負載轉移時，各進程的負載如圖10所示。

圖10 不同燃燒階段進程上的負載對比

圖10的橫軸為各進程編號，縱軸為該階段的累計負載，淺色直方圖為無負載轉移的統(tǒng)計負載，深色直方圖為有負載轉移的累計負載。圖10（a）顯示了劇烈化學反應開始之前，在各進程上無負載轉移和有負載轉移的負載，有負載轉移的部分負載更為平均。圖10（b）顯示了劇烈化學反應開始之后，無負載轉移和有負載轉移的負載累計，有負載轉移的部分仍然更為平均。從圖10 還可以看出，有負載轉移的最大負載接近無負載轉移的最大負載的一半，進程間的負載在有負載轉移時更加平衡。

（2）負載轉移時間分析

為了使負載轉移有效地執(zhí)行，即負載轉移的開銷不能大于負載轉移帶來的性能提升，因此需要度量負載轉移帶來的開銷。

由于每一個迭代步中，影響程序執(zhí)行時間的部分為負載最大的進程的執(zhí)行之間，因此使用每個迭代步負載最大的進程中，負載轉移開銷占總的求解常微分方程計算時間的比值作為負載轉移的開銷。如圖11 所示。

圖11 負載轉移所占時間比例

由圖11 中可以看出，負載轉移所占時間比例最大不超過總執(zhí)行時間的35%，而根據(jù)負載轉移策略，是將過載進程與欠載進程負載差值的一半轉移出去，即計算時間的50%，即節(jié)省了50%的時間，而最大的開銷僅為30%，特別是劇烈化學反應負載比較大的情況下，負載轉移策略的開銷不到計算時間的5%，由此可以得出負載轉移一定是有效的。

（3）程序總體性能的提升

為了分析負載轉移策略對主程序帶來的性能提升，需要分析求解常微分方程計算部分占整個程序運行時間的百分比。經(jīng)過實驗得出，在程序的核心計算過程中，解常微分方程計算時間占總執(zhí)行時間的5%～35%。

定義性能提升的度量為在并行情況下，不同核數(shù)執(zhí)行程序核心計算過程時，和有負載平衡的執(zhí)行時間與無負載平衡的執(zhí)行時間減少時間的比值。測試網(wǎng)格初始時為24×4，化學反應穩(wěn)定時的規(guī)模在30 000 左右。

從表1中可以看出，程序的核心計算性能總體隨著核數(shù)的增加而提高，12 核運行時程序性能提升達到了10.20%。

表1 程序的性能提升

6 總結和未來工作

本文對非結構網(wǎng)格自適應的燃燒并行數(shù)值模擬計算的特點進行了詳細的研究，提出了數(shù)據(jù)按列存儲的改進方式和根據(jù)不同化學反應情況預測負載的機制，還提出了將負載一對一地轉移到欠載進程上動態(tài)負載平衡策略，將最大負載與平均負載的比值降低了一半左右，程序總體性能提升在9.00%～10.20%。由于一對一負載轉移的限制，還有部分進程處于空閑狀態(tài)，下一步的工作是將負載轉移改為一對多，將過載進程上的負載盡可能多地分出去，使得所有進程上的負載更加平衡。

[1] Chen Z，Burke M P，Ju Y.Effects of compression and stretch on the determination of laminar flame speeds using propagating spherical flames[J].Combustion Theory and Modelling，2009，13（2）：343-364.

[2] Sun，Takayama.Conservative smoothing on an adaptive quadrilateral grid[J].Journal of Computational Physics，1999，150（1）：143-180.

[3] Hummel S F，Schonberg E，F(xiàn)lynn L E.Factoring：a method for scheduling parallel loops[J].Comm of the ACM，1992，35（8）：90-101.

[4] Kumar V，Grama A Y，RaoVempaty N.Scalable load balancing techniques for parallel computers[J].J of Parallel and Distr，1994，22（1）：60-79.

[5] Lain M S，Rothberg E E，Wolf M E.The cache performance and optimizations of block algorithms[C]//4th Int Conf on Architectural Support for Progr Lang and Operating Systems，1991：63-74.

[6] Orlando S，Perego R.A template for non-uniform parallel loops based on dynamic scheduling and prefetching techniques[C]//Proc of the 1996 ACM Int Conf on Supercomputing，1996：117-124.

[7] Orlando S，Perego H.SUPPLE：an efficient run-time support for non-uniform parallel loops[J].J of System Architecture，1999，45（15）：1323-1343.

[8] Oliker L，Biswas R.PLUM：Parallel Load balancing for adaptive Unstructured Meshes[J].Parallel and Distributed Computing，1998，52（2）：150-177.

[9] Bui T.Graph bisection algorithms with good average behavior[J].Combinatorica，1984，7（2）：181-191.

[10] Kernighan B，Lin S.An efficient heuristic procedure for partitioning graphs[J].Bell System Technical Journal，1970，49（2）：291-307.

[11] Chan T，Mathew T.Domain decomposition algorithms[J].Acta Numerica，1994，3（1）：61-143.

[12] 王學德，伍貽兆，夏健.動態(tài)負載平衡的二維非結構網(wǎng)格DSMC并行算法研究[J].空氣動力學學報，2007，25（3）：339-361.

[13] 劉鑫，陸林生，陳德訓.非結構網(wǎng)格并行計算預處理方法研究[J].計算機科學，2012，39（3）：308-311.

[14] 陳建軍.非結構網(wǎng)格的并行劃分和自適應非結構化網(wǎng)格生成及其并行化的若干問題研究[D].杭州：浙江大學，2006.

[15] 李桂波，楊國偉.基于多塊結構網(wǎng)格的并行計算及負載平衡研究[J].宇航學報，2011，32（6）：1224-1230.