CN102138146A

CN102138146A - 使用并行多级不完全因式分解求解储层模拟矩阵方程的方法

Info

Publication number: CN102138146A
Application number: CN200980133946.2A
Authority: CN
Inventors: O·迪严科夫; V·普瑞伍尼科夫; S·寇舍莱夫; N·库兹耐索娃; S·马利亚索夫
Original assignee: Exxon Production Research Co
Current assignee: ExxonMobil Upstream Research Co
Priority date: 2008-09-30
Filing date: 2009-07-17
Publication date: 2011-07-27
Also published as: CA2730149A1; WO2010039325A1; EP2350915A1; BRPI0919457A2; US20100082724A1; EP2350915A4

Abstract

本发明提供一种采用预处理器的并行计算迭代求解器，预处理器使用用于求解线性方程组的并行计算被处理。因此，使用预处理算法用于大型稀疏线性方程组(如，代数方程、矩阵方程等)系统的并行迭代求解，如通常出现在对真实世界系统的基于计算机的三维(3D)建模(如，油或气储层的3D建模)中的线性方程组。提出一种新技术，其将多级预处理策略应用到初始矩阵，初始矩阵被分割并转换成对角块加边形式。提供一种导出预处理器的方法，用于线性方程组的并行迭代求解。特别是，并行计算迭代求解器可导出和/或应用这样的预处理器用于通过并行处理求解线性方程组。

Description

使用并行多级不完全因式分解求解储层模拟矩阵方程的方法

相关申请的交叉参考

本申请要求2008年9月30号申请的美国临时专利申请61/101494的权益，其标题为METHOD FOR SOLVING RESERVOIR SIMULATION MATRIX EQUATION USING PARALLEL MULTI-LEVEL INCOMPLETE FACTORIZATION，该申请的全部内容通过引用包括在此。

技术领域

以下说明书一般涉及求解线性方程组的迭代求解器，更具体地涉及在高性能的并行计算系统上在求解线性方程组的并行迭代过程中执行预处理程序的系统和方法。

背景技术

在分析许多科学或工程应用时，通常需要同时求解大量的线性代数方程，这些方程可以以矩阵方程的形式表示如下：Ax＝b，(以下称为方程式(1))，其中A表示已知的n*n维的正方形系数矩阵，b表示通常称为“右手边”的已知n维向量，x表示通过求解该方程组得到的未知n维向量。已知各种技术用于求解这样的线性方程组。线性方程组通常在各种基于计算机的三维(3D)模拟或给定的真实世界系统的建模中遇到(并需要求解)。例如，地下含有碳氢化合物的储层(油或气储层)的现代三维(3D)模拟需要方程式(1)类型的线性代数系统的求解，一般地，在稀疏系数矩阵A中有数百万未知量以及数千万和数亿的非零元素。这些非零元素定义稀疏矩阵结构。

同样地，基于计算机的三维(3D)建模可以用来模拟真实世界系统，比如机械或电气系统(如可应用于汽车、飞机、轮船、潜艇、宇宙飞船等中)、人体(如整个人体或者人体的各部分，如重要器官等的建模)、天气形势以及各种其他要被建模的真实世界系统。通过这样的建模，可以分析或预测被建模系统的未来潜在性能。例如，提供给被建模系统的某些变化的条件对系统未来性能的影响可通过与基于计算机的模型的交互和分析进行评估。

作为一个例子来说，多孔介质中流体流动的建模是石油工业的主要焦点。不同的基于计算机的模型用于石油工业的不同领域，但大多数包括用偏微分方程系统或方程组(PDE’s)来描述该模型。一般地，这样的建模通常要求在给定网格上在空间和时间上对偏微分方程系统进行离散化，并为每个时间步长或时步执行计算直到达到规定的时间。在每个时步，离散方程被求解。通常离散方程是非线性的，并且求解过程是迭代的。非线性迭代方法的每个步骤通常包括非线性方程系统或方程组的线性化(如，雅可比构造)，求解该线性方程组以及用于计算下一个雅可比构造的性质计算。

图1示出一般的工作流程，其通常用于地下含碳氢化合物储层中随时间的流体流动的基于计算机的模拟(或建模)。内部循环101是求解非线性系统的迭代方法。再者，经过内部循环101的每次操作一般都包括非线性方程组的线性化(如，雅可比构造)11，求解线性系统(或方程组)12以及用于计算下一个雅可比(当循环返回到框11时)的性质计算13。外部循环102是时步循环(time step loop)。如所示的，针对每一个时步循环的边界条件可在框10中定义，然后在执行针对时步的内部循环101后，为时步计算的结果可在框14中输出(如结果可存储到数据存储介质或/和提供给软件应用程序以便生成显示，该显示表示为相应时步建模的地下含碳氢化合物储层中的流体流动)。如上述所示，除了地下含碳氢化合物储层中流体流动的建模之外，真实世界系统的基于计算机的3D建模可以以类似方式执行，即可采用求解线性方程组的迭代方法(如图1的框12中所示)。

线性方程组的求解是一种计算强度非常高的任务，因此需要高效的算法。有两类一般的线性求解器：1)直接法和2)迭代法。所谓的“直接法”基于高斯消元法(Gaussian elimination)，该方法中对矩阵A进行因式分解，并表示为下三角形矩阵L和上三角形矩阵U(因式)的积：A＝LU(以下称为方程式(2))。但对大稀疏矩阵A而言，计算三角形矩阵L和U非常耗费时间，且那些因子中的非零元素的数量会很大，因此它们可能不适合甚至现代高性能计算机的存储器。

迭代法是基于简单而通常低廉操作的反复应用，如矩阵向量积，其提供具有给定精确度的近似解。通常，对于出现在科学或工程应用中的方程式(1)类型的线性代数问题而言，系数矩阵的性质导致收敛到一个解的大量迭代。

为了减少迭代次数，并因此减少运用迭代方法求解矩阵方程的计算开销，经常使用预处理技术，在该技术中方程式(1)类型的初始矩阵方程与合适的预处理矩阵M(也可简称为预处理器(preconditioner))相乘，如：M^-1Ax＝M^-1b(以下称为方程式(3))。这里，M^-1表示矩阵M的逆矩阵。采用不同的预处理方法(矩阵M)能极大减少以足够高的准确度计算方程式(1)的合适解的运算开销。预处理技术的主要例子是代数多网格方法(algebraic multi-grid methods)和不完全上下三角因式分解(incomplete lower-upper factorizations)。

在第一种方法(即多网格方法)中，构造一系列减维的系数矩阵，且建立从较精细维到较粗维的一些数据转移方法。然后，近似求解矩阵方程(1)(所谓的“平滑”)，计算余项r＝Ax-b，并将得到的向量r传递到较粗维(所谓的“约束”)。因此可在较粗维中近似求解类似方程式(1)的方程式，计算余项并传递到较粗维等等。当在最粗维上计算该问题之后，此粗解会被传递回原始维(所谓的“延长”)，以得到将被添加在原始精细维上的近似解的缺陷。

另外一个预处理技术的例子是不完全上下三角因式分解(ILU型)，该方法取代了完全因式分解(如方程式(2)所示)，计算稀疏因子L、U使得其乘积近似等于初始系数矩阵：A≈LU(以下称为方程式(4))。

上述的两种预处理技术本质上是有顺序的，且不能直接应用在并行处理计算机上。随着出现在科学和工程应用中的代数问题的维数增加，对适合并行处理计算机的求解方法的需求变得越来越重要。因此，研发高效的并行线性求解器成为日益重要的任务，这对许多3D建模应用诸如石油储层建模而言尤其重要。尽管很多不同的用大稀疏系数矩阵求解矩阵方程的方法有实质进步，例如多网格或直接求解器，但是在过去的数十年中，采用基于不完全上下三角因式分解的预处理的迭代方法仍然是求解大型稀疏线性系统最普及的方法。如上所述，这些预处理技术实质上是有顺序的，且不能直接应用在并行处理计算机上。

