用于在并行处理器上再因式分解方形矩阵的系统和方法
【专利摘要】一种用于在并行处理器上再因式分解方形输入矩阵的系统和方法。在一个实施例中，该系统包括：（1）矩阵生成器，其可操作为在下和上三角矩阵组合的归零稀疏模式中嵌入一个置换形式的输入矩阵，以便产生中间矩阵，其中所述下和上三角矩阵源于对在先矩阵的在先LU因式分解，所述在先矩阵与所述输入矩阵具有相同的稀疏模式，并且通过重排序处理而将填充和枢纽策略最小化，以及（2）与矩阵生成器关联的再因式分解器，其可操作为使用并行线程来将用零填充的不完全LU因式分解应用于所述中间矩阵。
【专利说明】用于在并行处理器上再因式分解方形矩阵的系统和方法
【技术领域】
[0001]本申请主要涉及并行处理器(也就是具有至少两个能够通过协作来执行并行处理的处理器的计算机)，尤其涉及的是用于在并行处理器上将方形矩阵因式分解(refactorizing)成下和上三角矩阵的系统和方法。
【背景技术】
[0002]如果可以有效地对图形处理器(GPU)之类的并行处理器进行编程，那么它们有可能会非常擅于处理数值算法，尤其是用于直接求解大型稀疏线性系统的算法。
[0003]稀疏线性系统是具有系数系数矩阵的线性方程系统。在计算力学、地球物理学、生物学、电路仿真的上下文以及计算科学和工程领域的其他众多的上下文中都会出现此类系统。
[0004]用于求解稀疏线性系统的最常见的通用直接的技术是将其系数矩阵分解成下三角矩阵L和上三角矩阵U的乘积，这种处理被称为“因式分解”。然后，通过使用常规的前向和反向替换技术，可以使用L和U三角矩阵来求解该线性系统，并且由此获取稀疏线性系统的解。
【发明内容】
[0005]这里的一个方面提供的是一种用于在并行处理器上再因式分解方形输入矩阵的系统。在一个实施例中，该系统包括:(I)矩阵生成器，其可操作为在下和上三角矩阵组合的归零稀疏模式中嵌入一个置换形式的输入矩阵，以便产生中间矩阵，其中所述下和上三角矩阵源于对在先矩阵的LU因式分解，所述在先矩阵与所述输入矩阵具有相同的稀疏模式，并且通过重排序处理而将填充和枢纽策略最小化，以及(2)与矩阵生成器关联的再因式分解器(refactorizer)，其可操作为使用并行线程来将用零填充的不完全LU因式分解(ILUO)应用于所述中间矩阵。
[0006]这里的另一个方面提供了一种在并行处理器上再因式分解方形矩阵的方法。在一个实施例中，该方法包括:(1)在下和上三角矩阵组合的归零稀疏模式中嵌入一个置换形式的输入矩阵，以及通过重排序来将作为输入矩阵的填充和枢纽策略最小化，从而产生中间矩阵，其中所述下和上三角矩阵源于对具有相同稀疏模式的在先矩阵的LU因式分解，以及(2)使用并行线程来将用零填充的不完全LU因式分解(ILUO)应用于所述中间矩阵。
[0007]这里的另一个方面提供了一种SIMD(单指令多数据)处理器。在一个实施例中，SIMD处理器包括:(I)可操作为处理并行线程的通道，(2)管线控制单元系统，其可操作为控制通道中的线程处理，以及(3)通过使用通道和管线控制单元来再因式分解方形输入矩阵的系统，该系统具有:(3a)矩阵生成器，其可操作为在下和上三角矩阵组合的归零稀疏模式中嵌入一个置换形式的输入矩阵，并且通过重排序来将作为输入矩阵的填充和枢纽策略最小化，从而产生中间矩阵，其中所述下和上三角矩阵源于对具有相同稀疏模式的在先矩阵所进行的LU因式分解，以及(3b)与矩阵生成器关联的再因式分解器(refactorizer)，其可操作为使用并行线程来将用零填充不完全LU因式分解(ILUO)的应用于所述中间矩阵。
【专利附图】
【附图说明】
[0008]现在将结合附图来参考以下描述，其中:
[0009]图1是可操作为包含或者执行用于在并行处理器上将方形矩阵再因式分解成下和上三角矩阵的系统或方法的SMD处理器的框图。
[0010]图2是用于在并行处理器上将方形矩阵再因式分解成下和上三角矩阵的系统的一个实施例的流程图。
[0011]图3是用于在并行处理器上将方形矩阵再因式分解成下和上三角矩阵的方法的一个实施例的框图。
【具体实施方式】
[0012]如上所述，用于求解稀疏线性系统的最常见的通用方法是将其系数矩阵分解因式分解成下和上三角矩阵L和U的乘积。在这里可以认识到，现有技术并不适合利用GPU之类的并行处理器的架构。[0013]相应地，在这里引入了用于将方形矩阵再因式分解成下和上三角矩阵的系统和方法的不同实施例。通常，这里的不同实施例可用于加速用于求解如下形式的稀疏线性系统集合的应用:
[0014]Ai xi=fi for i=l，...，k,
[0015]其中系数矩阵Ai，右侧的f i，并且解Xi。
[0016]这里的系统和方法的某些实施例是当(I)系数矩阵Ai的稀疏模式，(2)用于最小化填充的重排序处理以及(3)用在因数分解过程中的枢纽策略在所有线性系统中全都保持相同的时候适用的。在这种情况下，所产生的用于每个线性系统的下(Li)和上(Ui)三角因数的稀疏模式同样会保持相同。在使用众所周知的重点应用于集成电路的仿真程序(SPICE)来仿真集成电路的过程中，这些状况将会频繁出现。
[0017]在这些特定的实施例中，LU因式分解仅仅需要在第一次(i=l)执行。随后(i=2,……，k)，所要执行的仅仅是LU因式分解。由于将LU因式分解应用于Li和Ui因数的稀疏模式不会产生附加填充，因此可以静态地分配所述再因式分解所需要的存储器。由此，LU再因式分解通常要远远快于LU因式分解。在一些实施例中，在现代的SMD处理器中可以在数十秒或数分钟的时间里执行再因式分解，而因数分解则有可能需要花费好多个小时。
[0018]由此，在这里可以认识到，对于i=2，……，k，在源自第一(i=l)因式分解的归零下和上三角因数的稀疏模式中可以嵌入系数矩阵Ai，SP
[0019]Mi=Ll (z)+Ul (z)+Ai,
[0020]其中LI (Z)=Ll以及Ul (Z)=Ul是用零填充的。由于将LU因式分解应用于矩阵Li和Ui的处理不会产生附加填充，因此，在这里可以认识到，在这个新产生的中间矩阵Mi中可以应用ILU0，以便产生系数矩阵Ai的LU再因式分解。
[0021]在某些实施例中，通过嵌入置换矩阵而不是系数矩阵，可以说明用于将填充和枢纽减至最少的重排序处理，由此:[0022]Mi=Ll (z) +Ul (z) +PT*Ai*Q,
[0023]其中PT和Q是与用于将第一(i=l) LU因式分解中使用的填充和枢纽减至最少的重排序处理相对应的置换矩阵。
[0024]由此，关于系数矩阵Ai的LU再因式分解的问题被改写成了对于中间矩阵Mi的ILUO的再因式分解，其中i=2，……，k。可以认识到的是，后者可以在诸如GPU之类的并行处理器上得到有效执行。在更详细地描述新颖的系统和方法之前，现在将对包含GPU的典型计算系统进行描述。
[0025]图1是可操作为包含或执行用于在并行处理器上将方形矩阵再因式分解成下和上三角矩阵的系统或方法的SMD处理器100的框图。SMD处理器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。关于SMD处理器100的替换实施例可以仅仅将核心106组织成线程组104，由此忽略线程块组织等级。
[0026]SMD处理器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的。关于SMD处理器的某些实施例将共享存储器110的共享部分分配给了每一个线程块102，并且允许线程块102内部的所有线程组104访问共享存储器110的共享部分。某些实施例包括只使用本地存储器112的线程组104。很多其他的实施例包括用于平衡本地存储器112和共享存储器110的使用的线程组104。
[0027]在描述了包含GPU的典型计算系统之后，现在将更详细地描述关于这里的新颖系统和方法的不同实施例。所述系统和方法的不同实施例使用了不同的新颖技术。
[0028]图2是用于在并行处理器上将方形矩阵再因式分解成下和上三角矩阵的系统的一个实施例的框图。在很多实施例中，因式分解器210并不是系统的一部分，所述因式分解器可操作为接收作为输入矩阵的方形矩阵Al，并且通过对该矩阵执行LU因式分解来产生下三角矩阵LI和上三角矩阵Ul，此外它还执行用于将因式分解中使用的填充和枢纽置换P和Q减至最少的重排序处理。
[0029]矩阵生成器220与因式分解器210相关联，并且可操作为通过在源于第一 LU因式分解的下和上三角矩阵组合的归零稀疏模式LI (z)+Ul (z)中嵌入置换输入矩阵PT*Ai*Q来产生中间矩阵Mi，其中i=2，...，k。再因式分解器230与矩阵生成器220相关联，其可操作为使用并行线程来对所产生的中间矩阵Mi应用ILU0，其中i=2，...，k。
[0030]在这里应该认识到，行更新本质上是一个高斯消去步骤，其中参考(枢纽)行被用一个乘数缩放并添加至当前行中，以便创建低于主对角线的零值。这些乘数是代替该更新创建的零值而被计算和保存的。在关于新颖系统和方法的一些实施例中，为了探测可用的并行性，列被分组成了等级而不是行。在这里应该认识到，将行分组成等级的常规方法不但更为直观，而且还会使发现数据依存性的处理更为简单。然而，在这里应该认识到，将行分组成等级的处理极大阻碍了处理速度。尤其应该认识到的是，与其他的行相比，接近矩阵底部的行往往具有多出很多的非零元素，由此会使这些行的更新相对较慢。相比之下，如果将列分组成行，那么往往允许我们通过将与底部的行相关联的计算均匀分布在多个列上来实现更好的工作负载平衡。通过将列分组成等级，可以显著提升处理速度，并且有可能将其提升好几个数量级。
[0031]通过将列分组成等级，还可以产生另一个潜在的优点。矩阵的右下角可以封装到密集存储中，并且可以使用密集LU因式分解(没有枢纽)来对其进行处理。在一个实施例中，只有在右下角的每一行都包含了足以证明采用密集格式的额外计算的正当性的元素的情况下，所述封装处理才会发生。
[0032]在系统和方法的其他实施例中，矩阵是用合并的压缩稀疏行/压缩稀疏列(CSR-CSC)的格式保持的。所合并的CSR-CSC格式同时包含了标准的CSR和CSC格式。其中，CSR格式是常规格式，并且包含了矩阵值阵列。CSC格式也是常规格式，但其包含的是指向CSR格式中包含的值的指针，而不是实际矩阵值。合并的CSR-CSC格式允许在ILUO期间的适当位置对采用CSR格式的单个数值阵列进行更新。而CSC格式指针则不需要更新。
[0033]在系统和方法的其他实施例中，由于矩阵Mi中的每一行的元素的数量在早期和晚期会发生很大变化，因此，分析阶段会为每一个等级计算单独的网格和块启动参数。
[0034]在系统和方法的其他实施例中，所述分析可以在每一级调度两个内核启动，而不是使用单个内核启动来处理一个等级(在每一行的X维度具有单个线程块)。如果启动两个内核，那么第一内核可以更新乘数(在每一行的X维度上具有单个线程块)，并且第二内核可以更新行中的剩余元素(每一行的X维度中具有多个线程块)。
[0035]在某些实施例中，在LU再因式分解的行更新期间会在当前行中搜索参考(枢纽)行的元素，而不是在参考(枢纽)行中搜索当前行的元素。以这种方式进行的搜索将会产生两个潜在的优点。第一，在给出了 Mi的构造的情况下，如果参考(枢纽)行中的元素的数量为“n”，并且当前行中的元素的数量为“m”，那么m≥η。由此，前一种和后一种方法分别包括0(m*log(n))和0(n*log(m))个步骤。由此,后一种方法包含较少的步骤,并且计算成本相对较低。此外，通过构造Mi，可以知道当前行中始终存在参考(枢纽)行中的元素，由此可以将线程分散度减至最小。
[0036]图3是用于在并行处理器上将方形矩阵再因式分解成下和上三角矩阵的方法的一个实施例的流程图。该方法始于开始步骤310。在步骤320，对输入矩阵Al执行LU因式分解，以便产生下三角矩阵LI和上三角矩阵U1，以及通过对其执行重排序处理来将因式分解中使用的填充和枢纽置换P和Q最小化。在步骤330，对于i=2，……，k，通过在源自在先LU因式分解的下和上三角矩阵组合的归零稀疏模式LI (z) +Ul (z)中嵌入置换形式的输入矩阵PT*Ai*Q来产生一个中间矩阵Mi。在步骤340，通过使用并行线程来将ILUO应用于所产生的中间矩阵Mi。在结束步骤350，该方法结束。
[0037]本申请的领域的技术人员将会认识到，针对所描述的实施例的其他更进一步的补充、删除、替换和修改都是可行的。
【权利要求】
1.一种用于在并行处理器上再因式分解方形输入矩阵的系统，包括: 矩阵生成器，其可操作为在下和上三角矩阵的组合的归零稀疏模式中嵌入所述输入矩阵的置换形式，以便产生中间矩阵，其中所述下和上三角矩阵源自对在先矩阵的LU因式分解，所述在先矩阵与所述输入矩阵具有相同的稀疏模式，并且通过重排序处理而将填充和枢纽策略最小化；以及 与所述矩阵生成器关联的再因式分解器，其可操作为使用并行线程来将用零填充的不完全LU因式分解应用于所述中间矩阵。
2.如权利要求1所述的系统，其中所述中间矩阵生成器可操作为将所述置换形式嵌入所述稀疏模式。
3.如权利要求1所述的系统，其中所述再因式分解器还可操作为对所述输入矩阵执行符号分析，以便将代表其列的独立节点分组成代表所述节点之间的数据依存性的等级。
4.如权利要求1所述的系统，其中所述再因式分解器还可操作为使用合并的压缩稀疏行-压缩稀疏列格式来执行所述用零填充的不完全LU因式分解。
5.如权利要求1所述的系统，其中所述再因式分解器还可操作为在定义所述并行线程的配置中计算特别针对每一个等级的网格和块启动参数。
6.如权利要求1所述的系统，其中所述再因式分解器还可操作为启动所述并行线程，以此作为所述后续矩阵的每等级中的两个同时内核，其中一个所述内核可操作为更新所述等级中的乘数，并且另一个所述内核可操作为更新所述等级中的剩余元素。
7.如权利要求6所述的系统，其中所述再因式分解器还可操作为使用SMD处理器来启动所述并行线程。
8.一种SMD处理器，包括: 可操作为处理并行线程的通道； 可操作为控制所述通道中的线程处理的管线控制单元系统；以及 通过使用所述通道和所述管线控制单元来再因式分解方形输入矩阵的系统，包括: 矩阵生成器，其可操作为在下和上三角矩阵组合的归零稀疏模式中嵌入所述输入矩阵的置换形式，以便产生中间矩阵，其中所述下和上三角矩阵源于对在先矩阵的LU因式分解，所述在先矩阵与所述输入矩阵具有相同的稀疏模式，并且通过重排序处理而将填充和枢纽策略最小化；以及 与所述矩阵生成器关联的再因式分解器，其可操作为将并行线程传送到所述管线控制单元，以便将用零填充的不完全LU因式分解应用于所述中间矩阵。
9.如权利要求8所述的处理器，其中所述中间矩阵生成器可操作为将所述置换形式嵌入所述稀疏模式。
10.如权利要求8所述的处理器，其中所述再因式分解器还可操作为对所述方形矩阵执行符号分析，以便将代表其列的独立节点分组成代表所述节点之间的数据依存性的等级。
【文档编号】G06F17/16GK103914433SQ201310746800
【公开日】2014年7月9日 申请日期:2013年12月30日 优先权日:2013年1月9日
【发明者】马克西姆·瑙莫夫, 莎兰彦·柴特勒尔, 简龙升, 罗伯特·斯特佐得卡, 菲利普·范德尔莫尔斯赤 申请人:辉达公司