用于在并行处理器上再因式分解方形矩阵的系统和方法

摘要

一种用于在并行处理器上再因式分解方形输入矩阵的系统和方法。在一个实施例中，该系统包括：（1）矩阵生成器，其可操作为在下和上三角矩阵组合的归零稀疏模式中嵌入一个置换形式的输入矩阵，以便产生中间矩阵，其中所述下和上三角矩阵源于对在先矩阵的在先LU因式分解，所述在先矩阵与所述输入矩阵具有相同的稀疏模式，并且通过重排序处理而将填充和枢纽策略最小化，以及（2）与矩阵生成器关联的再因式分解器，其可操作为使用并行线程来将用零填充的不完全LU因式分解应用于所述中间矩阵。

说明书

技术领域

本申请主要涉及并行处理器（也就是具有至少两个能够通过协作来执行 并行处理的处理器的计算机），尤其涉及的是用于在并行处理器上将方形矩阵 因式分解（refactorizing）成下和上三角矩阵的系统和方法。

背景技术

如果可以有效地对图形处理器（GPU）之类的并行处理器进行编程，那 么它们有可能会非常擅于处理数值算法，尤其是用于直接求解大型稀疏线性 系统的算法。

稀疏线性系统是具有系数系数矩阵的线性方程系统。在计算力学、地球 物理学、生物学、电路仿真的上下文以及计算科学和工程领域的其他众多的 上下文中都会出现此类系统。

用于求解稀疏线性系统的最常见的通用直接的技术是将其系数矩阵分解 成下三角矩阵L和上三角矩阵U的乘积，这种处理被称为“因式分解”。然 后，通过使用常规的前向和反向替换技术，可以使用L和U三角矩阵来求解 该线性系统，并且由此获取稀疏线性系统的解。

发明内容

这里的一个方面提供的是一种用于在并行处理器上再因式分解方形输入 矩阵的系统。在一个实施例中，该系统包括：（1）矩阵生成器，其可操作为 在下和上三角矩阵组合的归零稀疏模式中嵌入一个置换形式的输入矩阵，以 便产生中间矩阵，其中所述下和上三角矩阵源于对在先矩阵的LU因式分解， 所述在先矩阵与所述输入矩阵具有相同的稀疏模式，并且通过重排序处理而 将填充和枢纽策略最小化，以及（2）与矩阵生成器关联的再因式分解器 （refactorizer），其可操作为使用并行线程来将用零填充的不完全LU因式分 解（ILU0）应用于所述中间矩阵。

这里的另一个方面提供了一种在并行处理器上再因式分解方形矩阵的方 法。在一个实施例中，该方法包括：（1）在下和上三角矩阵组合的归零稀疏 模式中嵌入一个置换形式的输入矩阵，以及通过重排序来将作为输入矩阵的 填充和枢纽策略最小化，从而产生中间矩阵，其中所述下和上三角矩阵源于 对具有相同稀疏模式的在先矩阵的LU因式分解，以及（2）使用并行线程来 将用零填充的不完全LU因式分解（ILU0）应用于所述中间矩阵。

这里的另一个方面提供了一种SIMD（单指令多数据）处理器。在一个 实施例中，SIMD处理器包括：（1）可操作为处理并行线程的通道，（2）管 线控制单元系统，其可操作为控制通道中的线程处理，以及（3）通过使用通 道和管线控制单元来再因式分解方形输入矩阵的系统，该系统具有：（3a）矩 阵生成器，其可操作为在下和上三角矩阵组合的归零稀疏模式中嵌入一个置 换形式的输入矩阵，并且通过重排序来将作为输入矩阵的填充和枢纽策略最 小化，从而产生中间矩阵，其中所述下和上三角矩阵源于对具有相同稀疏模 式的在先矩阵所进行的LU因式分解，以及（3b）与矩阵生成器关联的再因 式分解器（refactorizer），其可操作为使用并行线程来将用零填充不完全LU 因式分解（ILU0）的应用于所述中间矩阵。

附图说明

现在将结合附图来参考以下描述，其中：

图1是可操作为包含或者执行用于在并行处理器上将方形矩阵再因式分 解成下和上三角矩阵的系统或方法的SIMD处理器的框图。

图2是用于在并行处理器上将方形矩阵再因式分解成下和上三角矩阵的 系统的一个实施例的流程图。

图3是用于在并行处理器上将方形矩阵再因式分解成下和上三角矩阵的 方法的一个实施例的框图。

具体实施方式

如上所述，用于求解稀疏线性系统的最常见的通用方法是将其系数矩阵 分解因式分解成下和上三角矩阵L和U的乘积。在这里可以认识到，现有技 术并不适合利用GPU之类的并行处理器的架构。