近年来在科学界，一种新型的使用多级块ILU因式分解技术的并行预处理策略被开发用于求解大型稀疏线性系统。这个新方法的大体思想是对未知量和对应的方程式重排序，并且将初始矩阵拆分成2*2的块结构，该拆分方式使得第一对角块成为块对角矩阵。该块矩阵会被并行进行因式分解。通过消去因式分解后的矩阵块形成舒尔补(Schur complement)之后，反复执行该步骤以得到舒尔补。新方法的效率依赖于初始矩阵和舒尔补被拆分成块的方式。在传统的方法中，多级因式分解是基于稀疏矩阵结构原始图的多着色或块独立集分割。这些技术进一步详述在以下文献中：a)C.Shen，J.Zhang和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，pARMS：A parallel version of the algebraic recuisive multilevel solver，Number.Linear Algebra Appl.，10(2003)，pp.485-509，该文献的公开内容通过引用包括进本发明。这些方法的缺点是它们改变矩阵的初始排序，这在很多情况下会导致预处理器的质量更糟，或迭代求解器的收敛速度更慢。另一个缺点是这样的并行重排序结构不能良好地扩展，即随着处理单元(处理器)数量的增加，该结构的质量和效率明显下降。

另外一类基于ILU因式分解的平行预处理策略使用源于域分解(domain decomposition)的思想。给定大型稀疏线性方程式系统(1)，首先使用分割软件(例如G.Karypis和V.Kumar在文献METIS：Unstructured Graph Partitioning and Sparse Matrix Ordering System，version 4.0，1998年九月中叙述的)，该矩阵A被分割成给定数量的子矩阵p，几乎每个子矩阵的行数都相同，且子矩阵间的连接很少。在分割步骤之后，进行局部重排序，首先要对每个子矩阵的内部行进行排序，然后是子矩阵的接口行，即那些与其他子矩阵有连接的行。然后，已分割的以及序列改变的初始矩阵A可表示在如下的加边对角块(BBD)形式中：

其中Q是具有局部置换矩阵Q1的置换矩阵，矩阵B是全局接口矩阵，

其包含所有接口行和所有子矩阵的外部连接，矩阵B具有如下结构：

这样的矩阵表示形式广泛应用在科学计算中，参看文献例如：a)D.Hyosom和A，Pothen，A scalable parallel algorithm for incomplete factor preconditioning，SIAM J.Sci.Comput，2001年22期，2194-2215页(以下称为文献“Hysom”；)b)G.Karypis和V.Kumar.Parallel Threshold-based ILU Factorization。AHPCRC，明尼阿波利斯，MN55455，Technical Report#96-061(以下称为文献“Karypis”)；以及c)Y.Saad，Iterative Methods for Sparse Linear System，2nd ed，SIAM，2003年，费城(以下称为文献“Saad”)。

基于BBD格式的并行预处理的下一步是因式分解过程。有几种方法进行因式分解。一种方法在例如：文献Hysom和Karypis中考虑。在文献Hysom中，首先，对内部行并行进行因式分解。如果对于一些处理元件而言没有更低序的连接，那么也对边缘行进行因式分解。否则，处理单元等待要接收到行结构和更低序连接的值，只有在此之后，才对边缘行进行因式分解。因此，该方案的时间平衡性不是很好，因为有更高指标的处理单元要等待来自具有更低指标的相邻处理单元的已因式分解的边缘行。因此，随着处理单元数量的增加，该方法的可扩展性就变糟糕了。

在文献Karypis中，BBD格式矩阵的上部的因式分解是并行执行的，而下矩形部分

的因式分解是使用块B的并行最大独立组重排序执行的，该操作可以执行多次。之后，经过修改的不完全因式分解过程的并行版本应用到矩阵的整个下部分。此外，使用独立组重排序的矩阵的部分置换可能导致收敛和可扩展性更糟糕。

美国专利5655137号(以下称为“137”专利)中描述另外一种方法，该专利的公开内容通过引用包括进本发明。一般来说，该“137专利”建议对对角线块A₁到A_p以

(不完全Cholesky因式分解)的形式并行进行因式分解，并使用这些局部因式分解计算矩阵B的舒尔补。该方法仅可以应用于对称正定矩阵。

文献Saad中描述的不同方法应用ILU因式分解的截断变型来以如下的方式对全部子矩阵进行因式分解，子矩阵包括边缘行，该方式使得对每个第i个子矩阵计算局部舒尔补S_i，且对这些局部舒尔补求和就得到全局舒尔补。因此，得到的舒尔补矩阵为如下形式：

该类方法有两个主要缺陷。首先是当矩阵部分数增加时，舒尔补S的大小会急剧增大。第二个问题是矩阵S的高效因式分解问题。

存在对改进的迭代求解方法的需求，该方法能够并行处理多级不完全因式分解。

发明内容

本发明针对的是采用并行计算迭代求解器的方法和系统。因此，本发明的实施例一般涉及并行高效能计算领域。本发明的实施例更具体涉及预处理算法，这些算法适用于大型稀疏线性方程组(如，代数方程、矩阵方程等)系统的并行迭代求解，如通常出现在基于计算机的对真实世界系统进行三维(3D)建模中的线性方程组(如油或气储层的三维建模)。

根据某些实施例，提出一种新技术，将多级预处理策略应用到初始矩阵，该初始矩阵被分割并转换成加边对角块形式。

根据一个实施例，提供一种导出预处理器的方法，预处理器用于线性方程组的并行迭代求解。特别是，并行计算迭代求解器可通过并行处理导出或应用这样的求解器，用于求解线性方程组。如本文进一步论述的，这种并行计算迭代求解器提高通过并行执行各种操作求解这种线性方程组的计算效率。

根据一个实施例，不交迭的域分解(non-overlapping domain decomposition)被应用到初始矩阵，通过使用p-way(p路)多级分割将初始图分割成p个部分。然后应用局部重排序。在局部重排序中，根据一个实施例，首先对每个子矩阵的内部行进行排序，然后才对其“接口”行(即那些与其他子矩阵有连接的行)进行排序。因此，局部第i个子矩阵将具有如下形式：

A_{i} = [\begin{matrix} A_{ii}^{I} & A_{ii}^{IB} \\ A_{ii}^{BI} & A_{ii}^{B} \end{matrix}] + \underset{j &NotEqual; i}{Σ} A_{ij} = [\begin{matrix} A_{i} & F_{i} \\ C_{i} & B_{i} \end{matrix}] + \underset{i &NotEqual; j}{Σ} A_{ij} = A_{ii} + \underset{i &NotEqual; j}{Σ} A_{ij}

(以下为方程(5))，

其中，A_i是内部行之间有连接的矩阵，F_i和C_i是内部行和接口行之间有连接的矩阵，B_i是接口行之间有连接的矩阵，且A_ij是子矩阵i和j之间有连接的矩阵。应当明白，矩阵A_ii对应于第i个子矩阵的对角块。

在一个实施例中，该过程对对角块执行并行截断因式分解，为每个子矩阵B_i的接口部分形成局部舒尔补。全局接口矩阵由对角块的局部舒尔补和在非对角块的子矩阵之间的连接形的。从结构来看，结果矩阵具有块结构。

然后上述过程反复执行，从接口矩阵的重新分割开始直到接口矩阵足够小(如，与一个预定义大小的最大值相比)。在某些实施例中执行接口矩阵的重新分割以使得子矩阵之间的连接数目最少。当确定该接口矩阵的维数足够小时，就可直接或使用并行迭代方法(如，Block-Jacoby)进行因式分解。

