在GPU上对大型矩阵进行端到端LU因子分解

2880在GPU上进行大矩阵的端到端LU分解0YangXia，俄亥俄州立大学计算机科学与工程系，美国俄亥俄州哥伦布，xia.425@osu.edu0PengJiang，爱荷华大学计算机科学系，美国爱荷华州艾奥瓦市，peng-jiang@uiowa.edu0GaganAgrawal，奥古斯塔大学计算机与网络科学学院，美国乔治亚州奥古斯塔，gagrawal@augusta.edu0RajivRamnath，俄亥俄州立大学计算机科学与工程系，美国俄亥俄州哥伦布，ramnath@cse.ohio-state.edu0摘要0对于稀疏矩阵的LU分解是许多工程和科学问题（如电路模拟）的重要计算步骤。有许多工作致力于并行化和扩展该算法，其中包括针对GPU的最近工作。然而，由于高内存需求和数据依赖性，要在GPU上部署完整的稀疏LU分解工作流仍然具有挑战性。在本文中，我们提出了第一个完整的稀疏LU分解GPU解决方案。为了实现这个目标，我们提出了一个符号执行阶段的外存储实现，从而消除由于大型中间数据结构而引起的瓶颈。接下来，我们提出了在GPU上对Kahn算法进行拓扑排序的动态并行性实现。最后，对于数值分解阶段，我们通过与现有的实现方法相比，增加了大矩阵的内存限制，从而提高了并行度。实验结果表明，与从GLU3.0修改的实现相比，我们的外存储版本的加速比为1.13-32.65倍。此外，我们的外存储实现与GPU上优化的统一内存实现相比，加速比为1.2-2.2。最后，我们证明了我们为数值分解引入的优化是有效的。0未经收费，可以制作本作品的全部或部分的数字或实体副本，但不得为了盈利或商业优势制作或分发副本，并且副本必须带有本通知和第一页的完整引用。必须尊重他人拥有本作品组成部分的版权。允许带有署名的摘要。要进行其他复制、再版、发布到服务器或者分发到列表，需要事先获得明确的许可和/或支付费用。请向permissions@acm.org申请权限。PPoPP ’23,2023年2月25日至3月1日，加拿大魁北克蒙特利尔 ©2023版权归所有者/作者所有。出版权由ACM许可。ACM ISBN979-8-4007-0015-6/23/02...$15.00https://doi.org/10.1145/3572848.35774860CCS概念：  •  计算方法学  →  并行算法；  线性代数算法。0关键词：GPU加速，LU分解，内存限制01 引言0广泛的工程和科学计算应用需要解决大型线性代数方程组，即  �  ×  � =  �。解决此问题的直接方法涉及将矩阵  �  转换为两个矩阵：下三角矩阵  �  和上三角矩阵  �，使得  �  =  �  ×  �。这个转换是通过常被称为LU分解的方法完成的。在调用该过程后，可以通过解方程组来轻松获得解  �，这些方程涉及这两个三角矩阵，计算上更容易。对于许多应用（例如电路模拟[16,24, 30]），矩阵  �可能非常大且稀疏。因此，对于这些应用来说，大稀疏矩阵的LU分解成为一个关键的计算核。稀疏LU分解通常引入新的填充元素（即在  �  或  �中不是非零元素的元素）。这些新填充元素的位置在运行时之前是未知的。因此，LU分解涉及符号分解阶段，其中确定这些填充元素的数量（和可能的位置）。这一步之后是数值分解阶段，其中计算实际的非零值。由于LU分解的重要性，加速稀疏LU分解的研究已经进行了几十年[3-5, 8-11, 15, 19, 27,32]。随着异构和加速器计算在高性能计算（HPC）领域的主导地位，最近的几项研究特别考虑了GPU上的稀疏LU分解[6, 13,29]。然而，这些工作中的每一个只使用GPU加速计算的一部分，具体来说，只加速数值分解[19, 28, 32,36]或只加速符号分解（并且进一步限制为计数非零元素而不是它们的位置）[11]。因此，这些工作无法充分利用GPU的并行性。在本文中，我们提供了第一个高效的实现，在GPU上执行稀疏LU求解器的所有步骤，同时还能处理中间内存需求可能超过GPU可用内存的大矩阵。为了实现这个目标，我们解决了三个主要问题。第一个问题是符号分解步骤中的大内存使用。为了将实现扩展到大矩阵，我们提出了一种在GPU上进行符号分解的外存储实现，通过迭代执行符号分解。第二个问题是调度过程的并行化。在背景下，数值分解的实现是逐列进行的，并且需要检测列之间的依赖关系并确定其用于数值分解的顺序。先前的工作要么使用消除树[10, 38]，要么使用层次化算法[32,33]，但是在CPU上部署该过程。在本工作中，我们观察到该调度步骤本质上是一种拓扑排序，并提出了一种基于Kahn算法的高效GPU实现，其中使用了动态并行性功能。最后，以前的GPU实现[19,32]使用密集数据格式来加快在数值分解步骤中的数据访问。在这种情况下，每列需要O( �)的内存空间，并且可以存储在GPU上的列数（从而可以并行处理）受到限制。在处理大矩阵时（即当  �增大时），这可能成为一个重大挑战。为了解决这个问题，我们提出在矩阵大小增大时切换到稀疏数据格式。尽管稀疏数据格式上的数据访问（即为给定的列ID搜索行ID）效率较低，但我们通过使用基于二分搜索的访问算法仍然可以实现更好的性能，因为对并行列数的限制被移除了。我们使用SuiteSparse矩阵集合中的一组矩阵对我们的实现进行评估。我们总结我们的观察结果如下。首先，我们的GPU外存储符号分解可以支持大矩阵，并且与最近的一篇有效的CPU实现（即GLU3.0[32]）相比，加速比在1.13到32.65之间。其次，与用于符号分解的优化统一内存解决方案相比，我们的GPU外存储解决方案在1.2到2.2的速度提升范围内，因为它可以更有效地访问GPU设备内存。然后，我们证明了我们在数值分解上的优化进一步提高了大矩阵的性能，最多可提高3.33倍，这是由于内存效率导致的并行性增加。2890PPoPP ’23, 2023年2月25日至3月1日，加拿大魁北克蒙特利尔  杨霞，蒋鹏，Gagan Agrawal和Rajiv Ramnath0(a) 一个示例矩阵A0(b) 矩阵A的图表示0(c) 矩阵A的依赖图0(d) 不同级别的列的列ID0图1. 示例矩阵。 (a)对稀疏矩阵A进行符号分解，以第9行为例，红色圆圈表示新的填充项 (b) 矩阵A的图表示，红色线代表新的填充项 (c)矩阵A的依赖图，虚线将不同级别的列分隔开 (d)用于不同级别的列的列ID的表。02 背景 LU分解一个� × �矩阵�的形式为� = � ×�，其中�是一个下三角矩阵，�是一个上三角矩阵。LU分解实现通常涉及预处理、符号分解和数值分解的主要步骤。在预处理过程中，对行和列进行排列以提高数值稳定性并减少�和�矩阵中的填充项的数量[9, 10,27]。因为我们的实现主要涉及数值和符号分解阶段，所以我们在这里对它们进行回顾。具体而言，我们首先详细介绍符号分解，然后概述数值分解的GPU实现：GLU3.0[19,32]。最后，我们简要解释我们的工作的动机。02.1 符号分解和GPU实现0LU分解涉及一系列的基本行操作，其中枢轴行乘以非零标量并更新为另一行，该行在枢轴列中有非零项。因此，对于稀疏矩阵，LU分解通常会引入新的非零项，这被称为填充项。图1显示了一个示例矩阵�，该矩阵将在本文中使用。考虑图1(a)中的第9行，由于第9行在第5列有一个非零项，第5行需要对第9行进行行操作。结果产生了一个新的填充项(9,8)。如前所述，需要进行符号分解步骤来确定�和�矩阵的非零结构。新填充项的位置由以下定理确定：2900在GPU上对大型矩阵进行端到端的LU分解 PPoPP ’23, 2023年2月25日-3月1日，加拿大蒙特利尔 Yang Xia、Peng Jiang、Gagan Agrawal和Rajiv Ramnath0定理1.当且仅当存在一条从�到�的有向路径，并且路径中的中间顶点小于�和� [34]时，索引(�, �)处引入了一个填充项。0根据定理1，已经提出了几种算法来实现符号分解[14, 15,34]。在本节中，我们简要总结了fill2算法，该算法具有很高的并行性，因此适用于GPU。详细的过程如算法1所示。它使用数组����来指示已经访问过的顶点，将����(��������)设，执行初始化操作，然后对于被视为前沿的每个阈值，检查其邻居，更新邻居的状态，并将新的填充项添加到�(���,  :)或�(���:)，以及�������������(:)。随后，它将在第11行继续下一个阈值顶点。最近，已经提出了fill2算法的GPU实现来加速符号分解[11]。如算法1所示，从每个源行开始的符号分解是独立的。因此，我们可以从所有列进行并行遍历。但是，对于这种独立的遍历，每个节点需要�(�)的内存，这导致了总共需要�(�2)的内存。虽然使用分布式资源集合可以增加可用内存总量（之前的研究[11]使用了44个节点和264个GPU），但更好的内存效率将是可取的。还应该注意的是，他们的解决方案只计算每行中的新填充项的数量，这对于后续的数值分解步骤（他们在GPU上不实现）来说是不够的。02.2 数值分解和GPU实现0传统的数值分解方法是右向LU分解：� � 11 � 12 �  220� � � 11 � 220� = � � 11 �  � 21� 220�0其中，� 11、� 11是标量，� 11 = 1，� 21是大小为(� - 1) ×1的列向量，� 12是大小为1 × (� - 1)的行向量，� 22和�  22是(� 1) × (� - 1)的子矩阵。要计算矩阵�和�，我们可以首先找出� 11 =� 11，� 12 = � 12，� 21 = � 21 / � 11。然后，我们可以解出� �  22  = � 22  -  � 21 × � 12。如图所示，传统的右向0算法1基于定理1的填充2算法。该算法显示了src��行的过程。0输入：src-矩阵的第src行，�(:,:) -原始矩阵输出：填充后的矩阵L和U01: ����(:) = 0; 2: ����(���) = ���; 3:对于�(���,:)中的每个v，执行4: ����(�)= ���；05: 如果v < src，则06: 将v添加到�(���,  :)；07: 否则08: 将v添加到�(���,  :)；09: 结束条件判断010: 结束循环11: 对于阈值从0到src-1，满足����(���������) ==���的情况，将阈值添加到������������( :)；013: 对于�������������(:)中的每个前沿：014: 对于�(��������,  :)中的每个邻居进行循环：015: 如果����(��������) < src，则016: ����(��������) = ���017: 如果neighbor > threshold，则018: 将邻居添加到�(���,  :)或�(���,  :)；019: 否则020: 将邻居添加到��������������( :)；021: 结束条件判断022: 结束条件判断；023: 结束循环024: 结束循环；025: 交换(�������������(:), ���������������(:))026: 转到第13行；027:  结束循环0该方法先解出矩阵�的一行，然后解出矩阵�的一列，然后用�次迭代递归地解出矩阵。然而，这种方法存在顺序数据依赖，这会限制并行性。因此，为了解决这个问题，之前的工作[19,23,32]提出了一种混合基于列的右向算法，可以利用列级并行性。该算法的过程如算法2所示。对于每一列�，第一步是计算当前列的�部分，如2-6行所示。然后，向右查找所有满足条件��(�,�) ≠0的列�，这些列被称为列�的子列。然后，可以并行地分解列�的所有子列，如8-14行所示。为了调度不同列的分解顺序，算法需要检查不同列之间的依赖关系。例如，对于任何�(�, �) ≠0，我们可以得出结论列�依赖于列�。具体来说，这是因为列�是列�的子列，算法在分解列�时会读取��(�,�)，如算法2的8-14行所示。还有其他的依赖关系，但为了简洁起见，我们将读者引用到早期的出版物[32]。GLU3.0然后推导出2910PPoPP '23, 2023年2月25日-3月1日，加拿大蒙特利尔 Yang Xia、Peng Jiang、Gagan Agrawal和Rajiv Ramnath0算法2 混合基于列的右向算法。0输入： �� - 符号分析后的非零填充矩阵A01:  对于  �  = 1;  �  < =  � ;  � ++  do 2:对于  �  = j+1;  �  < =  � ;  � ++  do03: 如果  �� ( �, � )  ≠  0  则04: �� ( �, � )  =  �� ( �, � )/ �� ( �, � )05: 结束如果06: 结束循环07: 对于  �  = j+1;  �  < =  � ;  � ++  do08: 如果  �� ( �,� )  ≠  0  则09: 对于  � =j+1;  �  < =  � ;  � ++  do010: 如果  �� ( �, � )  ≠  0  则011: �� ( �,� )  =  �� ( �,� ) −  �� ( �, � ) ×  �� ( �,� )012: 结束如果013: 结束循环014: 结束如果015: 结束循环016: 结束循环0列的所有依赖信息并构建依赖图。图1（b）显示了矩阵�的依赖图。图中的边（ �, �）表示列�依赖于列�。基于依赖图，算法将列分组为级别，使得级别内的列彼此独立，因此可以并行因式分解。图1（c）显示了矩阵�的级别信息。例如，列1、2、3、6和7相互独立，它们的处理可以并行进行。确定这样的级别的过程本质上是一个拓扑排序，但也被称为分层化。GLU3.0的研究发现，潜在的并行性在不同级别之间不断变化。一般来说，他们将级别分为三类。在因式分解的开始阶段，级别是“类型A"级别。这些级别通常具有大量可并行化的列，而每列只有很少的关联子列。因此，他们使用一个线程块来因式分解一列，一个warp分配给一个子列。相反，类型C级别位于因式分解过程的末尾。在这个阶段，级别只有有限数量的列，而每列通常有大量的子列。为了利用子列的并行性，每个子列分配一个线程块，每个列分配一个内核调用而不是块。过渡阶段中的B级别既有大量的列，同时列也有许多子列。02.3 目前工作的局限性和挑战0先前的研究工作证明了GPU可用于数值因式分解阶段[19，32]或部分符号因式分解阶段[11]。然而，还没有提出完整的GPU解决方案来解决LU分解问题。更具体地说，最相关的先前工作[11]提供了基于GPU的符号执行的部分解决方案，即调度阶段仅在CPU上执行。此外，该工作没有集成符号和数值执行阶段。基于GPU的设备内存越来越大，并且越来越多的应用程序被移植到GPU上的事实，我们提出开发第一个端到端的稀疏矩阵LU分解实现。实现这一目标的主要挑战是：1）在中间阶段的大内存需求，2）涉及数据依赖的计算。0图2. 用于大稀疏矩阵的GPU LU分解的总体框架。0在本节中，我们首先介绍我们端到端GPU解决方案的框架。然后，我们展示了我们的带外GPU实现，用于执行符号因式分解以及用于调度的拓扑排序的GPU实现。最后，我们提供了用于在数值因式分解阶段增加并行性的优化。03 GPU LU分解的端到端实现0在GPU上端到端大矩阵LU分解PPoPP’23，2023年2月25日至3月1日，加拿大蒙特利尔03.1总体框架我们稀疏LU分解的总体框架如图2所示。按照惯例，我们首先执行某些预处理步骤，即行和列的排列，目的是减少填充和提高数值稳定性。然后，我们分两个阶段进行符号因式分解。之后，在GPU上进行并行化版本的层次化——这一步的输出用于在数值因式分解阶段调度并行计算。最后，进行数值因式分解阶段。2920在实现中，我们的新贡献是在行数变大时切换到使用稀疏数据格式，从而减少内存需求并增加并行性。0����� _ ����  =  � /( �  ×  � ) .0因此，迭代次数��� _ ���� ，表示为  � / ����� _ ����。我们还注意到，在先前的基于GPU的符号因式分解工作[11]中存在两个问题（该工作仅执行了符号因式分解的一部分），对于我们的端到端带外GPU实现：首先，它只计算了LU分解期间产生的新填充的总数。这对于后续的数值因式分解是不够的，即，如算法2所示，数值因式分解算法需要：1）每行的新填充数量，2）每个新填充的确切位置。其次，[11]使用了固定的值来0（a）pre20（b）audikw_10图3. 迭代次数（x轴）对应的边界大小（y轴）0����� _ ���� 。为了确保中间数据结构驻留在GPU上，它将使用保守的����� _����值，并相应地限制并行度。0算法3 用于符号因子化的GPU外存储器实现。01:  num_iter = n / chunk_size 2:  �计算每行中的填充项数量 3:  for  iter = 0; iter >>  ( ����� _ ����)05:  end for 6:  � ��� _ ����� 记录每行中的非零元素数量 7:  prefix_sum( � ��� _����� _ � ) 8:  为因子化的矩阵分配内存。 9:  �在内核symbolic_2中计算填充项的位置 10:  for  iter = 0; iter >>  ( ����� _ ���� )012:  end for0为了解决[11]中的第一个问题，我们的实现包含两个阶段。对于第一个阶段，我们仅计算被因子化矩阵中每行的非零元素的数量。在每次迭代中，我们并行地启动内核��������_1来计算�����_����行的非零元素的数量，这在算法3的第3-5行中显示。��������_1的实现是从最近的出版物[11]修改的。每行中的非零元素的数量存储在数组����_����中。由于我们使用压缩稀疏行（CSR）数据格式来存储被因子化的矩阵，我们对数组����_����应用了GPU实现的前缀求和，以获得每行的起始位置和非零元素的总数（第7行）。有了这些信息，我们就能够为被因子化的矩阵分配设备内存（第8行）。最后，对于第二阶段（第10-12行），对于每次迭代，我们启动内核��������_2来处理�����_����行。��������_2的实现与��������_1相比的主要区别在于，一旦我们找到一个填充项，我们还会存储其位置。这也是为什么需要提前执行��������1来找到空间需求的原因。对于第二个问题，我们观察到随着源行标识符值的增加，每个源行的内存需求也增加。具体来说，这是结果的2930PPoPP ’23，2023年2月25日至3月1日，加拿大蒙特利尔，魁北克省，加拿大杨侠，姜鹏，Gagan Agrawal和Rajiv Ramnath0根据定理1，即随着源行标识符值的增加，需要考虑的中间顶点的数量会更多，因为中间顶点需要具有较小的值。我们通过示例矩阵pre2和audikw_1进行验证，并在图3中显示了结果 -具有比源行标识符较小的标识符并且需要主动考虑的中间顶点称为frontiers。从图中可以观察到，对于最后几次迭代，frontiers的数量通常很大，否则很小。因此，我们提出了一种动态并行性分配实现，总结如算法4所示。我们首先将行分成两部分：第一部分包含�1行，第二部分包含剩余的（�2）行 -这里，�1是在我们看到“大”数量的frontiers之前的行数，我们将其定义为我们遇到的最大frontiers数量的50％。关键区别在于我们计算并使用这两部分的不同����� _����。具体来说，对于第一部分，每行对应的存储器需求较小，因此我们分配一个较大的����� _����以增加并行性。然后，我们计算每部分的迭代次数，这在算法的第1-2行中显示。最后，对于每个部分，我们迭代地启动内核��������_1来计算非零元素的数量，这在第4-6行和第8-9行示。第二阶段的过程类似，省略在算法中。需要注意的是，在进行这种优化时，可以尝试使用多于2个阶段，但这也意味着更多的内核启动。0算法4用于符号因子化的动态并行性分配的GPU外存储器实现。01:  num_iter_1= n1 / chunk_size_1 2:  num_iter_2= n2 / chunk_size_2 3:  �  计算每行中填充项的数量4:  for  iter = 0; iter < num_iter_1; iter++  do 5:symbolic_1 <<<>>>  ( ����� _ ���� _1 )06:  end for 7:  ... 8:  for  iter = 0; iter >>  ( ����� _ ���� _2 )010:  end for 11:  �  第二阶段类似..12:  ...03.3 基于GPU的并行调度过程根据依赖图，之前的工作计算每列的级别编号如下：0����� ( � )  =  ��� (− 1 ,����� ( � 1 ) ,����� ( � 2 0这里 �1，�2，...是节点�的子节点。该过程的本质是串行的，因为不同列之间存在依赖关系：列�的级别编号取决于列�1，�2等的级别编号。因此，0先前在LU分解方面的努力都是在CPU上执行levelization，并且因此无法在GPU上实现端到端的LU分解。为了在GPU上并行化levelization，即本质上是一种拓扑排序，我们首先注意到先前将此内核映射到GPU的努力的限制。已经有一些通用的GPU图处理系统[21,41]报告说他们可以支持拓扑排序。然而，它们中没有一个明确地为这个算法进行优化。一些出版物[37]明确报告了GPU拓扑排序的实现，但他们使用CPU来启动内核，因此不能充分利用GPU提供的并行性。0算法5 在GPU上的并行levelization实现。01:  __global__ void Topo( ){ 2:  level_num = 0 3:  �  � _ �����是所有没有入边的节点的队列 4:  cons_queue <<<>>> (...) 5:level_num++ 6:  while  qsize > 0  do 7: update <<<>>>(...);08: qsize = 0;09: cons_queue <<<>>> (...)010: level_num++;011:  end while 12:  } 13:procedure  Levelization 14:cons_graph <<<>>> (...)015: cnt_indegree <<<>>> (...)016: Topo <<<>>> (...);017:  end procedure0在本工作中，我们提出了一种纯GPU实现，使用了cuda的动态并行性特性。与[37]相比，我们的动态并行性实现具有以下优点：首先，它避免了CPU和GPU之间的同步和数据传输。其次，使用在GPU内部调用的函数，大大减少了内核启动的开销。详细的过程显示在算法5中。该方法假设已经构建了一个依赖图�，基于第2节中提到的方法（第14行）。然后，我们使用内核cnt_indegree（第15行）计算每个节点的入度。实际的拓扑排序过程开始于我们启动内核Topo_Sort（第16行）。在内核中，我们首先创建d_queue数据结构，它表示要处理的节点。最初，这个集合包括所有没有入度的节点 -这是通过使用子内核cons_queue（第4行）来完成的。更具体地说，此内核检查每个节点的入度，如果一个节点的入度为零，该方法将节点放入d_queue并将级别号设置为0。在循环的每次迭代中，我们启动一个子内核update来更新d_queue中节点的邻居的入度值。然后，我们为这样的邻居构建一个新的d_queue（cons_queue过程，第9行）。我们增加级别2940大规模矩阵在GPU上的端到端LU分解 PPoPP '23，2023年2月25日至3月1日，加拿大蒙特利尔，魁北克省，加拿大0在第10行的每次迭代结束时的数字。我们重复这个过程，直到没有没有入度为零的节点，即 � _ ����变为零。与我们所知道的关于这个问题的现有工作[37]相比，我们的改进在于在GPU内部调用函数，而不是使用CPU来启动内核。虽然由于基准代码不可用，直接比较是不可能的，但我们可以期望显著的改进，因为内核启动的开销被消除了。顺序拓扑排序的计算复杂度是 � ( �  +  � )，其中 � 是节点的数量，�是边的数量。并行拓扑排序的跨度（最长执行步骤）是层级的数量。0算法6 二分搜索以访问As(i,j)。0输入: ��� _ �� � ���  - CSC格式的列偏移。��� _ ��  -CSC格式的行ID。���  - CSC格式的值。01：... 2：fs =col_offset[j] 3：fe =col_offset[j+1]4：while fe >= fs do5：mid = (fs + fe) / 206：if row_id[mid] == i)  then07：As(i,k) = As(i, k)- val[mid]  ×  As(j,k)08：跳出循环09：else if row_id[mid] > i then010：fe = mid - 1011：else012：fs = mid + 1013：end if014：结束循环03.4 通过消除数值因子分解的内存限制来增加并行性0之前在GPU上的数值因子分解实现-特别是GLU实现[19，32，33]-都使用了矩阵�� ( �, � )的稠密格式（参见算法2）。进一步解释一下，在算法2中，我们需要搜索大于列ID  � 的行ID  �。因此，当我们使用密集格式时，我们可以高效地访问数据，因为位置是直接的�。然而，我们观察到这会增加总内存需求，限制可以存储在块中的行数，从而减少并行性的数量。假设总的可用设备内存为�，并且最大可并行化的列�可以计算为：0�  = �0�  ×  ������  ( ���� ���� )0，其中�是顶点数量。由于我们使用一个线程块来执行一列的数值因子分解，�表示最大可能的并发线程块。从这个表达式可以看出，当�变大时，�可能会小于最大的并发线程数（表示为�� _���）。这会导致实现无法利用足够的GPU并行性。为了解决这个问题，当我们确定�大于�时，我们建议采用压缩的稀疏列（CSC）数据格式来表示�� ( �, � )。0表1 Nvidia Tesla V100的规格0GPU Tesla V1000#SM 800FP32 CUDA核心/ GPU 51200内存接口4096位HBM20寄存器文件大小/ SM（KB）655360每个线程的最大寄存器数2550共享内存大小/ SM（KB）可配置最多96 KB0最大线程块大小10240不断增大。最终，当�变得非常大时，�可能会小于最大并发线程数（表示为�� _���）。这将导致实现无法利用足够的GPU并行性。为了解决这个问题，我们建议在�� ( �, �)中采用压缩稀疏列（CSC）数据格式，当我们确定�大于�时。0�� _ ��� × ������  ( ���� ���� )。然而，这种更改的挑战在于对于给定的列ID  �，我们无法直接获得大于列ID  � 的行ID  �。因此，我们利用CSC格式中的升序属性，并执行二分搜索以找到行ID  � 大于列ID  � 1的位置。详细过程如算法6所示。在该算法中，�表示列ID，�表示行ID，并且��� _ �� � ��� [j]和��� _ �� � ���[j+1]之间的索引是有序的。��表示要搜索的最小可能索引，初始化为��� _ �� � ��� [j]，��表示要搜索的最大位置索引，初始化为��� _ �� � ���[j+1]。在每次迭代中，我们将索引的中间值���与�进行比较。如果���中的行ID与�相同，则我们找到了索引，即���。否则，如果���中的行ID大于�，则��更新为���−1。或者，如果���中的行ID小于�，则��更新为���−1，并继续搜索。04评估和性能研究在本节中，我们展示了一组大型矩阵上的实验结果，以证明我们的端到端GPU实现的有效性。我们首先展示了我们的实验环境以及所选输入矩阵的特点。然后，为了展示我们的用于符号因子分解的基于外存的GPU实现的有效性，我们将其与从GLU3.0[32]修改的并行实现和优化的统一内存实现进行比较。接下来，我们评估了引入的优化所带来的好处，包括在最大矩阵的数值因子分解期间使用稀疏数据表示。01 请注意我们的CSC格式是有序的。2950PPoPP'23，2023年2月25日至3月1日，加拿大蒙特利尔，QC，Yang Xia，Peng Jiang，Gagan Agrawal和Rajiv Ramnath0表2.符号因子分解的内存要求超过GPU内存大小的输入矩阵。0矩阵缩写 n nnz nnz/n0g7jac200sc G7 59310 837936 14.10rma10 RM 46835 2374001 50.70pre2 PR 659033 5959282 9.00inline_1 IN 503712 18660027 37.00crankseg_2 CR2 63838 7106348 111.30bmwcra_1 BMC 148770 5396386 36.30crankseg_1 CR1 52804 5333507 101.00bmw7st_1 BM7 141347 3740507 26.50apache2 AP 715176 2766523 3.90s3dkq4m2 S34 90449 2455670 27.10s3dkt3m2 S33 90449 1921955 21.20onetone2 OT2 36057 227628 6.30rajat15 R15 37261 443573 11.90bbmat BB 38744 1771722 45.70mixtank_new MI 29957 1995041 66.60Goodwin_054 GO 32510 1030878 31.70onetone1 OT1 36057 341088 9.50windtunnel_evap3d WI 40816 2730600 66.904.1 实验设计0环境：我们在Nvidia TeslaV100上进行了实验。GPU的规格如表1所示。GPU连接到运行频率为2.4 GHz的Intel(R) Xeon(R) CPUE5-2680（2013年的Ivy Bridge）-CPU包含14个物理核心，并为每个核心提供2个超线程，用于我们的基准实现。我们的实验主机内存大小为128GB。我们实验的主机操作系统是CentOS Linux release7.4.1708（Core）。我们的GPU实现基于CUDA11.2工具包，使用NVCC V11.2.152编译我们的程序。输入矩阵：我们从SuiteSparse矩阵集合[7]中选择了18个矩阵进行详细研究和分析。选择这些矩阵是因为它们可以进行LU分解，并且根据我们的验证，中间数据结构的内存需求超过了Nvidia TeslaV100的设备内存大小。换句话说，对于这些矩阵的每一个，符号因子分解不能在GPU上执行，除非进行显式的数据移动或使用统一内存。矩阵的规格如表2所示。我们的实验使用float作为数据类型。基准选择：我们主要将我们的基于外存的GPU实现的性能与从GLU3.0[32]修改的并行实现进行比较，GLU3.0是一种最近的高效实现。对于符号部分，因为我们的实现直接从Gaihre等人的代码开始，并在功能上进行改进（即计算数值阶段的非零位置并防止较大数据的内存不足），而不是性能，性能的比较将没有意义。其他用于数值阶段优化的工作并没有始终提供代码，此外，数值阶段是完整代码的执行时间的一小部分。最后，我们0与另一种设计选择进行了广泛比较，即使用统一内存。04.2 与修改后的GLU3.0实现的比较此比较的结果显示在图4中。执行时间按符号分解和数值分解的时间分解。我们注意到加速比（整个执行过程）在1.13-32.65的范围内。可以看出，我们的基于外部存储器的GPU实现与GLU3.0实现之间的差异主要来自符号分解阶段。尽管多核CPU和GPU之间的相对性能差异很大，但GPU的加速比似乎取决于每行非零元素的数量 ��� / �。当这个比值较大时，加速比往往较大。例如，矩阵WI和MI都具有最高的比值��� / �，并且拥有最高的加速比，而AP和OT2则处于相反的情况。这与一般观察结果一致，即随着计算变得（相对）密集，GPU变得更高效。04.3 与统一内存解决方案的比较我们进一步将我们的基于外部存储器的GPU实现与统一内存实现进行比较。首先，我们注意到即使统一内存解决方案也受到CPU主内存大小的限制。因此，对于这个实验，我们选择了18个矩阵中的7个矩阵，这些矩阵的中间数据大小可以适应CPU主内存，但不能适应GPU设备内存。具体来说，这些矩阵是表2中 �最小的7个矩阵，所有这些矩阵的行数少于41000。对于我们的实验，我们进一步调整了统一内存实现并尝试了不同的优化。我们发现预取中间数据结构可以提高效率。将我们的实现与统一内存实现（启用预取）进行比较的结果显示在图5中。与前一个图中一样，执行时间按符号执行和数值分解阶段分别分解。我们可以看到，我们的基于外部存储器的GPU实现比优化的统一内存实现高效1.06-2.22倍。进一步的分析显示，在相对密度较高的矩阵（例如WI和MI）中，统一内存实现非常有竞争力。另一方面，对于密度最低的R15和OT2，统一内存开销较大。这符合我们的预期，即计算量越小，页面错误的影响越大。我们还创建了一个没有任何预取的统一内存实现版本-这个版本仅限于符号分解。我们将我们的基于外部存储器的GPU实现的符号分解阶段的执行时间与统一内存实现进行了比较，结果显示在图6中。如图中所示，没有预取的情况下，统一内存实现的性能较差。他们2960大规模矩阵在GPU上的端到端LU分解 PPoPP ’23，2023年2月25日至3月1日，加拿大蒙特利尔0图4. 基于外部存储器的GPU实现和修改后的GLU3.0基线的归一化端到端执行时间（符号执行和数值分解分开计时）。0图5. 基于外部存储器的GPU实现和统一内存GPU基线的归一化端到端执行时间（符号执行和数值分解分开计时）。0图6. 我们的基于外部存储器的GPU实现、带有预取和不带预取的统一内存实现的归一化符号执行阶段时间归一化结果。0相对性能对于密度较低的矩阵，如R15和OT2，变差。对三个版本之间性能差异的进一步阐述显示在表3中。主要0统一内存实现中按需分页的性能缺陷是GPU页面错误的开销。如表中所示，统一内存版本花费了大量时间来服务页面错误，而2970PPoPP ’23，2023年2月25日至3月1日，加拿大蒙特利尔，Yang Xia，Peng Jiang，Gagan Agrawal和Rajiv Ramnath0表3.无预取和有预取情况下的GPU页面错误组数和服务GPU页面错误时间百分比的比较。wp表示有预取，wo p表示无预取。0矩阵 # GPU错误数 wo p错误数 wp pc. wo p(%) pc. wp(%) pc. ooc(%)0OT2 16734 4638 78.37 56.60 0.060R15 17322 4392 86.21 65.46 0.150BB 19753 5798 4