專利名稱：使用并行多級(jí)不完全因式分解求解儲(chǔ)層模擬矩陣方程的方法
技術(shù)領(lǐng)域：
以下說明書一般涉及求解線性方程組的迭代求解器，更具體地涉及在高性能的并 行計(jì)算系統(tǒng)上在求解線性方程組的并行迭代過程中執(zhí)行預(yù)處理程序的系統(tǒng)和方法。
背景技術(shù)：
在分析許多科學(xué)或工程應(yīng)用時(shí)，通常需要同時(shí)求解大量的線性代數(shù)方程，這些方 程可以以矩陣方程的形式表示如下Ax = b，(以下稱為方程式⑴)，其中A表示已知的n*n 維的正方形系數(shù)矩陣，b表示通常稱為“右手邊”的已知η維向量，χ表示通過求解該方程組 得到的未知η維向量。已知各種技術(shù)用于求解這樣的線性方程組。線性方程組通常在各種 基于計(jì)算機(jī)的三維(3D)模擬或給定的真實(shí)世界系統(tǒng)的建模中遇到(并需要求解)。例如， 地下含有碳?xì)浠衔锏膬?chǔ)層(油或氣儲(chǔ)層)的現(xiàn)代三維(3D)模擬需要方程式(1)類型的 線性代數(shù)系統(tǒng)的求解，一般地，在稀疏系數(shù)矩陣A中有數(shù)百萬未知量以及數(shù)千萬和數(shù)億的 非零元素。這些非零元素定義稀疏矩陣結(jié)構(gòu)。同樣地，基于計(jì)算機(jī)的三維(3D)建?？梢杂脕砟M真實(shí)世界系統(tǒng)，比如機(jī)械或電 氣系統(tǒng)(如可應(yīng)用于汽車、飛機(jī)、輪船、潛艇、宇宙飛船等中)、人體(如整個(gè)人體或者人體的 各部分，如重要器官等的建模)、天氣形勢(shì)以及各種其他要被建模的真實(shí)世界系統(tǒng)。通過這 樣的建模，可以分析或預(yù)測(cè)被建模系統(tǒng)的未來潛在性能。例如，提供給被建模系統(tǒng)的某些變 化的條件對(duì)系統(tǒng)未來性能的影響可通過與基于計(jì)算機(jī)的模型的交互和分析進(jìn)行評(píng)估。作為一個(gè)例子來說，多孔介質(zhì)中流體流動(dòng)的建模是石油工業(yè)的主要焦點(diǎn)。不同的 基于計(jì)算機(jī)的模型用于石油工業(yè)的不同領(lǐng)域，但大多數(shù)包括用偏微分方程系統(tǒng)或方程組 (PDE' s)來描述該模型。一般地，這樣的建模通常要求在給定網(wǎng)格上在空間和時(shí)間上對(duì)偏 微分方程系統(tǒng)進(jìn)行離散化，并為每個(gè)時(shí)間步長或時(shí)步執(zhí)行計(jì)算直到達(dá)到規(guī)定的時(shí)間。在每 個(gè)時(shí)步，離散方程被求解。通常離散方程是非線性的，并且求解過程是迭代的。非線性迭代 方法的每個(gè)步驟通常包括非線性方程系統(tǒng)或方程組的線性化(如，雅可比構(gòu)造)，求解該線 性方程組以及用于計(jì)算下一個(gè)雅可比構(gòu)造的性質(zhì)計(jì)算。

圖1示出一般的工作流程，其通常用于地下含碳?xì)浠衔飪?chǔ)層中隨時(shí)間的流體流 動(dòng)的基于計(jì)算機(jī)的模擬(或建模)。內(nèi)部循環(huán)101是求解非線性系統(tǒng)的迭代方法。再者，經(jīng) 過內(nèi)部循環(huán)101的每次操作一般都包括非線性方程組的線性化(如，雅可比構(gòu)造)11，求解 線性系統(tǒng)(或方程組)12以及用于計(jì)算下一個(gè)雅可比(當(dāng)循環(huán)返回到框11時(shí))的性質(zhì)計(jì) 算13。外部循環(huán)102是時(shí)步循環(huán)(time step loop)。如所示的，針對(duì)每一個(gè)時(shí)步循環(huán)的邊界條件可在框10中定義，然后在執(zhí)行針對(duì)時(shí)步的內(nèi)部循環(huán)101后，為時(shí)步計(jì)算的結(jié)果可在 框14中輸出(如結(jié)果可存儲(chǔ)到數(shù)據(jù)存儲(chǔ)介質(zhì)或/和提供給軟件應(yīng)用程序以便生成顯示，該 顯示表示為相應(yīng)時(shí)步建模的地下含碳?xì)浠衔飪?chǔ)層中的流體流動(dòng))。如上述所示，除了地下 含碳?xì)浠衔飪?chǔ)層中流體流動(dòng)的建模之外，真實(shí)世界系統(tǒng)的基于計(jì)算機(jī)的3D建?？梢砸?類似方式執(zhí)行，即可采用求解線性方程組的迭代方法(如圖1的框12中所示)。線性方程組的求解是一種計(jì)算強(qiáng)度非常高的任務(wù)，因此需要高效的算法。有兩類 一般的線性求解器1)直接法和2)迭代法。所謂的“直接法”基于高斯消元法(Gaussian elimination)，該方法中對(duì)矩陣A進(jìn)行因式分解，并表示為下三角形矩陣L和上三角形矩陣 U(因式)的積A = LU (以下稱為方程式O))。但對(duì)大稀疏矩陣A而言，計(jì)算三角形矩陣 L和U非常耗費(fèi)時(shí)間，且那些因子中的非零元素的數(shù)量會(huì)很大，因此它們可能不適合甚至現(xiàn) 代高性能計(jì)算機(jī)的存儲(chǔ)器。迭代法是基于簡單而通常低廉操作的反復(fù)應(yīng)用，如矩陣向量積，其提供具有給定 精確度的近似解。通常，對(duì)于出現(xiàn)在科學(xué)或工程應(yīng)用中的方程式(1)類型的線性代數(shù)問題 而言，系數(shù)矩陣的性質(zhì)導(dǎo)致收斂到一個(gè)解的大量迭代。為了減少迭代次數(shù)，并因此減少運(yùn)用迭代方法求解矩陣方程的計(jì)算開銷，經(jīng)常使 用預(yù)處理技術(shù)，在該技術(shù)中方程式(1)類型的初始矩陣方程與合適的預(yù)處理矩陣M(也可 簡稱為預(yù)處理器(preconditioner))相乘，如Hx = M^b (以下稱為方程式⑶)。這 里，M—1表示矩陣M的逆矩陣。采用不同的預(yù)處理方法(矩陣M)能極大減少以足夠高的 準(zhǔn)確度計(jì)算方程式(1)的合適解的運(yùn)算開銷。預(yù)處理技術(shù)的主要例子是代數(shù)多網(wǎng)格方法 (algebraic multi-grid methods)禾口不完全上下三角因式分角軍(incomplete lower-upper factorizations)0在第一種方法(即多網(wǎng)格方法)中，構(gòu)造一系列減維的系數(shù)矩陣，且建立從較精細(xì) 維到較粗維的一些數(shù)據(jù)轉(zhuǎn)移方法。然后，近似求解矩陣方程(1)(所謂的“平滑”)，計(jì)算余 項(xiàng)^· = “+，并將得到的向量^·傳遞到較粗維(所謂的“約束”)。因此可在較粗維中近似 求解類似方程式(1)的方程式，計(jì)算余項(xiàng)并傳遞到較粗維等等。當(dāng)在最粗維上計(jì)算該問題 之后，此粗解會(huì)被傳遞回原始維(所謂的“延長”)，以得到將被添加在原始精細(xì)維上的近似 解的缺陷。另外一個(gè)預(yù)處理技術(shù)的例子是不完全上下三角因式分解(ILU型)，該方法取代了 完全因式分解(如方程式( 所示)，計(jì)算稀疏因子L、U使得其乘積近似等于初始系數(shù)矩 陣A LU (以下稱為方程式G))。上述的兩種預(yù)處理技術(shù)本質(zhì)上是有順序的，且不能直接應(yīng)用在并行處理計(jì)算機(jī) 上。隨著出現(xiàn)在科學(xué)和工程應(yīng)用中的代數(shù)問題的維數(shù)增加，對(duì)適合并行處理計(jì)算機(jī)的求解 方法的需求變得越來越重要。因此，研發(fā)高效的并行線性求解器成為日益重要的任務(wù)，這對(duì) 許多3D建模應(yīng)用諸如石油儲(chǔ)層建模而言尤其重要。盡管很多不同的用大稀疏系數(shù)矩陣求 解矩陣方程的方法有實(shí)質(zhì)進(jìn)步，例如多網(wǎng)格或直接求解器，但是在過去的數(shù)十年中，采用基 于不完全上下三角因式分解的預(yù)處理的迭代方法仍然是求解大型稀疏線性系統(tǒng)最普及的 方法。如上所述，這些預(yù)處理技術(shù)實(shí)質(zhì)上是有順序的，且不能直接應(yīng)用在并行處理計(jì)算機(jī) 上。近年來在科學(xué)界，一種新型的使用多級(jí)塊ILU因式分解技術(shù)的并行預(yù)處理策略
6被開發(fā)用于求解大型稀疏線性系統(tǒng)。這個(gè)新方法的大體思想是對(duì)未知量和對(duì)應(yīng)的方程式 重排序，并且將初始矩陣拆分成2*2的塊結(jié)構(gòu)，該拆分方式使得第一對(duì)角塊成為塊對(duì)角矩 陣。該塊矩陣會(huì)被并行進(jìn)行因式分解。通過消去因式分解后的矩陣塊形成舒爾補(bǔ)(Schur complement)之后，反復(fù)執(zhí)行該步驟以得到舒爾補(bǔ)。新方法的效率依賴于初始矩陣和舒爾 補(bǔ)被拆分成塊的方式。在傳統(tǒng)的方法中，多級(jí)因式分解是基于稀疏矩陣結(jié)構(gòu)原始圖的多著 色或塊獨(dú)立集分割。這些技術(shù)進(jìn)一步詳述在以下文獻(xiàn)中a)C.Shen，J. ^iang和K. Wang， Distributed block independent set algorithms and parallel multi-level ILU preconditioner. J. Parallel Distrib. Comput. 65(2005)，pp.331-346 ；禾口 b)Z. Li, Y. Saad 禾口 M. Sosonkina, ρARMS :A parallel version of the algebraic recuisive multilevel solver, Number. Linear Algebra App 1.，10(2003)，pp. 485-509，該文獻(xiàn)的公開內(nèi)容通過引 用包括進(jìn)本發(fā)明。這些方法的缺點(diǎn)是它們改變矩陣的初始排序，這在很多情況下會(huì)導(dǎo)致預(yù) 處理器的質(zhì)量更糟，或迭代求解器的收斂速度更慢。另一個(gè)缺點(diǎn)是這樣的并行重排序結(jié)構(gòu) 不能良好地?cái)U(kuò)展，即隨著處理單元(處理器)數(shù)量的增加，該結(jié)構(gòu)的質(zhì)量和效率明顯下降。另外一類基于ILU因式分解的平行預(yù)處理策略使用源于域分解(domain decomposition)的思想。給定大型稀疏線性方程式系統(tǒng)(1)，首先使用分割軟件(例如 G. Karypis 禾口 V. Kumar 在文獻(xiàn) METIS -Unstructured Graph Partitioning and Sparse Matrix Ordering System, version 4. 0，1998年九月中敘述的)，該矩陣A被分割成給定數(shù) 量的子矩陣P，幾乎每個(gè)子矩陣的行數(shù)都相同，且子矩陣間的連接很少。在分割步驟之后，進(jìn) 行局部重排序，首先要對(duì)每個(gè)子矩陣的內(nèi)部行進(jìn)行排序，然后是子矩陣的接口行，即那些與 其他子矩陣有連接的行。然后，已分割的以及序列改變的初始矩陣A可表示在如下的加邊 對(duì)角塊(BBD)形式中
權(quán)利要求
1.一種方法，其包括(a)使用多級(jí)分割將初始矩陣分割成多個(gè)子矩陣；(b)對(duì)對(duì)角塊并行執(zhí)行截?cái)嘁蚴椒纸?，為所述多個(gè)子矩陣的每一個(gè)的接口部分形成局 部舒爾補(bǔ)，；(c)由所述多個(gè)子矩陣的對(duì)角塊上的局部舒爾補(bǔ)和非對(duì)角塊上的所述多個(gè)子矩陣之間 的連接形成全局接口矩陣；(d)判定以下至少一個(gè)i)該全局接口矩陣是否足夠小，以滿足預(yù)定義的大小閾值，且 )是否達(dá)到最后的允許級(jí)；以及(e)當(dāng)判定該全局接口矩陣不是小得足以滿足所述預(yù)定義的大小閾值或者達(dá)到該最后 的允許級(jí)，則使用多級(jí)分割將所述全局接口矩陣分割成多個(gè)第二子矩陣和為所述多個(gè)第二 子矩陣重復(fù)操作(b)到(e)。
2.根據(jù)權(quán)利要求1所述的方法，進(jìn)一步包括(f)當(dāng)判定該全局接口矩陣足夠小以滿足所述預(yù)定義的大小閾值時(shí)，因式分解所述全 局接口矩陣。
3.根據(jù)權(quán)利要求1所述的方法，進(jìn)一步包括 對(duì)所述多個(gè)子矩陣的每一個(gè)進(jìn)行重排序。
4.根據(jù)權(quán)利要求3所述的方法，其中所述重排序在操作(a)之后且操作(b)之前執(zhí)行。
5.根據(jù)權(quán)利要求3所述的方法，其中所述重排序包括 首先重排序所述多個(gè)子矩陣的每個(gè)的內(nèi)部行；和 其次重排序所述多個(gè)子矩陣的接口行。
6.根據(jù)權(quán)利要求1所述的方法，進(jìn)一步包括在所述多個(gè)子矩陣上執(zhí)行局部縮放算法，以改進(jìn)所述子矩陣的數(shù)值性質(zhì)。
7.根據(jù)權(quán)利要求1所述的方法，其中所述形成所述全局接口矩陣包括 隱式形成所述全局接口矩陣。
8.根據(jù)權(quán)利要求7所述的方法，其中所述隱式形成所述全局接口矩陣包括 通過多個(gè)處理單元中的每一個(gè)，存儲(chǔ)所述接口矩陣的相應(yīng)部分。
9.根據(jù)權(quán)利要求1所述的方法，其中所述并行執(zhí)行包括 通過多個(gè)并行處理單元執(zhí)行所述操作(b)。
10.根據(jù)權(quán)利要求1所述的方法，其中所述使用多級(jí)分割將所述全局接口矩陣分割成 所述多個(gè)第二子矩陣包括使用所述接口矩陣的多級(jí)分割以避免在主處理單元顯式形成所述接口矩陣。
11.根據(jù)權(quán)利要求1所述的方法，進(jìn)一步包括 對(duì)最后一級(jí)接口矩陣的處理。
12.根據(jù)權(quán)利要求11所述的方法，其中所述對(duì)最后一級(jí)接口矩陣的處理包括 在所述最后一級(jí)，通過應(yīng)用預(yù)定義的預(yù)處理器對(duì)相應(yīng)的接口矩陣進(jìn)行因式分解。
13.根據(jù)權(quán)利要求12所述的方法，其中所述相應(yīng)的接口矩陣串行進(jìn)行因式分解。
14.根據(jù)權(quán)利要求12所述的方法，其中所述相應(yīng)的接口矩陣并行進(jìn)行因式分解。
15.根據(jù)權(quán)利要求11所述的方法，其中所述對(duì)最后一級(jí)接口矩陣的處理包括在所述最后一級(jí)通過應(yīng)用串行高質(zhì)量ILU因式分解對(duì)相應(yīng)的接口矩陣進(jìn)行因式分解。
16.根據(jù)權(quán)利要求11所述的方法，其中所述對(duì)最后一級(jí)接口矩陣的處理包括 在所述最后一級(jí)通過應(yīng)用并行迭代松弛的Block-Jacoby預(yù)處理器對(duì)相應(yīng)的接口矩陣進(jìn)行因式分解，并且對(duì)角塊進(jìn)行高質(zhì)量的ILU因式分解。
17.一種方法，包括(a)使用多級(jí)分割將初始矩陣分割成多個(gè)子矩陣；(b)對(duì)對(duì)角塊并行執(zhí)行截?cái)嘁蚴椒纸?，并且形成所述多個(gè)子矩陣的每一個(gè)的接口部分 的局部舒爾補(bǔ)；(c)由所述多個(gè)子矩陣的對(duì)角塊上的局部舒爾補(bǔ)和非對(duì)角塊上的所述多個(gè)子矩陣之間 的連接形成全局接口矩陣；(d)判定該全局接口矩陣是否足夠小，以滿足預(yù)定義的大小閾值；以及(e)當(dāng)判定該全局接口矩陣不是足夠小以滿足所述預(yù)定義的大小閾值時(shí)，使用多級(jí)分割將所述全局接口矩陣分割成多個(gè)第二子矩陣和為所述多個(gè)第二子矩陣重復(fù)操作(b)到 ⑷。
18.根據(jù)權(quán)利要求17所述的方法，進(jìn)一步包括 對(duì)所述多個(gè)子矩陣中的每個(gè)執(zhí)行局部重排序。
19.根據(jù)權(quán)利要求18所述的方法，其中所述局部重排序是在步驟(a)之后且步驟(b) 之前執(zhí)行的。
20.一種方法，包括使用Piay多級(jí)分割，對(duì)初始矩陣應(yīng)用不交迭的域分解以將所述初始矩陣分割成ρ個(gè) 部分，因此形成多個(gè)子矩陣預(yù)定義允許的最大遞歸級(jí)數(shù)；預(yù)定義最小大小閾值，其規(guī)定相對(duì)于所述初始矩陣大小的接口矩陣的允許的最小行數(shù)；遞歸執(zhí)行操作(a)到(d)(a)通過多個(gè)并行處理單元對(duì)所述多個(gè)子矩陣中的每一個(gè)并行執(zhí)行i)對(duì)角塊的并行 截?cái)嘁蚴椒纸?；和ii)為每個(gè)子矩陣的接口部分形成局部舒爾補(bǔ)；(b)由對(duì)角塊上的局部舒爾補(bǔ)和非對(duì)角塊上子矩陣間的連接隱式形成全局接口矩陣；(c)判定是否達(dá)到預(yù)定義的允許的最大遞歸級(jí)數(shù)或該全局接口矩陣的大小是否小于預(yù) 定義的最小大小閾值；(d)當(dāng)在操作(c)中判定沒有達(dá)到預(yù)定義的允許的最大遞歸級(jí)數(shù)并且該全局接口矩陣 的大小不小于預(yù)定義的最小大小閾值時(shí)，使用多級(jí)分割將該全局接口矩陣分割成另外的多 個(gè)子矩陣和針對(duì)所述另外的多個(gè)子矩陣重復(fù)操作(a)-(d)。
21.根據(jù)權(quán)利要求20所述的方法，進(jìn)一步包括當(dāng)在操作(c)中判定達(dá)到預(yù)定義的允許的最大遞歸級(jí)數(shù)或該全局接口矩陣的大小小 于預(yù)定義的最小大小閾值，則結(jié)束該遞歸處理。
22.根據(jù)權(quán)利要求21所述的方法，進(jìn)一步包括當(dāng)在操作(c)中判定達(dá)到預(yù)定義的允許的最大遞歸級(jí)數(shù)或該全局接口矩陣的大小小 于預(yù)定義的最小大小閾值，則對(duì)所述全局接口矩陣進(jìn)行因式分解。
23.根據(jù)權(quán)利要求20所述的方法，進(jìn)一步包括對(duì)所述多個(gè)子矩陣的每一個(gè)執(zhí)行局部重排序。
24.根據(jù)權(quán)利要求23所述的方法，其中所述重排序包括 首先對(duì)所述多個(gè)子矩陣的每一個(gè)的內(nèi)部行進(jìn)行重排序；以及 其次對(duì)所述多個(gè)子矩陣的接口行進(jìn)行重排序。
25.根據(jù)權(quán)利要求20所述的方法，其中所述隱式形成包括 通過所述多個(gè)并行處理單元的每一個(gè)存儲(chǔ)所述接口矩陣的相應(yīng)部分。
全文摘要
本發(fā)明提供一種采用預(yù)處理器的并行計(jì)算迭代求解器，預(yù)處理器使用用于求解線性方程組的并行計(jì)算被處理。因此，使用預(yù)處理算法用于大型稀疏線性方程組(如，代數(shù)方程、矩陣方程等)系統(tǒng)的并行迭代求解，如通常出現(xiàn)在對(duì)真實(shí)世界系統(tǒng)的基于計(jì)算機(jī)的三維(3D)建模(如，油或氣儲(chǔ)層的3D建模)中的線性方程組。提出一種新技術(shù)，其將多級(jí)預(yù)處理策略應(yīng)用到初始矩陣，初始矩陣被分割并轉(zhuǎn)換成對(duì)角塊加邊形式。提供一種導(dǎo)出預(yù)處理器的方法，用于線性方程組的并行迭代求解。特別是，并行計(jì)算迭代求解器可導(dǎo)出和/或應(yīng)用這樣的預(yù)處理器用于通過并行處理求解線性方程組。
文檔編號(hào)G06G7/48GK102138146SQ200980133946
公開日2011年7月27日 申請(qǐng)日期2009年7月17日 優(yōu)先權(quán)日2008年9月30日
發(fā)明者N·庫茲耐索娃, O·迪嚴(yán)科夫, S·寇舍萊夫, S·馬利亞索夫, V·普瑞伍尼科夫 申請(qǐng)人:?？松梨谏嫌窝芯抗?br>