根据某些实施例，该算法是反复(递归)执行以上所述步骤，同时隐式形成大小大于某个预定义的大小阈值的接口矩阵或目前级数小于允许最大级数。同时，上述步骤在下级应用之前，通过某个分割器

(如这里进一步说明的并行多级分割器)对接口矩阵进行重新分割。此外，在并行截断因式分解之前使用局部对角缩放，以在某些实施例中改进局部因式分解对角块的数值性质。如以下所述，很多更复杂的局部重排序可应用在某些实施例中。大体来说，基于对已知算法序列的反复(递归)应用，某个实施例的算法与一个通用的框架中的算法(它们是本领域中广泛已知的)合并，以形成减维的矩阵序列(多级方法)。

上述使用多级方法的方法可以作为预处理器应用于迭代求解器中。此外，可以应用特定的局部缩放和局部重排序算法，以改进预处理器的质量。该算法可以应用在共享存储和分布式存储并行体系结构中。

前面已经相当宽泛地概括了本发明的特征和技术优点，以便可更好地理解下面本发明的详细说明。本发明的额外特征和优点将在下面说明，这形成本发明权利要求的主题。本领域技术人员应该理解，公开的概念和特定实施例可易于用作修改或设计执行本发明相同目的的其他结构的基础。本领域技术人员还应该了解这样的等效构造不偏离权利要求中限定的本发明的精神和范畴。被认为是本发明特征的新颖特征，关于其组织和运算的方法，以及本发明的进一步的目的和优点可更好地从下面的说明中结合附图理解。然而应该明确理解，提供的每个图的目的仅仅是例示和说明，而非用作限制本发明的定义。

附图说明

为了更完整地理解本发明，参考结合附图的以下说明，其中：

图1示出一般工作流程，其通常用于对地下含碳氢化合物储层中流体流动随时间的基于计算机的模拟(或建模)。

图2示出根据本发明的一个实施例实现并行计算迭代求解器的基于计算机的示例性系统的框图。

图3示出根据本发明的一个实施例实现并行计算迭代求解器的基于计算机的另一个示例性系统框图。

图4示出根据本发明某些实施例的示例性计算机系统，该系统能实现全部的或部分的并行计算迭代求解器。

具体实施方式

本发明的实施例一般涉及并行高性能计算领域。本发明的实施例更具体地针对预处理算法，其适合于并行迭代求解大型稀疏系统的线性方程组(如，代数方程、矩阵方程等)，例如通常出现在基于计算机的对真实世界系统的三维(3D)建模(如，油或气储层的三维建模)中的线性方程组。

根据某些实施例，提出一种新技术，用于将多级次预处理策略应用到被分割并且转换到加边对角块形式的初始矩阵。

根据一实施例，图2显示导出用于线性方程组的并行迭代求解的预处理器的方法。图2示出根据本发明一实施例的基于计算机的示例性系统200的框图。如所示的，系统200包括基于处理器的计算机221，例如个人计算机(PC)、膝上型(即便携式)计算机、服务器计算机、工作站计算机、多处理器计算机、计算机集群等。此外，并行迭代求解器(如应用软件)222运行在计算机221上。计算机221可以是任何基于处理器的设备，本文将进一步说明，该设备能够运行并行迭代求解器222。优选的，计算机221是包括多个处理器的多处理器系统，其能执行并行迭代求解器222的并行操作。虽然出于方便说明，在图2中将并行迭代求解器222示为运行在计算机221上，但是应当明白这种求解器222可驻留和/或局部运行在计算机221或远程计算机(如服务器计算机)上，计算机221通过通信网络例如局域网(LAN)、因特网或其他广域网(WAN)等通信耦合到远程计算机。此外，应当明白计算机221可包括多个集群的或分布式的计算设备(如服务器)，如本领域已知的，并行计算求解器222可存储和/或运行在该设备上。

和许多常规的基于计算机的迭代求解器一样，并行迭代求解器222包括存储在计算机可读介质上的计算机可执行的软件代码，该介质是计算机221的(多个)处理器可读的，当这种处理器运行所述软件代码时使得计算机221执行本文进一步说明的针对这种并行迭代求解器222的各种操作。并行迭代求解器222是可操作的，以使用迭代过程来求解线性方程组，其中迭代过程的各部分被并行执行(例如在计算机221的多个处理器上)。如上所述，迭代求解器通常用于基于计算机的三维(3D)建模。例如，并行迭代求解器222可用在(图1中)常规工作流程的操作框12中，用于地下含碳氢化合物储层之中流体流动的基于计算机的三维建模。在图2所示的例子中，模型223(例如，包含各种有关要建模的真实世界系统的信息，如有关地下含碳氢化合物储层中要建模的随时间的流体流动的信息)存储在数据存储器224上，存储器224通信耦合到计算机221。数据存储器224可包括硬盘、光盘、磁盘和/或可操作用于存储数据的其他计算机可读的数据存储介质。

正如许多用于基于计算机的三维(3D)建模的常规迭代求解器一样，并行迭代求解器222可操作用于接收模型信息223，并执行迭代方法用于求解线性方程组以产生基于计算机的三维(3D)模型，例如，地下含碳氢化合物存储层中随着时间的流体流动的模型。如本文进一步论述的，并行迭代求解器222可通过并行执行各种操作提高求解这种线性方程组的计算效率。根据一实施例，并行迭代求解器能执行下文论述的操作201-209。

如在框201中所示，将不交迭的域分解(non-overlapping domain decomposition)应用于初始矩阵，以使用p-way多级分割将初始图分割成p个部分。应该认识到，该分割可能被认为是算法外部的操作，因为对原始数据的分割通常是任何并行计算的必要操作。

在框202中，应用局部重排序。如在子框203中所示，首先对每个子矩阵的内部行进行排序，然后在子框204中对子矩阵的“接口”行(即那些与其他子矩阵有连接的行)进行排序。因此，第i个局部子矩阵具有上述方程5的形式。除了(或替代)局部重排序，在某些实施例中也可执行局部缩放算法以改进子矩阵的数值性质，因此，提高独立截断因式分解的质量。在某些实施例中，框202中的局部重排序是该算法的一个选择，其可能在某些实施中被省略。局部重排序不仅仅是简单的重排序来以给定的自然顺序先移动内部结点，最后移动接口结点，而且实施为更复杂的算法，如图多级方式，以最小化重排序的对角块的轮廓(profile)，如下文进一步说明的。

在框205中，该过程执行对角块的并行截断因式分解，形成每个子矩阵B_i的接口部分的局部舒尔补。

在框206中，由对角块的局部舒尔补和非对角块的子矩阵之间的连接形成全局接口矩阵(参见方程(4))。通过构造，结果矩阵(合成的矩阵)具有块结构。应当认识到，在某些实施例中，全局接口矩阵没有在框206中显式形成(这可能是开销非常大的运算)，而是用于并行处理的多个处理单元的每一个可能存储接口矩阵的其相应部分。这样，在某些实施例中，可能隐式而不是显式形成该全局接口矩阵。

反复应用所有的子框202-206，从对接口矩阵的重新分割(在子框208)开始直到该接口矩阵足够小。术语“足够小”在该实施例中理解为以下意义。该方法有两个参数，它们限制应用多级算法：1)max levels规定允许的最大级数，2)min size是一个阈值，其确定相对于初始矩阵的大小而言以接口矩阵行数表示的允许的最小大小。根据本实施例，当递归级达到允许的最大级数值时或该接口矩阵的大小变得比min size和初始矩阵的大小的乘积值还小时，则该递归过程终止且执行最低级的预处理。

因此，在框207中，无论接口矩阵是否足够小，都要做出判定。如果判定接口矩阵不是足够小，那么操作推进至框208，以对接口矩阵进行重新分割(如在框201对初始矩阵进行分割那样)，并在框202-206中反复对接口矩阵执行重新分割。在框208中的重新分割很重要，从而最小化子矩阵之间的连接数目。当在框207中判定该接口矩阵的维数足够小时，就可在框209中直接或使用并行迭代法(例如Block-Jacoby)对其进行因式分解。