相应地，在这里引入了用于将方形矩阵再因式分解成下和上三角矩阵的 系统和方法的不同实施例。通常，这里的不同实施例可用于加速用于求解如 下形式的稀疏线性系统集合的应用：

Ai xi=fi for i=1,…,k,

其中系数矩阵Ai，右侧的fi，并且解xi。

这里的系统和方法的某些实施例是当（1）系数矩阵Ai的稀疏模式，（2） 用于最小化填充的重排序处理以及（3）用在因数分解过程中的枢纽策略在所 有线性系统中全都保持相同的时候适用的。在这种情况下，所产生的用于每 个线性系统的下（Li）和上（Ui）三角因数的稀疏模式同样会保持相同。在 使用众所周知的重点应用于集成电路的仿真程序（SPICE）来仿真集成电路 的过程中，这些状况将会频繁出现。

在这些特定的实施例中，LU因式分解仅仅需要在第一次（i=1）执行。 随后（i=2，……，k），所要执行的仅仅是LU因式分解。由于将LU因式分 解应用于Li和Ui因数的稀疏模式不会产生附加填充，因此可以静态地分配 所述再因式分解所需要的存储器。由此，LU再因式分解通常要远远快于LU 因式分解。在一些实施例中，在现代的SIMD处理器中可以在数十秒或数分 钟的时间里执行再因式分解，而因数分解则有可能需要花费好多个小时。

由此，在这里可以认识到，对于i=2，……，k，在源自第一（i=1）因式 分解的归零下和上三角因数的稀疏模式中可以嵌入系数矩阵Ai，即

Mi=L1(z)+U1(z)+Ai，

其中L1(z)=L1以及U1(z)=U1是用零填充的。由于将LU因式分解应 用于矩阵Li和Ui的处理不会产生附加填充，因此，在这里可以认识到，在 这个新产生的中间矩阵Mi中可以应用ILU0，以便产生系数矩阵Ai的LU再 因式分解。

在某些实施例中，通过嵌入置换矩阵而不是系数矩阵，可以说明用于将 填充和枢纽减至最少的重排序处理，由此：

Mi=L1(z)+U1(z)+PT*Ai*Q,

其中PT和Q是与用于将第一（i=1）LU因式分解中使用的填充和枢纽 减至最少的重排序处理相对应的置换矩阵。

由此，关于系数矩阵Ai的LU再因式分解的问题被改写成了对于中间矩 阵Mi的ILU0的再因式分解，其中i=2，……，k。可以认识到的是，后者可 以在诸如GPU之类的并行处理器上得到有效执行。在更详细地描述新颖的系 统和方法之前，现在将对包含GPU的典型计算系统进行描述。

图1是可操作为包含或执行用于在并行处理器上将方形矩阵再因式分解 成下和上三角矩阵的系统或方法的SIMD处理器100的框图。SIMD处理器 100包括被组织成线程组104或“包（wrap）”的多线程处理器或核心106。 SIMD处理器100包括J个线程组104-1到104-J，其中每一个线程组都具有 K个核心106-1到106-K。在某些实施例中，线程组104-1到104-J可以进一 步被组织成一个或多个线程块102。在一个具体的实施例中，每一个线程组 104具有32个核心106。其他实施例可以在线程组中只包含4个核心以及包 含多达数万个核心。某些实施例将核心106组织成单个线程组104，而其他 实施例则可以具有数百乃至数千个线程组104。关于SIMD处理器100的替 换实施例可以仅仅将核心106组织成线程组104，由此忽略线程块组织等级。

SIMD处理器100还包括管线控制单元108，共享存储器110以及与线程 组104-1到104-J相关联的逻辑存储器阵列112-1到112-J。线程运行控制单 元108通过数据总线114来将任务分发给不同的线程组104-1到104-J。线程 组内部的核心106以相互并行的方式运行。线程组104-1到104-J通过存储器 总线116来与共享存储器110进行通信。所述线程组104-1到104-J分别通过 本地总线118-1到118-J来与本地存储器112-1到112-J进行通信。作为示例， 线程组104-J是通过在本地总线118-J上进行通信来使用本地存储器112-J的。 关于SIMD处理器的某些实施例将共享存储器110的共享部分分配给了每一 个线程块102，并且允许线程块102内部的所有线程组104访问共享存储器 110的共享部分。某些实施例包括只使用本地存储器112的线程组104。很多 其他的实施例包括用于平衡本地存储器112和共享存储器110的使用的线程 组104。

在描述了包含GPU的典型计算系统之后，现在将更详细地描述关于这里 的新颖系统和方法的不同实施例。所述系统和方法的不同实施例使用了不同 的新颖技术。

