GPU负载均衡的编程模型

790GPU负载平衡的编程模型0Muhammad Osamamosama@ucdavis.edu加利福尼亚大学戴维斯分校，加利福尼亚州戴维斯市，美国0Serban D. Porumbescusdporumbescu@ucdavis.edu加利福尼亚大学戴维斯分校，加利福尼亚州戴维斯市，美国0John D. Owens0jowens@ucdavis.edu加利福尼亚大学戴维斯分校，加利福尼亚州戴维斯市，美国0摘要我们提出了一种GPU细粒度负载平衡抽象，将负载平衡与工作处理分离，旨在支持静态和动态调度，并提供可编程接口以实现新的负载平衡调度。在我们的工作之前，在不规则问题上释放GPU的潜力的唯一方法是通过特定于应用程序的紧密耦合负载平衡技术进行工作负载平衡。通过我们的开源负载平衡框架，我们希望提高程序员在开发GPU上的不规则并行算法时的生产力，并通过允许快速尝试各种现有负载平衡技术来改进此类应用程序的整体性能特性。因此，我们还希望通过在我们的抽象中将负载平衡与工作处理的相关问题分离，更容易管理和扩展现有代码到未来的架构。0CCS概念：•计算方法→共享内存算法。0关键词：负载平衡，稀疏计算，GPU，调度01 引言图形处理单元（GPU）在常规细粒度并行问题（如一般矩阵乘法（GEMM））方面表现出色，通常是为此类问题而设计的。在常规问题中，相邻的线程具有类似或相同的工作负载，并且通常可以达到接近100%的峰值GPU理论性能。更具挑战性的是具有丰富细粒度并行性但不规则并行性的应用。在这种应用中，相邻的线程0分发声明“A”（经批准，发布无限制）。0PPoPP’23，2023年2月25日至3月1日，加拿大蒙特利尔，魁北克省，加拿大 © 2023版权由所有者/作者所有。ACM ISBN979-8-4007-0015-6/23/02。https://doi.org/10.1145/3572848.35774340在类似GPU这样的高并行机器上，以锁步方式运行的线程具有不同的工作负载，可能具有不同数量的工作，这使得在此类机器上实现高效率的工作成为一项重大挑战。以稀疏矩阵向量乘法（SpMV）为例，输入为稀疏矩阵A和密集向量x。SpMV计算输出向量y =Ax，是非规则细粒度并行的一个示例。与GEMM不同，SpMV中的稀疏矩阵可以在矩阵的行内含有不规则性：矩阵的行可以具有不同数量的非零条目。将每个GPU线程映射到一行的简单映射可以暴露这种不规则性，其中相邻的线程可能被分配不同数量的非零条目进行处理，导致同一线程束中的线程等待具有大量非零条目的线程。由此不规则性造成的不平衡问题，即当工作在并行执行者之间分布不均时，一些执行者处于空闲状态，而其他执行更多的工作，被定义为负载不平衡问题。目前的实现通过应用程序特定的负载平衡技术来解决GPU上的负载不平衡问题，旨在均匀分布工作，使每个线程获得相同数量的工作项以实现最大性能（例如，Merrill和Garland的负载平衡SpMV实现[20]）。这些负载平衡技术通常与应用程序本身紧密耦合。这些实现中的负载平衡组件既复杂又经常是应用程序性能的最重要贡献者。我们的工作将当今特定于应用程序的负载平衡算法概括为一个清晰、模块化、强大的抽象，可应用于许多复杂的非规则工作负载。在构建我们的抽象的过程中，我们确定了目前在GPU上稀疏、不规则应用程序中部署的常见负载平衡方法：如GraphIt [5]、Gunrock[29]和GraphBLAST[31]等应用程序特定框架；来自低级CUDA库（如ModernGPU[3]和CUB[24]）的技术；以及在应用程序中手工编码的负载平衡算法实现，如SpMV/SpMM [10, 14, 20]、三角形计数 [13,16]和广度优先搜索 [6,21]。我们展示了以简单、直观、强大的抽象方式，这些负载平衡调度可以扩展以支持比它们设计用途更广泛的不规则工作负载。我们通过使用基于稀疏线性代数的负载平衡来遍历数据中心图的核心来证明这一点。编写高性能的负载平衡代码是复杂的，主要是因为这段代码必须承担许多任务。除其他任务外，它必须从特定的数据结构中接收数据，对该数据执行用户定义的计算，并以负载平衡的方式调度该计算。我们抽象的关键洞察力是将工作负载映射（负载平衡任务）与工作执行（用户定义的计算）分开，其中我们将稀疏格式（如压缩稀疏行（CSR））映射到称为工作原子、瓷砖和集合的简单抽象组件。这些基本组件以可组合的C++范围和基于范围的for循环的形式表示，并用于构建负载平衡调度。然后，程序员可以使用这些API来构建负载平衡、高性能的应用程序和基元。以这种方式表达，我们可以重建现有的应用程序相关负载平衡技术，以解决不规则性，使其更加通用、可移植和可编程。我们的工作的贡献如下：01CUDA线程束（warp）是执行指令的32个线程的集合。线程束中的线程是无分歧的，并以单指令多线程（SIMT）的方式运行。0本作品根据知识共享署名4.0国际许可协议获得许可。800PPoPP ’23，2023年2月25日至3月1日，加拿大蒙特利尔，魁北克省，加拿大 M. Osama，S. D. Porumbescu和J. D. Owens0通过一个简单、直观、强大的抽象，可以将这些负载平衡调度扩展到支持比它们设计用途更广泛的不规则工作负载。我们通过使用基于稀疏线性代数的负载平衡来遍历数据中心图的核心来证明这一点。编写高性能的负载平衡代码是复杂的，主要是因为这段代码必须承担许多任务。除其他任务外，它必须从特定的数据结构中接收数据，对该数据执行用户定义的计算，并以负载平衡的方式调度该计算。我们抽象的关键洞察力是将工作负载映射（负载平衡任务）与工作执行（用户定义的计算）分开，其中我们将稀疏格式（如压缩稀疏行（CSR））映射到称为工作原子、瓷砖和集合的简单抽象组件。这些基本组件以可组合的C++范围和基于范围的for循环的形式表示，并用于构建负载平衡调度。然后，程序员可以使用这些API来构建负载平衡、高性能的应用程序和基元。以这种方式表达，我们可以重建现有的应用程序相关负载平衡技术，以解决不规则性，使其更加通用、可移植和可编程。我们的工作的贡献如下：01.我们提出了一种用于GPU上的非规则并行工作负载的新颖抽象。在高层次上，我们的抽象允许程序员以最少的代码开发稀疏的、非规则的并行算法，并实现高性能。02.我们设计并实现了一组直观的API，这些API可用于我们的开源GPU负载平衡框架，该框架基于使用CUDA-C++范03.通过实现一种新颖的基于协作组的负载平衡调度，我们展示了实现新的负载平衡调度的简易性，该调度在第 5.2节中描述，是之前的线程、warp和块级负载平衡调度[30]04.我们提供了与SuiteSparse矩阵集[11]相比，使用简单的启发式算法和3个GPU负载平衡调度的SpMV性能的最新结果，几何平均加速比为2.7×。02 设计目标我们的编程模型专注于细粒度嵌套数据并行性的广泛类别。负载平衡的任务级并行性需要不同的方法，超出了本工作的范围。本节重点介绍我们负载平衡抽象的设计目标：0实现高性能。首要目标是实现现有不规则应用程序的高性能负载平衡算法。我们的抽象不能以显著的开销或性能下降为代价。通过将我们的抽象的性能与现有硬连线实现的性能进行比较，我们来衡量我们在实现高性能方面的成功。0可组合和可编程的接口。重要的是，我们不希望将用户限制在接管较大系统的库接口上。程序员更喜欢采用适应其控制结构的新软件组件，而不是要求他们采用新的控制结构。我们希望允许用户（1）保持对GPU内核边界（内核启动）的控制，（2）能够添加新的负载平衡算法，以及（3）通过现有的负载平衡API组合新的负载平衡原语。我们通过比较我们的抽象的代码行数（LOC）与现有实现来衡量我们的工作的可编程性，并通过使用我们现有的API来实现新的负载平衡算法来展示组合性。0可扩展到新的应用。我们的目标是将特定于应用程序的负载平衡技术与新的非规则并行领域分离并扩展。我们的抽象力求促进将现有的负载平衡技术重用于新的应用程序。我们使用SpMV作为基准应用程序，使用三种不同的负载平衡技术来实现，其中一些技术以前用于实现并行图分析内核[5, 6, 10,29]。0促进优化探索。我们的抽象的一个关键目标是通过切换用于平衡工作的底层负载平衡算法来促进对给定应用程序的优化探索。我们希望鼓励用户尝试启发式算法和新的负载平衡技术，以发现对其应用程序需求最有效的方法。我们通过优化SpMV在多个不同负载平衡技术下的性能响应来衡量这个目标的成功。0非目标除了上述设计目标外，我们还定义了我们的非目标：0支持其他并行体系结构。尽管我们相信所学到的经验应该适用于其他并行体系结构，但我们明确针对NVIDIA的CUDA体系结构和编程模型[23]。我们的抽象的许多组成部分利用了CUDA的线程、瓦片和块的计算层次结构，这些层次结构映射到物理流多处理器，以及分配更多的工作量来充分饱和底层硬件的超订阅模型，这些在第5.2节中描述，以实现高性能。810GPU负载平衡的编程模型 PPoPP '23,2023年2月25日至3月1日，加拿大蒙特利尔0针对其他并行体系结构。尽管我们相信所学到的经验应该适用于其他并行体系结构，但我们明确针对NVIDIA的CUDA体系结构和编程模型[23]。我们的抽象的许多组成部分利用了CUDA的线程、瓦片和块的计算层次结构，这些层次结构映射到物理流多处理器，以及分配更多的工作量来充分饱和底层硬件的超订阅模型，这些在第5.2节中描述，以实现高性能。0多GPU支持。这项工作侧重于单个GPU的负载不平衡问题，不考虑多GPU单节点或多节点系统，尽管这些是未来工作的有趣方向。03 我们的负载平衡抽象0我们GPU负载平衡抽象的关键见解是将工作项映射到处理单元和工作执行之间的关注点分离。我们将我们的抽象分为三个关键概念（如图1所示），每个概念描述实现的不同方面：（1）定义工作；（2）在GPU线程、瓦片或块之间定义工作负载平衡；（3）定义平衡工作中每个线程的工作执行和计算。这种分离使我们能够清晰地划分工作在应用程序开发人员和负载平衡库开发人员之间的分工，并通过混合不同的负载平衡技术和稀疏非规则算法来促进优化的探索。Sidebar1给出了我们负载平衡抽象的动机的一个实际示例。03.1稀疏数据结构的输入我们从某种稀疏数据结构中开始表达我们的输入数据。这些数据结构的示例包括但不限于压缩稀疏行（CSR）和坐标（COO）格式。我们抽象的第一阶段的目标是将输入数据格式映射到一种通用的数据框架和词汇，该框架和词汇是下一阶段的输入。这个词汇有三个简单的组成部分，共同表达输入数据：01.工作原子，要安排到处理器上的单个工作单元（例如，稀疏矩阵的非零元素）。我们假设所有工作原子在执行过程02.工作瓷砖，逻辑实体，表示为一组工作原子（例如，稀疏矩阵的一行）。工作瓷砖在执行过程中可能具有不同的成本。正如我们在介绍中强调的那样，工作在逻辑上最适合并行化工作瓷砖，但通常在工作原子上并行化效率最高，并且工作瓷砖和工作原子之间的映射可能是昂贵且复杂的。03.瓷砖集，一组共同构成整个工作问题的工作瓷砖（例如，稀疏矩阵）。在我们的抽象中，瓷砖集中的瓷砖必须是独立的（因此可以在多个处理器上并行运行）。0稀疏格式和原子/瓷砖/瓷砖集之间的映射由用户定义。虽然我们尚未实现0Sidebar 1稀疏非规则工作负载负载平衡的现有主流方法的实际示例。0考虑在开源CUDACUB库[24]中提供的GPU上的SpMV实现。CUB实现并维护了Merrill和Garland在论文[20]中提出的SpMV算法。基于合并的SpMV在第5.2.1节中详细解释，是一种基于CSR的完全负载平衡的SpMV，其中每个线程获得相等份额的工作，工作量由矩阵行的总数和非零元素的总数之和定义。在参考文献中，这个高效的、最先进的实现需要1100行代码（或503行内核代码）跨越3个文件（不包括额外的234行代码用于分段修复步骤的第4个文件）。相比之下，该参考实现中SpMV的实际计算在一个for循环中表达，仅需4-5行代码！在将工作项以负载平衡方式映射到处理单元所需的代码行数与表达所需计算之间存在差异，这是我们工作的关键动机。此外，CUB实现专门用于SpMV算法，将其应用于其他算法，甚至在相同的计算域内，都需要进行重写。其中一个精确重写的例子是由杨等人完成的，他们将基于合并路径的负载平衡从SpMV扩展到稀疏矩阵稠密矩阵乘法（SpMM）实现[30]。这两部分的负载平衡算法相同，但应用于不同的计算，这促使了重用的需求。0虽然我们相信我们的映射抽象足够灵活，可以以适合负载平衡抽象的方式表达文献中各种各样的现有稀疏数据格式[12]，但我们尚未实现所有这些格式，以便它们适用于我们抽象的下一阶段的负载平衡。此外，我们的负载平衡库实现中已经包含了几种常见的稀疏格式（CSR、CSC、COO），用户可以简单地选择和使用它们，而无需实现它们。给定对原子/瓷砖/瓷砖集的映射，我们可以进一步实现一个负载平衡算法，该算法可以透明地并行化工作原子或瓷砖的计算。03.2通过在不同粒度（例如，瓷砖集、原子和瓷砖）上捕获工作负载的抽象来定义负载均衡，可以更容易地将计算均匀地分布在GPU的可用资源上。根据用户定义的输入瓷砖集及其关联的原子和瓷砖序列，以及用户选择的分区算法，我们的负载均衡阶段输出分配给处理器ID的原子和瓷砖的子序列（即，原子或瓷砖将被处理的位置）。136820PPoPP '23，2023年2月25日至3月1日，加拿大蒙特利尔0负载平衡调度器20稀疏数据结构迭代器表示负载平衡0工作执行0工作原子和瓦片0原子迭代器 = 0,1,2,3瓦片迭代器 = 0,1,2,3原子/瓦片 = 0,1,3,0 值 =1,3,6,20图1.负载平衡作为我们抽象的三个关键概念的简单流水线：（1）稀疏数据结构表示为迭代器，（2）负载平衡算法将工作分配到线程上，（3）用户定义的计算消耗平衡的工作并在每个线程上执行。0将子序列分配给处理器ID的结果分配对于有效地平衡处理单元之间的工作负载至关重要，并且通常是特定于问题和数据集的。用户必须指定必要的序列。理想情况下，一个oracle会获取这些序列并选择每个处理元素的最优子序列。寻找这样的oracle是一个开放的问题，因此我们提供了下一个最好的选择：用户可以从一组预定义的调度中选择并进行实验，并实现自己的调度。一般而言，负载平衡算法设计者必须在调度的成本和更好的调度的收益之间取得平衡。调度可以非常简单，只需将处理元素分配给具有任意数量原子的瓦片（例如，带有稀疏矩阵中任意数量非零元素的行）；或者可以更复杂/昂贵，以更全面的方式处理工作（例如，考虑具有不同数量非零元素的多个行的工作）。03.3 定义工作执行我们负载平衡抽象的最后一个组成部分是表达不规则并行计算本身。前一个阶段输入负载不平衡的工作并对其进行负载平衡；然后，该阶段通过在负载平衡的工作上执行计算来消耗该工作。计算可以表达的范围非常广泛，只受限于如何在CUDA内核中消耗表示为序列的负载平衡工作。由于该框架不控制内核，因此您可以在CUDA内核中编写的任何内容在我们的框架中也可以工作。例如，程序员可以表达对每个原子或每个工作瓦片进行的数学操作，或者构建合作算法，不仅消耗分配给每个线程的工作，还将结果与相邻线程组合以实现更复杂的算法，如并行归约或扫描。我们在我们的框架中实现的实际示例（请参见第4.3节和第5.3节）使用此抽象包括但不限于稀疏线性代数核，例如稀疏矩阵和稀疏张量收缩，以及以数据为中心的并行图0算法，例如基于邻域遍历内核的单源最短路径（SSSP）和广度优先搜索（BFS）。我们预计我们库的典型用户只会为抽象的这个阶段编写自己的代码，并使用已经是我们库的一部分的标准数据结构和负载平衡调度。但是，这些用户也可以实现自定义数据格式和负载平衡调度。04 高级框架实现0我们的GPU负载平衡框架使用C ++17和CUDA实现了第3节中描述的抽象。在我们的系统中，程序员使用CUDA / C++开发不规则并行算法，并实现新的负载平衡调度。根据我们的设计目标，即可组合的API，可扩展性和重用性，本节和下一节介绍了我们API的实现细节，以及如何使用它来开发促进框架内高性能负载平衡技术重用的新应用。我们还探讨了一种基于CUDA的合作组模型的新的负载平衡方法（第5.2节）。此外，我们还确定了我们的工作如何用于促进给定应用程序（如SpMV）的优化探索。04.1 实现稀疏数据结构0我们的框架将稀疏数据结构（例如COO，CSR，CSC）转换为工作原子，工作瓦片和瓦片集（第3.1节）使用简单的C++迭代器。 C++迭代器是指向一系列元素中的某个元素的对象，并使用一组运算符对该范围的元素进行迭代。例如，计数迭代器是一个迭代器，表示指向连续值范围的指针[1]。我们的框架要求用户使用C++定义三个重要的迭代器：（1）遍历所有工作原子的迭代器；（2）遍历工作瓦片的迭代器；和（3）遍历每个工作瓦片中原子数量的迭代器。（我们的库已经支持几种常见的稀疏数据结构。）使用这些迭代器，负载平衡调度可以确定并分发负载平衡的工作到11});20}27}28};830GPU负载平衡的编程模型PPoPP '23，2023年2月25日至3月1日，加拿大蒙特利尔01 // 原子和瓦片的简单迭代器。02 counting_iterator< int > atoms_iter( 0 , nnz);03 counting_iterator< int > tile_iter( 0 , rows);04 // 在tile i中迭代原子。05 auto  atoms_per_tile = make_transform_iterator(06 tile_iter,07 [tile_iter, row_offsets]08 __host__ __device__( const int & i) {09 return  (row_offsets[tile_iter[i +  1 ]] -010 row_offsets[tile_iter[i]]);0Listing 1.使用C++17在我们的框架中表达的压缩稀疏行（CSR）格式。CSR格式使用三个数组描述矩阵：（1）非零值的列索引；（2）行的范围（行偏移量）；（3）非零值。由于CSR数据结构不包含指向原子和瓦片（非零值和行）的数组，在上述列表中，我们将原子和瓦片迭代器定义为从0到非零值（nnz）的总数和从0到矩阵的总行数（行）的简单计数迭代器（第2–3行）。使用变换迭代器来表示每个瓦片id的原子-工作-瓦片（atoms-per-work-tile）迭代器，该迭代器对于CSR来说只是当前瓦片的行偏移量减去下一个瓦片的偏移量（第5–11行）。0底层硬件。列表1显示了我们的抽象如何将常用的CSR格式作为我们框架中的瓦片集来表示。04.2实现负载平衡调度也许最直接的调度是将每个工作瓦片调度到一个GPU线程上。尽管这种策略在存在显著的负载不平衡时效果不佳[3，10，21，26，30]，但我们在这里使用它作为例子来说明如何在我们的框架中定义负载平衡。输入是来自上一个阶段的三个迭代器以及一个原子和瓦片计数。然后，负载平衡算法开发人员实现tiles（）和atoms（）过程调用，返回要由当前线程处理的C++范围的瓦片和原子，从而创建分配的处理器ID和工作负载部分之间的映射。列表2显示了线程映射调度的完整示例。尽管这是一个简单的算法，但对于具有粗粒度并行性（每个瓦片的原子数量很小）的平衡工作负载（例如将稀疏向量乘以稠密向量）可以提供高性能。此外，我们的抽象不仅限于简单的调度算法，第5.2节提供了更复杂的负载平衡算法的示例。01 class  schedule_t {02 // 构造一个基于线程映射的调度。03 __host__ __device__04 schedule_t(atoms_it_t atoms_it,05 tiles_it_t tiles_it,06 atoms_it_t atoms_per_tile_it,07 size_t  num_atoms,  size_t  num_tiles) :08 m_atoms_it(atoms_it), m_tiles_it(tiles_it),09 m_atoms_per_tile_it(atoms_per_tile_it),010 m_num_atoms(num_atoms),011 m_num_tiles(num_tiles) {}012 // 在“this”线程中处理的瓦片范围。013 // 通过网格维度进行步幅。014 __host__ __device__  auto  tiles() {015 auto  begin = m_tiles_it(blockDim.x * blockIdx.x016 + threadIdx.x);017 auto  end = m_tiles_it(m_num_tiles);018 return  range(begin, end)019 .step(gridDim.x * blockDim.x);021 // 在“this”线程中处理的原子范围。022 __host__ __device__  auto  atoms(023 const  std:: size_t & tile) {024 auto  begin = m_atoms_per_tile_it[tile];025 auto  end = m_atoms_per_tile_it[tile +  1 ];026 return  range(begin, end).step( 1 );029 using  schedule_t = thread_mapped_schedule_t;0列表2.一个基于线程映射的负载均衡算法，使用从列表1中定义的原子和瓦片作为迭代器的C++范围表示。每个瓦片都映射到一个线程，其中线程ID对应于瓦片在瓦片集中的索引。一个瓦片内的所有原子都由线程顺序处理。在处理完一个瓦片后，线程映射到下一个瓦片，该瓦片通过将索引按内核的网格大小进行步幅。04.3 实现工作执行020}21}840PPoPP '23, 2023年2月25日至3月1日，加拿大蒙特利尔，M. Osama，S. D. Porumbescu和J. D. Owens0为了能够表达的多样性，用户现在可以指定多个负载平衡的工作域、基于范围的for循环，甚至融合多个计算以构建更复杂的算法在单个内核中。（3）可以使用更高级别的API来构建更简单的高级抽象，这些抽象拥有内核边界并提供更简单的API，但代价是灵活性。作为这一阶段的输入，用户使用负载平衡的C++范围来实现他们的计算。可以通过多种方式实现这一点，但最常见的模式之一是嵌套的基于范围的for循环，循环遍历所有分配的瓦片和原子范围。列表3展示了一个简单的示例，其中实现了使用CSR格式和基于线程映射的负载平衡算法描述的SpMV算法的CUDA内核。在这个示例中，每个线程内的外部for循环迭代稀疏矩阵（瓦片）的分配行，内部循环按顺序处理分配的非零元素（原子）在每行内。在第5.3节中，我们实现并讨论了更复杂的内核和计算。05 实施细节05.1灵活的、可组合的基于CUDA的范围我们的API支持负载平衡原语和应用程序的组合能力是我们框架中的有意设计选择，通过使用基于CUDA的C++范围来实现。我们的框架不拥有内核边界（内核启动），这迫使我们的API集中在内核中并且被包含在内核中。这使得程序员能够构建和维护自己的内核，同时仍然从我们框架的负载平衡能力中受益。这主要是使用带有CUDA的__device__关键字标记的设备范围C++函数和类来实现的。我们实现并暴露了几种不同类型的专门化范围，这些范围在实现负载平衡调度时特别有用。0•  step_range：迭代从开始到结束的范围0以步长step迭代的步骤范围。对于定义需要自定义步进范围或每个线程处理固定数量的工作项的负载平衡调度非常有用（可以使用step来定义）。0•infinite_range：从开始到无穷大的范围。在持续内核模式下定义负载平衡调度非常有用[32]，其中内核持续运行，直到所有工作完成或算法收敛。0•grid_stride_range：使用CUDA内核的网格大小，从开始到结束以步长step迭代的步骤范围的特殊情况。还支持block和step。02使用__device__关键字修饰的方法允许CUDA编译器生成可从内核内部调用的入口点。这允许代码从内核内部调用[23]。01 // 实现负载平衡的SpMV内核。02 __global__  void  spmv(const size_t  rows,03 const size_t  cols,  const size_t  nnz,04 const int * offsets,  const int * indices,05 const float * values,  const float * x,06 float * y) {07 // 配置负载均衡。08 // 输入：为CSR格式定义的迭代器。09 schedule_t  config (010 atoms_iter, tile_iter,011 atoms_per_tile_it,012 nnz, rows);013 // 使用基于范围的for循环消耗行。014 for  ( auto  row : config.tiles()) {015 type_t sum =  0 ;016 // 使用基于范围的for循环消耗原子。017 for  ( auto  nz : config.atoms(row))018 sum += values[nz] * x[indices[nz]];019 y[row] = sum;022 // 启动SpMV内核。023 constexpr size_t  blocks =  256 ;024 size_t  grid = (rows + blocks -  1 ) / blocks;025 spmv<<>>(rows, cols, nnz,026 offsets, indices, values, x, y);0在我们的负载均衡抽象中实现的稀疏矩阵向量乘法（SpMV），使用基于范围的嵌套循环。稀疏矩阵使用基于CSR的格式表示，其中x是密集输入向量，y是密集输出向量（y =Ax）。第9-12行使用第2个示例中实现的负载均衡调度和第1个示例中定义的迭代器来构建要处理的负载均衡工作。第14行和第17行显示每个线程内的for循环，它们迭代分配给稀疏矩阵的行，并按顺序处理每行中分配的原子。第18行显示在每个工作原子（非零）上执行的实际计算，第19行将结果写入密集输出向量y。0warp stride变体，其迭代步长为块或warp的大小。05.2 实现非平凡的负载平衡正如我们在第5.1节中所描述的，我们可以将现有的负载平衡技术解耦并表示为一组C++范围。为了说明这种抽象的潜力，我们首先解耦并表示一个称为merge-path的最先进的负载平衡算法，该算法先前用于平衡基于CSR的SpMV和SpMM[20,30]，并实现了三个额外的负载平衡算法（warp-，block-和group-mapped），所有这些算法都可以在我们的库中供程序员使用。我们的新的组映射算法是一个基于组的瓦片每组调度，其中一个组被定义为一个组850GPU负载平衡的编程模型PPoPP '23年2月25日至3月1日，加拿大蒙特利尔0作为任意大小线程集合（不限于warp或block大小）的集合的抽象。我们的组映射调度是使用CUDA的合作组编程模型[18]将瓦片每组一个的调度的泛化版本，其中组被定义为任意大小的组05.2.1 合并路径负载平衡。从稀疏矩阵的角度来看，合并路径假设矩阵中的每个非零元素和每个新行都是等量的工作量，然后将 nnzs + rows的工作均匀分配给一组工作线程。然后，每个线程在CSR矩阵的非零索引和行偏移内进行二维二分搜索，以找到行和非零元素的起始位置。然后，线程顺序处理从起始位置开始直到达到其分配的工作的末尾的行和非零元素[20]。我们将此算法表达为我们抽象中的一个负载平衡调度，分为两个步骤：（1）设置：C++调度类的初始化步骤计算每个线程的工作单位数，按上述描述进行二分搜索，并将每个瓦片和原子的起始位置存储在线程局部变量中。（2）范围：算法的第二个步骤构建每个线程要处理的范围，作为“完整”瓦片和“部分”瓦片[20]。如果线程的原子范围完全位于一个瓦片内，则它是“完整的”，并在一个简单的嵌套循环中处理。如果线程的范围跨越瓦片边界，则线程在一个单独的嵌套循环中处理其工作。由于我们将负载平衡方法（第4.2节和上面）与工作执行（第4.3节）解耦，因此我们可以使用这个合并路径实现来实现不仅SpMV还包括任何其他工作可以划分为瓦片和原子的算法，例如用于实现广度优先搜索的图邻域遍历算法[29]。更重要的是，合并路径调度现在不再局限于基于CSR的稀疏格式。支持其他格式只需要构建必要的稍微复杂的迭代器，这些迭代器能够计算每个瓦片的原子数（CSR实现在列表1中使用行偏移数组进行计算）。05.2.2 Warp-和block级负载平衡。warp级和block级负载平衡调度的目标是将每个warp或block分配给相等份额的瓦片，然后按顺序处理。每个warp或block内的可用线程并行处理瓦片中的工作原子。每个线程按照warp或block的大小跨越步长来处理新的工作原子，直到达到工作的末尾。不同处理单元之间的不平衡由硬件调度器处理。此调度器依赖于CUDA的超订阅模型，其中程序员可以启动比GPU在任何给定时间物理调度的warp或block数量更多的数目。当处理单元完成其工作时，新的处理单元将从超订阅池中调度出来[5,21]。0处理单元完成其工作后，将从超订阅池中调度新的工作单元[5,21]。05.2.3 组级负载平衡。组级负载平衡推广了warp级和block级调度。与以上要求组大小为warp或block大小不同，该方法利用CUDA的合作组（CG）编程模型[18]允许程序员指定动态大小的任意大小的组。在这些组内，CG模型允许对组的同步行为进行详细控制，以及简单的并行组级集合操作，如reduce或scan。我们利用这个强大的工具来实现一个通用的组级负载平衡调度，当组大小等于warp或block的大小时，我们可以免费使用上述warp级和block级调度。我们的调度将工作瓦片分配给一个组，每个组查看其平均份额的瓦片并计算每个瓦片的原子数量，并将其存储在一个临时内存（CUDA的共享内存）中。然后，组执行并行前缀和，这是一种广泛使用的并行算法，它输入一个数组并生成一个新数组，其中任何位置的元素是所有前面元素的总和[4]。我们使用这个前缀和数组有两个目的：（1）前缀和数组的最后一个元素表示组必须处理的工作原子的聚合数量，（2）前缀和数组中每个和的位置对应于正在处理的原子所属的工作瓦片。调度的设置阶段在临时内存中为每个组建立前缀和数组，并且调度的范围循环以每个线程返回要处理的原子。如果需要，通过一个简单的get_tile(atom_id)操作获取对应的瓦片，该操作在前缀和数组中执行二分搜索以找到对应于正在处理的原子的瓦片。依赖于CG模型的这种负载平衡调度具有配置组大小的独特优势（实际上是直接映射到硬件的软件结构），可以根据问题的形状和底层硬件体系结构进行配置。例如，针对GPU的目标，其中warp大小不是32个线程（AMD的GPU体系结构支持64个线程的warp大小[2]），现在可以通过一个简单的编译时常数实现，或者将组大小配置为与问题的结构完全对齐。05.3 应用空间我们的工作定义（第3.1节），可组合的API（第5.1节）和多种复杂的高性能负载平衡调度（第5.2节）共同提供了一个多功能和可扩展的框架，为应用程序特定的优化提供了充足的空间。在列表3中，我们已经展示了如何使用我们的框架实现SpMV算法13}14}17}181101001,00010,000100,0001,000,00010,000,000100,000,0000.0010.0020.010.020.10.2121020100860PPoPP '23年2月25日至3月1日，加拿大蒙特利尔 M. Osama, S. D. Porumbescu和J. D. Owens01 // ... 在CUDA内核内部。02 // 循环遍历所有已分配的行。03 for  ( auto  row : config.tiles()) {04 // 循环遍历矩阵B的所有列。05 for  ( auto  col : range( size_t ( 0 ), B.cols)06 .stride( size_t ( 1 ))) {  /// < 新循环07 float  sum =  0 ;08 // 循环遍历所有已分配的非零元素。09 for  ( auto  nz : config.atoms(row))010 sum += values[nz] * B(nz, col);011 // 将总和输出到Matrix-C。012 C(row, col) = sum;0列表4.在列表3中引入的简单循环包围SpMV，允许我们表示稍微复杂的SpMM负载平衡计算。0我们的框架。一个简单而自然的扩展是实现稀疏矩阵矩阵乘法（SpMM）。列表4显示了所需的微小更改，它在现有代码周围添加了对矩阵B的列的另一个循环，以实现SpMM。这种实现也可以扩展到支持Gustavson的通用稀疏矩阵矩阵乘法（SpGEMM），使用两个内核和一个分配阶段；第一个内核将计算用于为输出稀疏矩阵分配内存的输出行的大小，而第二个内核将执行乘法累加。除了稀疏线性代数，我们可以使用我们的框架来解决其他领域的应用。列表5使用我们的组级负载平衡调度实现了图形基元的单源最短路径（SSSP）。SSSP在GPU上的性能在很大程度上受到良好的负载平衡的限制[5,29]，但如果程序员从我们的库中选择一个负载平衡调度，负载平衡的详细信息将完全隐藏。而且，在一个应用领域（例如，稀疏线性代数）中使用的相同调度在这个不同的应用领域中很容易重复使用。06 评估0我们的目标是展示我们的基于负载平衡抽象构建的框架能够在稀疏不规则问题上实现高性能和更好的可编程性。我们下面的评估使用我们的SpMV实现作为与NVIDIA的(开源)CUB库和生产(闭源)cuSparse库提供的最先进实现相比的基准。我们考虑(并实现)了其他几个应用程序进行评估，包括SSSP、BFS和SpMM。我们发现它们得出了类似的高级结论。因此，我们这里的评估重点是SpMV。我们的测试集包括约整个SuiteSparse矩阵集合[11]，涵盖广泛的稀疏矩阵。01 // ... 在CUDA内核中。02 // 循环处理所有分配的边。03 for  ( auto  edge : config.atoms()) {04 auto  source = config.get_tile(edge);05 // G是图数据结构06 auto  neighbor = G.get_neighbor(source, edge);07 auto  weight = G.get_edge_weight(edge);08 float  source_dist = dist[source];09 float  neighbor_dist = source_dist + weight;010 // 检查目标节点是否已声明为某个节点的子节点。011 // 声明为某个节点的子节点。012 float  recover_distance =013 atomicMin(&(dist[neighbor]), neighbor_dist);014 // 将邻居添加到frontier。015 if  (neighbor_dist < recover_distance)016 out_frontier[neighbor