使用多级法的该方法作为预处理器应用于迭代求解器。此外，可应用特定的局部缩放和局部重排序算法，以提高预处理器的质量。该算法可以同时应用在共享存储和分布式存储并行体系结构中。

根据本发明的实施例，图3示出示例性的基于计算机的系统300的另一框图。如以上关于图2所讨论的，系统300也包括基于处理器的计算机221，并行迭代求解器的一个示例性实施例，如图3示出的并行迭代求解器222A运行在计算机221上，执行下文所述操作。根据该实施例，并行迭代求解器222A使用多级法，如下文关于框301-307的说明。

传统上而言，该多级预处理器MLPrec包括以下参数：MLPrec(l，A，Prec1，Prec2，l_max，τ)，其中l是当前级数，A是不得不进行因式分解的矩阵，Prec1是用于独立子矩阵A_ii＝L_iU_i的因式分解的预处理器，Prec2是对最后一级的舒尔补S进行因式分解的预处理器，l_max是允许的最大级数，且τ是阈值，用于定义相对于矩阵A的大小的允许的最小值S。

在这个示例实施例的操作中，该并行迭代求解器于框301中以MLPrec(0，A，Prec1，Prec2，l_max，τ)开始。在框302中，该迭代求解器判定是否满足|S|＞τ·|A|以及l＜l_max。当在框302中判定满足|S|＞τ·|A|以及l＜l_max时，在框303中为修正的舒尔补矩阵S’：MLPrec(l+1，S’，Prec1，Prec2，l_max，τ)递归重复上述(图2的)并行方法。例如，这样的递归重复操作可包括在子框304中对修正的舒尔补矩阵进行分割(如图2中的框208)，在子框305中对分割的舒尔补子矩阵进行局部重排序(如图2中的框202)，并在子框306中对对角块执行并行截断因式分解(如图2中的框205)。因此，在一个实施例中，修正的矩阵S’是对矩阵S应用某个分割器之后得到的(如图2中的框208)，其试图最小化矩阵S中的连接数目。该分割器可以和在第一级(即图2的框201中)用于初始矩阵分割的分割器一样，或者分割器在某些实施中有所不同。

当在框302中判定满足|S|＜τ·|A|或l＞l_max，在框307中使用预处理器Prec2在最后一级对舒尔补矩阵S_i进行因式分解。如本文进一步说明的，作为例子，可以使用针对很小S_i的高质量串行(或连续的)ILU预处理器或带对对角块进行ILU因式分解的并行块雅各比(Jacoby)预处理器。

为了提高预处理器的质量，某些实施例也使用两种额外的局部预处理技术。第一种是从矩阵A₁₁到A_pp的局部缩放。第二种技术是特殊的局部重排序，其最后移动接口行，接着以图多级方式对内部行进行排序，最小化该重排序的对角块的轮廓A_ii＝Q_iA_iiQ_i ^T。

除了局部重排序，也可在某些实施例中执行局部缩放算法，以改进子矩阵的数值性质，并因此提高独立截断因式分解的质量。进一步，局部重排序并不要求用于所有的实施例，但可作为实施该算法的实施例的一种选择。局部重排序可以不仅仅包括简单的重排序，按照给定的自然顺序先移动内部结点最后移动接口结点，还能作为更复杂的算法，如上述最小化重排序对角块的轮廓的图多级方式。

因此，根据本发明的某些实施例，并行迭代求解器使用基于域分解法的多级方法，将初始矩阵转换到2乘2的块形式。进一步，在某些实施例中，该并行迭代求解器使用局部对角块的ILU型因式分解的截断变型，以得到作为局部舒尔补矩阵之和的全局舒尔补矩阵。而且在某些实施例中，在针对得到的全局舒尔补矩阵重复多级过程之前，并行迭代求解器对得到的全局舒尔补矩阵进行重新分割，以最小化已分割矩阵的连接数目。在该多级方法的最后一级，在某些实施例中，该并行迭代求解器使用串行ILU预处理器或并行块雅各比预处理器。此外，在某些实施例中，该并行迭代求解器应用局部缩放和减少局部重排序的矩阵轮廓的特殊变型。

并行迭代求解器的一个说明性实施例将在下面通过运行在具有若干独立处理器的分布式存储体系结构上的并行求解方法示例进一步阐述。实施例可能也应用于共享存储和混合型的体系结构。可用于共享存储体系结构(SMP)和混合体系结构的算法，与下文针对说明性实施例描述的示例性算法很相似，除了在某些实施细节上不同。这些细节能够很容易被普通技术人员认识到(这些将在下文适用的地方分别论述)。

这个说明性实施例的并行多级预处理器是基于不完全因式分解，并且可以在下面简称为PMLILU。

预处理器构造。在这个说明性实施例中，PMLILU预处理器是基于域分解法的不交迭(overlapping)形式。域分解法假定整个问题的解决方案可通过在子问题之间接口上的解决方案聚合的特定过程，由以某种方式分解的子问题的解决方案得到。

初始矩阵A的稀疏结构图G_A被分割成给定的p个不交迭的子图G_i，使得

G_k∩G_m＝φ：k≠m。

这样的分割对应于矩阵A逐行分割成p个子矩阵：

其中是行带，可以以块的形式表示如下：

分割成行带对应于处理单元间的矩阵分布。应当注意，在这个说明性实施例中，向量是按相同的方式分布的，即对应于子图G_i的元素的向量的那些元素存储在行带

也存储在其中的相同处理单元中。

这种分割在下面表示为

其中，l是级数(有时因简明性被省略)，p是部分数。第i部分的大小(行数)表示为N_i，该部分对于第一行的偏移量(以行表示)记为O_i。因此，

一般的矩阵分割块形式可以写为：

在下面的论述中，为了简洁使用该矩阵符号。术语“矩阵行”通常用于代替更传统的术语“图结点”，但是这两个术语都可以在以下说明中交替地使用。因此，图结点对应于矩阵行，而图边缘对应于矩阵非零非对角项，它们是行之间的连接。符号

意思是第k个矩阵行属于第i个行带。标准的图符号m∈adj(k)用来说明a_km≠0，意思是在第k和第m个矩阵行之间存在着连接。在下文的说明中，术语“部分”对应于术语“行带”。术语“块”用于定义对应于分割的行带的一部分。具体是，A_ii对应于对角块，其中

一般来说，根据该说明性实施例，该预处理器构造算法的主要步骤可以简洁表述如下：

1、矩阵被(连续或者并行)分割成给定的p个部分(如图2中框201所示)。在这样的分割之后，该矩阵在处理器中作为行带分布。

2、分割之后，

的行划分为两组：1)内部行，即与其他部分的行没有连接的行，2)接口(边缘)行，与其他部分有连接的行。应用局部重排序(如图2的框202中所示)到每个行带，从而先移动内部行，最后移动接口结点。对每个行带独立应用该重排序(并行的)。

3、计算对对角块的并行截断因式分解，并且计算针对对应接口对角块的舒尔补(如图2的框205中所示)。

4、形成接口矩阵(如图2的框206中所示)。在这个说明性实施例中，接口矩阵包括接口对角矩阵的舒尔补和非对角连接矩阵的舒尔补。

a、如果判定(如图2的框207中所示)接口矩阵的大小足够小或达到允许的最大级数时，就可对接口矩阵进行因式分解(如图2的框209中所示)。

b、否则，上面步骤1-4中所述的同样算法以和应用到初始矩阵相同的方式应用到接口矩阵。应当注意，在构造过程的步骤1中，接口矩阵应该再次被分割，以最小化各部分之间的连接数(例如在图2的框208中接口矩阵被分割，然后反复执行图2的框202-207中的操作，用于处理已分割的接口矩阵)。