图2是用于在并行处理器上将方形矩阵再因式分解成下和上三角矩阵的 系统的一个实施例的框图。在很多实施例中，因式分解器210并不是系统的 一部分，所述因式分解器可操作为接收作为输入矩阵的方形矩阵A1，并且通 过对该矩阵执行LU因式分解来产生下三角矩阵L1和上三角矩阵U1，此外 它还执行用于将因式分解中使用的填充和枢纽置换P和Q减至最少的重排序 处理。

矩阵生成器220与因式分解器210相关联，并且可操作为通过在源于第 一LU因式分解的下和上三角矩阵组合的归零稀疏模式L1(z)+U1(z)中嵌入置 换输入矩阵PT*Ai*Q来产生中间矩阵Mi，其中i=2,...,k。再因式分解器230 与矩阵生成器220相关联，其可操作为使用并行线程来对所产生的中间矩阵 Mi应用ILU0，其中i=2,...,k。

在这里应该认识到，行更新本质上是一个高斯消去步骤，其中参考（枢 纽）行被用一个乘数缩放并添加至当前行中，以便创建低于主对角线的零值。 这些乘数是代替该更新创建的零值而被计算和保存的。在关于新颖系统和方 法的一些实施例中，为了探测可用的并行性，列被分组成了等级而不是行。 在这里应该认识到，将行分组成等级的常规方法不但更为直观，而且还会使 发现数据依存性的处理更为简单。然而，在这里应该认识到，将行分组成等 级的处理极大阻碍了处理速度。尤其应该认识到的是，与其他的行相比，接 近矩阵底部的行往往具有多出很多的非零元素，由此会使这些行的更新相对 较慢。相比之下，如果将列分组成行，那么往往允许我们通过将与底部的行 相关联的计算均匀分布在多个列上来实现更好的工作负载平衡。通过将列分 组成等级，可以显著提升处理速度，并且有可能将其提升好几个数量级。

通过将列分组成等级，还可以产生另一个潜在的优点。矩阵的右下角可 以封装到密集存储中，并且可以使用密集LU因式分解（没有枢纽）来对其 进行处理。在一个实施例中，只有在右下角的每一行都包含了足以证明采用 密集格式的额外计算的正当性的元素的情况下，所述封装处理才会发生。

在系统和方法的其他实施例中，矩阵是用合并的压缩稀疏行/压缩稀疏列 （CSR-CSC）的格式保持的。所合并的CSR-CSC格式同时包含了标准的CSR 和CSC格式。其中，CSR格式是常规格式，并且包含了矩阵值阵列。CSC 格式也是常规格式，但其包含的是指向CSR格式中包含的值的指针，而不是 实际矩阵值。合并的CSR-CSC格式允许在ILU0期间的适当位置对采用CSR 格式的单个数值阵列进行更新。而CSC格式指针则不需要更新。

在系统和方法的其他实施例中，由于矩阵Mi中的每一行的元素的数量 在早期和晚期会发生很大变化，因此，分析阶段会为每一个等级计算单独的 网格和块启动参数。

在系统和方法的其他实施例中，所述分析可以在每一级调度两个内核启 动，而不是使用单个内核启动来处理一个等级（在每一行的x维度具有单个 线程块）。如果启动两个内核，那么第一内核可以更新乘数（在每一行的x维 度上具有单个线程块），并且第二内核可以更新行中的剩余元素（每一行的x 维度中具有多个线程块）。

在某些实施例中，在LU再因式分解的行更新期间会在当前行中搜索参 考（枢纽）行的元素，而不是在参考（枢纽）行中搜索当前行的元素。以这 种方式进行的搜索将会产生两个潜在的优点。第一，在给出了Mi的构造的情 况下，如果参考（枢纽）行中的元素的数量为“n”，并且当前行中的元素的 数量为“m”，那么m≥n。由此，前一种和后一种方法分别包括O(m*log(n)) 和O(n*log(m))个步骤。由此，后一种方法包含较少的步骤，并且计算成本相 对较低。此外，通过构造Mi，可以知道当前行中始终存在参考（枢纽）行中 的元素，由此可以将线程分散度减至最小。

图3是用于在并行处理器上将方形矩阵再因式分解成下和上三角矩阵的 方法的一个实施例的流程图。该方法始于开始步骤310。在步骤320，对输入 矩阵A1执行LU因式分解，以便产生下三角矩阵L1和上三角矩阵U1，以及 通过对其执行重排序处理来将因式分解中使用的填充和枢纽置换P和Q最小 化。在步骤330，对于i=2，……，k，通过在源自在先LU因式分解的下和上 三角矩阵组合的归零稀疏模式L1(z)+U1(z)中嵌入置换形式的输入矩阵 PT*Ai*Q来产生一个中间矩阵Mi。在步骤340，通过使用并行线程来将ILU0 应用于所产生的中间矩阵Mi。在结束步骤350，该方法结束。

本申请的领域的技术人员将会认识到，针对所描述的实施例的其他更进 一步的补充、删除、替换和修改都是可行的。