接口矩阵在最低(最后)级的因式分解可作为接口矩阵的完全ILU因式分解(它是更健壮的变型)串行执行，或使用迭代松弛的Block Jacoby方法并行执行对对角块进行ILU因式分解。因此，在某些实施例中，一些相对小的接口矩阵允许某个串行工作，但串行工作的优点就是随着部分数的增加达到一个稳定的迭代数。

应当注意，整个并行求解过程可从初始矩阵的分割开始(如图2的框201中所示)，任何算法(比如预处理器、迭代方法等等)都可使用该求解过程。因此，初始分割(图2的框201中)是相对于预处理器的外部操作。因此，该说明性实施例的PMLILU具有一个已分割(并且已发布的)矩阵作为输入参数。这在下面以算法1(用于预处理器的构造)的示例性伪代码进行说明：

已知一个矩阵A，右手边向量b，未知向量x，部分数p

//应用外部分割

//调用并行多级ILU。

在该说明性实施例中，根据以下的算法2的示例性伪代码，可以写出并行多级ILU算法(PMLILU)：

定义的算法：Truncated_ILU，Last_level_Prec，Local_Reordering，

Local_Scaling，Partitioner_IM(接口矩阵分割器)

参数：max_levels，min_size_prc

应当注意，上述算法2针对PMLILU使用的任何类型的基本算法定义，例如Truncated_ILU，Last_level_Prec，Local_Reordering，Local_Scaling，Partitioner。可以选择任何合适的算法，并使用该算法而不是PMLILU。

下面考虑根据该说明性实施例的该算法构造的步骤，并且与所考虑步骤相关的实施细节将针对该说明性实施例进一步说明。

矩阵分割。如上所述，系统的初始分割和分布在预处理器构造算法的外部执行如下：

对初始分割(“IPT”)而言，可以使用高质量的多级法，其与来自知名软件包METIS(如G.Karypis和V.Kumar所著的METIS：Unstructured Graph Partitioning and Sparse Matrix Ordering System，Version 4.0，1998年9月，所述，该书的公开内容通过引用包括进本申请)的方法相类似。下面将进一步说明接口矩阵分割。

还应该注意，通常任何分割器改变初始矩阵的顺序。取决于基本算法和质量约束，该已分割矩阵可以写成

因此，上述算法2在输入端得到已置换、已分割并且已分布的矩阵。

对SMP体系结构而言，在分布式的数据结构中存储行带

也是非常有优势的，它允许存储访问的成本明显下降。对于这点，

应该被并行分配，然后这允许任意线程使用最优位于存储区的矩阵部分。在允许将特殊线程绑定到某些处理单元的那些共享存储体系结构中，该绑定过程在性能上可提供附加改进。

局部重排序。在分割之后，行带

的行要被划分为两个子集：

1)内部行的集合

它们是与其他部分的行没连接的行：

以及

2)接口(边缘)行的集合，它们是与其他部分的行有连接的行：

{k &Element; A_{i^{*}}^{B} : &Exists; m &Element; adj (k), m &Element; A_{j^{*}}, j &NotEqual; i} . .

应用局部重排序，以首先计算内部行，然后是接口行：

P_{i} A_{ii} P_{i}^{T} = (\begin{matrix} A_{ii}^{I} & B_{ii}^{IB} \\ A_{ii}^{BI} & A_{ii}^{B} \end{matrix}) = (\begin{matrix} A_{i} & F_{i} \\ C_{i} & B_{i} \end{matrix})

由于该操作的局部性，所以在该说明性实施例中并行执行该操作。在对角块局部置换之后，集合所有来自相邻的处理器的局部置换向量，从而对非对角矩阵A_i进行置换(如果Aij≠0)。该说明性实施例的局部重排算法的总框架可以用伪代码写作如下(称作是算法3)：

使用各种局部重排序算法是可能的，但是为了简洁使用一种自然排序，诸如在该说明性实施例的以下示例性伪代码中(称作是算法4)：

因此，在应用局部置换之后，该矩阵可以写成如下：

[\begin{matrix} A_{1} & F_{1} \\ C_{1} & B_{1} & A_{12} & \cdot \cdot \cdot & A_{1 p} \\ A_{2} & F_{2} & \cdot \cdot \cdot \\ A_{21} & C_{2} & B_{2} & A_{2 p} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ A_{p} & F_{p} \\ A_{p 1} & A_{p 2} & \cdot \cdot \cdot & C_{p} & B_{p} \end{matrix}]

重排列所有接口(边缘)矩阵块B_i和A_ij，一直到矩阵的最后，得到矩阵划分的加边对角块形式(BBD)：

[\begin{matrix} A_{1} & F_{1} \\ A_{2} & F_{2} \\ \cdot \cdot \cdot & \cdot \cdot \cdot \\ A_{p} & F_{p} \\ C_{1} & B_{1} & A_{12} & \cdot \cdot \cdot & A_{1 p} \\ C_{2} & A_{21} & B_{2} & \cdot \cdot \cdot & A_{2 p} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ C_{p} & A_{p 1} & A_{p 2} & \cdot \cdot \cdot & B_{p} \end{matrix}],

其中，矩阵的右下角是接口矩阵，该接口矩阵集合矩阵各部分之间的所有连接：

(\begin{matrix} B_{1} & A_{12} & \cdot \cdot \cdot & A_{1 p} \\ A_{21} & B_{2} & \cdot \cdot \cdot & A_{2 p} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ A_{p 1} & A_{p 2} & \cdot \cdot \cdot & B_{p} \end{matrix})

根据该说明性实施例，接口矩阵的处理过程在下面进一步说明。

局部缩放。缩放能显著地提高预处理器的质量，因此，能提高该迭代求解器的整体性能。对于由每个网格单元具有若干未知量(自由度)的偏微分方程(PDE)的离散化得来的矩阵尤其如此。一般来说，应用一些考虑因素，缩放算法对两个对角矩阵D^L和D^R进行计算，这提高某些矩阵的缩放性质(如，均衡对角项或行/列模范数的大小)，这通常可使因式分解更稳定。应用全局缩放可导致增加处理元件之间的额外通信开销，而对一个部分的对角矩阵应用局部缩放只要求列缩放矩阵D^R的部分集合，而质量上不会有显著的损失。

应当注意，在该说明性实施例中，缩放应用在一个部分的整个对角块：

局部缩放算法可以记为如下：

并行截断因式分解。在该说明性实施例中，算法的下一步骤是对角块的并行截断因式分解与计算相应的接口对角块的舒尔补：

//Parallel truncated factorization

in_parallel i＝l:p{

//truncated factorization

}

ILU因式分解的截断(受限的)变型要计算不完全的因子并且模拟舒尔补，并可与2003年Y.Saad所著的Iterative Methods for Sparse Linear Systems，2nd ed，SIAM，Philadelphia，2003中描述的操作类似地实施，该书的公开内容通过引用并入本文。

因此，已因式分解的对角块具有以下结构：

S_i＝B_i-C_i(L_iU_i)^-1F_i

接口矩阵处理。在该说明性实施例中，算法的最后一步骤是对接口矩阵的处理。在执行上述并行截断因式分解之后，接口矩阵可以记为如下：

S_{L} = (\begin{matrix} S_{1} & A_{12} & \cdot \cdot \cdot & A_{1 p} \\ A_{21} & S_{2} & \cdot \cdot \cdot & A_{2 p} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ A_{p 1} & A_{p 2} & \cdot \cdot \cdot & S_{p} \end{matrix}),

其中，S_i是第i个接口矩阵的舒尔补。现在校验该矩阵的大小(行数)，如果其足够小，则应用最后一级因式分解；否则，反复递归执行针对已重新分割矩阵S_i的整个过程，例如在下面的示例性伪代码所示：

应当注意，一般来说，根据该说明性实施例，并不需要显式形成这个矩阵，两种情况除外：

1、如果使用串行分割对接口矩阵进行分割，那么执行显式形成这个矩阵的操作以集合其图；

2、如果针对最后一级因式分解定义串行预处理，那么执行显式形成这个矩阵的操作以将其集合作为一个整体。

一般来说，该接口矩阵分割器，Partitioner_IM可与初始分割器(在图2的框201中使用的)不同。如果使用相同结构的矩阵求解一系列线性代数问题，如在建模时间相关的问题，该初始分割器可以是串行的，并且也可以在整个多时步模拟中仅仅使用几次(或甚至一次)。与此同时，Partitioner_IM应该是并行的以避免串行分割(尽管该变型也是可能的，可用在某些实施中)的接口矩阵图聚集到一起。

最后一级因式分解有两个主要的变型，可以按照本说明性实施例来使用：

1、纯串行ILU；和

2、具有对角块中的ILU的迭代松弛的Block-Jacoby(IRBJILU)。

因此，根据某些实施例，该算法有利地使用接口矩阵的并行多级分割，以避免该接口矩阵明显地形成在主处理单元上，而串行多级分割这样要求。在处理接口矩阵的最后一级时，相应的接口矩阵可通过应用预定义的预处理器进行串行因式分解或并行因式分解。可用于接口矩阵最后一级的这种处理的可能变型包括：高质量的串行ILU因式分解或具有对角块的高质量ILU因式分解的并行迭代松弛的Block-Jacoby预处理器。

在大多数数值实验中，第一种变型(即纯串行ILU)使并行迭代求解器的整体性能更好，并保持和对应串行变型的迭代数量几乎相同的收敛所要求的迭代数量；而第二种变型(即IRBJILU)在矩阵部分数增加时，降低收敛。

现在考虑按照该说明性实施例，每个处理单元可存储哪些预处理器数据。对每级l而言，第i个处理器存储

其中

来自接口矩阵分割器的某些聚集信息，是第i个处理器所需的(置换向量、分割阵列以及其他例子中更多)。此外，主处理器存储最后一级因式分解的预处理矩阵M_i。应当注意，在该说明性实施例中，在因式分解过程中使用接口矩阵之后没有必要保持接口矩阵。

并行预处理器求解方法。在该说明性实施例的迭代方法的每次迭代中，用通过ILU型因式分解得到的预处理器解决线性代数问题。从结构来看，预处理器在该说明性实施例中表示为上三角矩阵和下三角矩阵的乘积；且可以将求解过程定义为向前和向后替代的方式：

LUt＝s

Lw＝s；Ut＝w

为使提出的方法在该说明性实施例中有效，需要开发该过程的有效并行公式。在下面的说明中，实际是下三角求解的向前替代表示为Lsolve，向后替代则表示为U solve。

三角求解器的并行公式采用由因式分解过程产生的因子L、U的多级结构。它实现多级求解法的并行变型。根据初始的分割使向量t、s以及w进行拆分：

t = (\begin{matrix} t_{1} \\ t_{2} \\ \cdot \cdot \cdot \\ t_{p} \end{matrix}),

s = (\begin{matrix} s_{1} \\ s_{2} \\ \cdot \cdot \cdot \\ s_{p} \end{matrix}),

w = (\begin{matrix} w_{1} \\ w_{2} \\ \cdot \cdot \cdot \\ w_{p} \end{matrix}),

其中每一部分t_i、s_i以及w_i被拆分成内部子部分和接口子部分：

t_{i} = (\begin{matrix} t_{i}^{I} \\ t_{i}^{B} \end{matrix}) = (\begin{matrix} x_{i} \\ y_{i} \end{matrix}),

s_{i} = (\begin{matrix} s_{i}^{I} \\ s_{i}^{B} \end{matrix}) = (\begin{matrix} r_{i} \\ q_{i} \end{matrix}),

w_{i} = (\begin{matrix} w_{i}^{I} \\ w_{i}^{B} \end{matrix}) = (\begin{matrix} u_{i} \\ v_{i} \end{matrix}) .

应该记得，第i块的因式分解具有如下结构：

然后，根据该说明性实施例，预处理器求解过程的算法可以写为如下：

假设来自PMLIL U.Construct()：ML structure of L，U factors，last_level

因此，根据该说明性实施例，求解过程包括：

1、在多级方法的最后一级中的串行ILU解法；

2、将该内部分割器应用于向量v、q，其暗示在并行处理中使用的处理器之间的一些数据交换；

3、应用逆内部分割器以恢复使用在并行处理中使用的处理器之间的向量v的初始分布。

针对P＝4，L＝2的例子。为了进一步说明，考虑一个例子的部分数等于4，级数等于2，且在最后一级进行ILU型因式分解。根据该说明性实施例，如下所述这样执行该构造过程。

1)第一级。在外部初始分割成4个部分之后，该系统具有如下形式：

(\begin{matrix} A_{11}^{1} & A_{12}^{1} & A_{13}^{1} & A_{14}^{1} \\ A_{21}^{1} & A_{22}^{1} & A_{23}^{1} & A_{24}^{1} \\ A_{31}^{1} & A_{32}^{1} & A_{33}^{1} & A_{34}^{1} \\ A_{41}^{1} & A_{42}^{1} & A_{43}^{1} & A_{44}^{1} \end{matrix}) (\begin{matrix} x_{1}^{1} \\ x_{2}^{1} \\ x_{3}^{1} \\ x_{4}^{1} \end{matrix}) = (\begin{matrix} b_{1}^{1} \\ b_{2}^{1} \\ b_{3}^{1} \\ b_{4}^{1} \end{matrix}) .

其中，上标1是级数。

应用局部重排序并转换到加边对角块(BBD)的形式导致以下结构：

(\begin{matrix} A_{1}^{1} & F_{1}^{1} \\ A_{2}^{1} & F_{2}^{1} \\ A_{3}^{1} & F_{3}^{1} \\ A_{4}^{1} & F_{4}^{1} \\ C_{1}^{1} & B_{1}^{1} & A_{12}^{1} & A_{13}^{1} & A_{14}^{1} \\ C_{2}^{1} & A_{21}^{1} & B_{2}^{1} & A_{23}^{1} & A_{24}^{1} \\ C_{3}^{1} & A_{31}^{1} & A_{32}^{1} & B_{3}^{1} & A_{34}^{1} \\ C_{4}^{1} & A_{41}^{1} & A_{42}^{1} & A_{43}^{1} & B_{4}^{1} \end{matrix})

对对角块应用该并行截断因式分解导出下面的LU因式分解：

(\begin{matrix} L_{1}^{1} \\ L_{2}^{1} \\ L_{3}^{1} \\ L_{4}^{1} \\ L_{1}^{C 1} & I_{1} \\ L_{2}^{C 1} & I_{2}^{1} \\ L_{3}^{C 1} & I_{3}^{1} \\ L_{4}^{C 1} & I_{4}^{1} \end{matrix}) \times (\begin{matrix} U_{1}^{1} & U_{1}^{F 1} \\ U_{2}^{1} & U_{2}^{F 1} \\ U_{3}^{1} & U_{3}^{F 1} \\ U_{4}^{1} & U_{4}^{F 1} \\ S_{1}^{1} & A_{12}^{1} & A_{13}^{1} & A_{14}^{1} \\ A_{21}^{1} & S_{2}^{1} & A_{23}^{1} & A_{24}^{1} \\ A_{31}^{1} & A_{32}^{1} & S_{3}^{1} & A_{34}^{1} \\ A_{41}^{1} & A_{42}^{1} & A_{43}^{1} & S_{4}^{1} \end{matrix})

假定该接口矩阵的大小被确定为足够大(如在图2的框207中)，过程进行到第二级，这将在下面论述。

2)第二级。首先，第一级接口矩阵被再次进行分割，如下所示：

S^{1} = (\begin{matrix} S_{1}^{1} & A_{12}^{1} & A_{13}^{1} & A_{14}^{1} \\ A_{21}^{1} & S_{2}^{1} & A_{23}^{1} & A_{24}^{1} \\ A_{31}^{1} & A_{32}^{1} & S_{3}^{1} & A_{34}^{1} \\ A_{41}^{1} & A_{42}^{1} & A_{43}^{1} & S_{4}^{1} \end{matrix})

通过分割器Partitioner_IM。应当注意，使用串行或并行分割，整个矩阵被再次分割，包括舒尔补。由于这个原因，在该说明性实施例中实现并行分割器，其中该并行分割器能够为接口矩阵的每个块行带构造高质量的并行分割。

在重新分割之后，得到下面的矩阵：

A^{2} = (\begin{matrix} A_{11}^{2} & A_{12}^{2} & A_{13}^{2} & A_{14}^{2} \\ A_{21}^{2} & A_{22}^{2} & A_{23}^{2} & A_{24}^{2} \\ A_{31}^{2} & A_{32}^{2} & A_{33}^{2} & A_{34}^{2} \\ A_{41}^{2} & A_{42}^{2} & A_{43}^{2} & A_{44}^{2} \end{matrix}),

并应用上面描述的过程来构造

因此得到第二级接口矩阵：

S^{2} = (\begin{matrix} S_{1}^{2} & A_{12}^{2} & A_{13}^{2} & A_{14}^{2} \\ A_{21}^{2} & S_{2}^{2} & A_{23}^{2} & A_{24}^{2} \\ A_{31}^{2} & A_{32}^{2} & S_{3}^{2} & A_{34}^{2} \\ A_{41}^{2} & A_{42}^{2} & A_{43}^{2} & S_{4}^{2} \end{matrix})

由于达到允许的最大级数，对上述的接口矩阵S²执行最后一级因式分解：

在该实施例中，允许的最大级数是该算法(见算法2)的一个参数。而且，在上面例子中，最大级数设为2。

现在，初始化步骤已经完成，该迭代求解器继续进行求解过程。

第一级的L solve矩阵可以写为如下：

(\begin{matrix} L_{1}^{1} \\ L_{2}^{1} \\ L_{3}^{1} \\ L_{4}^{1} \\ L_{1}^{C 1} & I_{1} \\ L_{2}^{C 1} & I_{2}^{1} \\ L_{3}^{C 1} & I_{3}^{1} \\ L_{4}^{C 1} & I_{4}^{1} \end{matrix}) (\begin{matrix} u_{1}^{1} \\ u_{2}^{1} \\ u_{3}^{1} \\ u_{4}^{1} \\ v_{1}^{1} \\ v_{2}^{1} \\ v_{3}^{1} \\ v_{4}^{1} \end{matrix}) = (\begin{matrix} r_{1}^{1} \\ r_{2}^{1} \\ r_{3}^{1} \\ r_{4}^{1} \\ q_{1}^{1} \\ q_{2}^{1} \\ q_{3}^{1} \\ q_{4}^{1} \end{matrix})

并行求解

并代入就得到第一级接口矩阵S¹的右手边向量为：然后，该迭代求解器以不同次序排列并重新分布向量v，使s²＝Part_IM(v¹)，并重复上述的过程：

和

以在第二级上执行L solve。

然后，该迭代求解器在最后一级执行完全求解：L_LLU_LLy²＝v²。之后，该说明性实施例的迭代求解器以从第二级开始的向后顺序递归执行U solve。

现在进一步详细考虑该说明性实施例在第二级U solve。由最后一级预处理器求解获得的

被用于并行修改右手边向量，然后求解该系统

(\begin{matrix} U_{1}^{2} \\ U_{2}^{2} \\ U_{3}^{2} \\ U_{4}^{2} \end{matrix}) (\begin{matrix} t_{1}^{2} \\ t_{2}^{2} \\ t_{3}^{2} \\ t_{4}^{2} \end{matrix}) = (\begin{matrix} u_{1}^{2} - U_{1}^{F 2} y_{1}^{2} \\ u_{2}^{2} - U_{2}^{F 2} y_{2}^{2} \\ u_{3}^{2} - U_{3}^{F 2} y_{3}^{2} \\ u_{4}^{2} - U_{4}^{F 2} y_{4}^{2} \end{matrix})

应用逆置换和重新分布y¹＝invPart_IM(t²)，该迭代求解器可以应用上述的算法在第一级上执行U solve：

(\begin{matrix} U_{1}^{1} \\ U_{2}^{1} \\ U_{3}^{1} \\ U_{4}^{1} \end{matrix}) (\begin{matrix} t_{1}^{1} \\ t_{2}^{1} \\ t_{3}^{1} \\ t_{4}^{1} \end{matrix}) = (\begin{matrix} u_{1}^{1} - U_{1}^{F 1} y_{1}^{1} \\ u_{2}^{1} - U_{2}^{F 1} y_{2}^{1} \\ u_{3}^{1} - U_{3}^{F 1} y_{3}^{1} \\ u_{4}^{1} - U_{4}^{F 1} y_{4}^{1} \end{matrix})

因此，上述的说明性实施例使用针对大型稀疏线性系统的并行求解的方法，该方法实现高度并行化的执行因式分解方案。最优的变型允许非常少的串行工作，其可能还不到整个工作的1％，但是允许获得并行预处理器，其具有与对应的串行预处理器几乎相同的质量(在收敛要求的迭代求解器的迭代次数方面)。此外，应用纯并行局部重排序和缩放可显著地提高预处理器的质量。

实施例或其部分可在程序或代码段中具体化，所述程序或代码段在基于处理器的系统(如计算机系统)可操作用以执行在这里描述的针对并行计算迭代求解器的功能和操作。组成不同的实施例的程序或代码段可以存储在计算机可读介质上，其可包括临时或永久存储这类代码的任何合适介质。计算机可读介质的例子包括物理计算机可读介质，如电子存储器电路、半导体存储器器件、随机存取存储器(RAM)、只读存储器(ROM)、可擦除只读存储器(EROM)、闪存、磁存储装置(例如软盘)、光存储装置(例如光盘(CD)、数字多功能盘(DVD))、硬盘等等。

图4示出根据本发明一个实施例的示例性计算机系统400，执行上述并行计算迭代求解器的处理操作的软件可在其上实现。中央处理单元(CPU)401耦合到系统总线402。虽然示出单独的CPU 401，但是应该认识到计算机系统400优选包括多个处理单元(如CPU 401)，其用在上述并行计算中。CPU 401可以是任何通用的CPU。本发明不受限于CPU 401的架构(或示例性系统400的其他组件)，只要CPU401(和系统400的其他组件)支持这里描述的本发明的操作。根据上述的实施例，CPU 401可运行不同逻辑指令。例如CPU 401可运行机器级指令，用于执行根据上面结合图2和图3所述的并行计算迭代求解器的实施例的示例性操作流的处理。

计算机系统400也优选包括随机存取存储器(RAM)403，其可以是SRAM、DRAM、SDRAM等。计算机系统400优选包括只读存储器(ROM)404，其可以是PROM、EPROM、EEPROM等。RAM 403和ROM 404保存用户数据、系统数据及程序，这是本领域公知的。

计算机系统400也优选包括输入/输出(I/O)适配器405、通信适配器411、用户接口适配器408以及显示适配器409。I/O适配器405、用户接口适配器408和/或通信适配器411在某些实施例中可使得用户能够与计算机系统400交互，以便输入信息。

I/O适配器405优选连接到(多个)存储装置406上，如一个或更多硬盘驱动器、光盘(CD)驱动器、软盘驱动器、磁带驱动器等连接到计算机系统400。存储装置可在RAM403不足以满足与存储用于本发明实施例的操作的数据关联的存储器要求时使用。计算机系统400的数据存储可用于存储这样的信息，如模型(例如图2-3的模型223)、由并行计算迭代求解器计算的中间或最终结果和/或根据本发明实施例的使用或生成的其他数据。通信适配器411优选适于耦合计算机系统400到网络412上，这使得信息能够经网络412(如因特网或其他广域网、局域网、公共或专用交换电话网络、无线网络、前述网络的任何组合)输入系统400和/或从其中输出。用户接口适配器408耦合输入装置，如键盘413、指示装置407、麦克风414和/或输出装置，如扬声器415到计算机系统400。显示适配器409由CPU驱动从而控制显示装置410上的显示，从而例如显示正在分析的和模型相关的信息，如显示生成的地下含碳氢化合物储层中流体流随时间的3D表示。

应当理解，本发明并不限于系统400的架构。例如，可以使用任何合适的基于处理器的装置，用于实施本发明所有或部分实施例，包括不限于个人计算机、便携式计算机、计算机工作站、服务器、和/或其他多处理器计算装置。此外，实施例可在专用集成电路上(ASIC)或超大规模集成电路(VLSI)上实施。事实上，本领域普通技术人员可使用任何数目的能够执行根据实施例的逻辑操作的合适结构。

尽管已经详细说明了本发明及其优势，但应当理解可以做出各种变化、替代和改变而不偏离权利要求限定的本发明的精神和范畴。而且，本申请的范畴不限于说明书中说明的过程、机器、制造、组成、手段、方法和步骤的特定实施例。本领域普通技术人员很容易从本发明的公开内容理解，可根据本发明利用执行与这里所述相应实施例基本相同的功能或实现基本相同结果的现有或以后开发的过程、机器、制造、组成、手段、方法或步骤。因此，权利要求在其范畴中包括这样的过程、机器、制造、组成、手段、方法或步骤。

Claims

1.一种方法，其包括：

(a)使用多级分割将初始矩阵分割成多个子矩阵；

(b)对对角块并行执行截断因式分解，为所述多个子矩阵的每一个的接口部分形成局部舒尔补，；

(c)由所述多个子矩阵的对角块上的局部舒尔补和非对角块上的所述多个子矩阵之间的连接形成全局接口矩阵；

(d)判定以下至少一个：i)该全局接口矩阵是否足够小，以满足预定义的大小阈值，且ii)是否达到最后的允许级；以及

(e)当判定该全局接口矩阵不是小得足以满足所述预定义的大小阈值或者达到该最后的允许级，则使用多级分割将所述全局接口矩阵分割成多个第二子矩阵和为所述多个第二子矩阵重复操作(b)到(e)。

2.根据权利要求1所述的方法，进一步包括：

(f)当判定该全局接口矩阵足够小以满足所述预定义的大小阈值时，因式分解所述全局接口矩阵。

3.根据权利要求1所述的方法，进一步包括：

对所述多个子矩阵的每一个进行重排序。

4.根据权利要求3所述的方法，其中所述重排序在操作(a)之后且操作(b)之前执行。

5.根据权利要求3所述的方法，其中所述重排序包括：

首先重排序所述多个子矩阵的每个的内部行；和

其次重排序所述多个子矩阵的接口行。

6.根据权利要求1所述的方法，进一步包括：

在所述多个子矩阵上执行局部缩放算法，以改进所述子矩阵的数值性质。

7.根据权利要求1所述的方法，其中所述形成所述全局接口矩阵包括：

隐式形成所述全局接口矩阵。

8.根据权利要求7所述的方法，其中所述隐式形成所述全局接口矩阵包括：

通过多个处理单元中的每一个，存储所述接口矩阵的相应部分。

9.根据权利要求1所述的方法，其中所述并行执行包括：

通过多个并行处理单元执行所述操作(b)。

10.根据权利要求1所述的方法，其中所述使用多级分割将所述全局接口矩阵分割成所述多个第二子矩阵包括：

使用所述接口矩阵的多级分割以避免在主处理单元显式形成所述接口矩阵。

11.根据权利要求1所述的方法，进一步包括：

对最后一级接口矩阵的处理。

12.根据权利要求11所述的方法，其中所述对最后一级接口矩阵的处理包括：

在所述最后一级，通过应用预定义的预处理器对相应的接口矩阵进行因式分解。

13.根据权利要求12所述的方法，其中所述相应的接口矩阵串行进行因式分解。

14.根据权利要求12所述的方法，其中所述相应的接口矩阵并行进行因式分解。

15.根据权利要求11所述的方法，其中所述对最后一级接口矩阵的处理包括：

在所述最后一级通过应用串行高质量ILU因式分解对相应的接口矩阵进行因式分解。

16.根据权利要求11所述的方法，其中所述对最后一级接口矩阵的处理包括：

在所述最后一级通过应用并行迭代松弛的Block-Jacoby预处理器对相应的接口矩阵进行因式分解，并且对角块进行高质量的ILU因式分解。

17.一种方法，包括：

(a)使用多级分割将初始矩阵分割成多个子矩阵；

(b)对对角块并行执行截断因式分解，并且形成所述多个子矩阵的每一个的接口部分的局部舒尔补；

(d)判定该全局接口矩阵是否足够小，以满足预定义的大小阈值；以及

(e)当判定该全局接口矩阵不是足够小以满足所述预定义的大小阈值时，使用多级分割将所述全局接口矩阵分割成多个第二子矩阵和为所述多个第二子矩阵重复操作(b)到(d)。

18.根据权利要求17所述的方法，进一步包括：

对所述多个子矩阵中的每个执行局部重排序。

19.根据权利要求18所述的方法，其中所述局部重排序是在步骤(a)之后且步骤(b)之前执行的。

20.一种方法，包括：

使用p-way多级分割，对初始矩阵应用不交迭的域分解以将所述初始矩阵分割成p个部分，因此形成多个子矩阵：

预定义允许的最大递归级数；

预定义最小大小阈值，其规定相对于所述初始矩阵大小的接口矩阵的允许的最小行数；

递归执行操作(a)到(d)：

(a)通过多个并行处理单元对所述多个子矩阵中的每一个并行执行：i)对角块的并行截断因式分解；和ii)为每个子矩阵的接口部分形成局部舒尔补；

(b)由对角块上的局部舒尔补和非对角块上子矩阵间的连接隐式形成全局接口矩阵；

(c)判定是否达到预定义的允许的最大递归级数或该全局接口矩阵的大小是否小于预定义的最小大小阈值；

(d)当在操作(c)中判定没有达到预定义的允许的最大递归级数并且该全局接口矩阵的大小不小于预定义的最小大小阈值时，使用多级分割将该全局接口矩阵分割成另外的多个子矩阵和针对所述另外的多个子矩阵重复操作(a)-(d)。

21.根据权利要求20所述的方法，进一步包括：

当在操作(c)中判定达到预定义的允许的最大递归级数或该全局接口矩阵的大小小于预定义的最小大小阈值，则结束该递归处理。

22.根据权利要求21所述的方法，进一步包括：

当在操作(c)中判定达到预定义的允许的最大递归级数或该全局接口矩阵的大小小于预定义的最小大小阈值，则对所述全局接口矩阵进行因式分解。

23.根据权利要求20所述的方法，进一步包括：

对所述多个子矩阵的每一个执行局部重排序。

24.根据权利要求23所述的方法，其中所述重排序包括：

首先对所述多个子矩阵的每一个的内部行进行重排序；以及

其次对所述多个子矩阵的接口行进行重排序。

25.根据权利要求20所述的方法，其中所述隐式形成包括：

通过所述多个并行处理单元的每一个存储所述接口矩阵的相应部分。