0BenchCouncil基准，标准和评估交易1（2021）1000070ScienceDirect提供目录列表0BenchCouncil基准，标准和评估交易0期刊主页：https://www.keaipublishing.com/en/journals/benchcouncil-transactions-on-benchmarks-standards-and-evaluations/0基于MPI和CUDA的并行稀疏近似逆预处理算法 �0Yizhou Wang，Wenhao Li，Jiaquan Gao �0江苏省计算机与电子信息学科重点实验室，南京师范大学计算机与电子信息学院，中国南京2100230A R T I C L E I N F O0关键词：稀疏近似逆预处理CUDA GPU MPI0A B S T R A C T0在本研究中，我们提出了一种基于MPI和CUDA的高效并行稀疏近似逆（SPAI）预处理算法，称为HybridSPAI。对于HybridSPAI，它优化了最新的静态SPAI预处理算法，并从单个GPU扩展到多个GPU，以处理大规模矩阵。我们做出了以下重大贡献：（1）提出了一个基于MPI和CUDA的优化静态SPAI预处理器的通用并行框架；（2）对于预处理器的每个组件，建立了一个决策树来选择计算它的最佳内核。实验结果表明，HybridSPAI是有效的，并且在两个公共库中优于流行的预处理算法和最新的并行SPAI预处理算法。01. Introduction0已经证明了稀疏近似逆（SPAI）预条件0tioners can effectively accelerate the convergence rate of Krylov sub- spacemethods, such as the generalized minimal residual method (GM- RES) [ 1 ]and the biconjugate gradient stabilized method (BiCGSTAB) [ 2 ].Moreover, compared with the incomplete factorization precon- ditioners [3 – 6 ] and the factorized sparse approximate inverse (FSAI)preconditioners [ 7 – 10 ], SPAI preconditioners neither require exces-sively sparse matrix–vector multiplication operations nor take care of therisk of breakdowns that can be encountered by FSAI precondi- tioners [ 11]. Consequently, SPAI preconditioners have attracted much attention [ 12 –17 ].0近年来，图形处理单元（GPU）已成为0important resource for scientific computing because of their many corestructures and powerful computation efficiency, and have been used astools for high-performance computation in a lot of fields [ 18 – 21 ]. As weknow, the cost of constructing SPAI preconditioners is commonly veryexpensive for large-scale matrices, because the memory requirements tostore them, and the computation requirements to calculate them areapproximately the scale with the square to third power of the number ofnonzeros in each row.0随着图形处理单元（GPU）的出现，许多研究0have been conducted to accelerate the construction of SPAI precon-ditioners on the GPU architecture, and many parallel preconditioningalgorithms [ 11 , 22 – 26 ] are proposed. Based on the degree of freedomused, SPAI preconditioner generation is classified as static (a priori)0� 该研究得到了中国自然科学基金（批准号61872422）和江苏省自然科学基金（批准号BK20171480）的支持。� 通讯作者。邮箱地址：1966224230@qq.com（Y.Wang），917339495@qq.com（W. Li），springf12@163.com（J. Gao）。0or adaptive. In this paper, we focus on optimizing a latest static SPAIpreconditioning algorithm and extend it from one GPU to multiple GPUs.There has existed some work about static SPAI preconditioners on GPU [11 , 27 ], but the detailed implementations never be given and the sourcecode is not public. Furthermore, He and Gao et al. propose two static SPAIpreconditioning algorithms on GPU, called SPAI-Adaptive [ 28 ] andGSPAI-Adaptive [ 29 ], and give their imple- mentation details. The twoalgorithms are verified to be effective for large-scale matrices. In thisstudy, inspired by Gao’s work, we further investigate how to highlyoptimize the static SPAI on multi-GPUs instead of only single GPU in thispaper. We propose an optimized SPAI preconditioning algorithm based onMPI and CUDA, called Hy- bridSPAI. Compared to a latest static SPAIpreconditioning algorithm, the proposed algorithm has the followingdistinct characteristics. First, a general parallel framework based on MPIand CUDA is presented to optimize the static SPAI preconditioner, and isextended from one GPU to multiple GPUs. For each GPU, it operates sameprocedures as shown in Section 3.3 , such as finding indices I and J ,constructing the local submatrix, decomposing the local submatrix intoQR , and solving the upper triangular linear systems. For MPI, it provides asimple and easy-to-use parallel controlling capability on multicore CPUs,which dedicates one thread for controlling one GPU. Second, when asparsity pattern of the preconditioner is given, we use the thread-adaptiveallocation strategy to choose the optimized number of threads for eachcolumn of the preconditioner, and construct the decision tree to choosethe optimization kernel to calculate each one of components0https://doi.org/10.1016/j.tbench.2021.100007 2021年8月6日收到；2021年10月11日修订后接受；2021年10月20日在线发布 2772-4859/ © 2021 The Authors. Publishing servicesby Elsevier B.V. on behalf of KeAi Communications Co. Ltd. 本文是根据CC BY许可（http://creativecommons.org/licenses/by/4.0/）的开放获取文章。Y. Wang, W. Li and J. GaoBenchCouncil Transactions on Benchmarks, Standards and Evaluations 1 (2021) 100007min 𝐴𝑀 − 𝐸 2𝐹 .(1)min𝑚𝐴𝑚𝑘 − 𝑒𝑘22 ,𝑘 = 1, 2, … , 𝑛(2)20图1. 基于MPI的CPU-GPU混合并行计算模型。0预处理器的实验结果表明，HybridSPAI是有效的，并且优于CUSPARSE库中流行的不完整LU分解算法[30]、ViennaCL库中的静态SPAI预处理算法[24]和最新的GSPAI-Adaptive[29]。0本文的主要贡献总结如下。0∙ 基于MPI和CUDA的一般并行框架被提出。0本文的主要贡献总结如下。0∙ 提出了一种选择最佳线程数的策略0对于预处理器的每一列。0∙ 基于并行框架和提出的策略，0基于MPI和CUDA的优化SPAI预处理算法，称为HybridSPAI。在HybridSPAI中，通过并行计算找到索引，构造局部子矩阵，将局部子矩阵分解为QR，并解决上三角线性系统，并且通过决策树优化选择计算它们的核心。0本文的其余部分组织如下。第 2 节描述了SPAI预处理算法，第 3节给出了HybridSPAI的详细实现，第 4 节介绍了实验分析和评估，第 5节包含我们的结论并指出我们未来的研究方向。02. SPAI算法0SPAI过程的基本思想[22]描述如下：0使用稀疏矩阵 M 来近似 A 的逆，并且 M 由以下公式计算：0由于矩阵 M 的列是独立的，上述方程可以分解为以下 n 个独立的最小二乘问题0其中 � � 是单位矩阵的第 k 列， � � 代表矩阵 M 中的第 k列。有关SPAI的实现细节的描述，我们参考文献[27]。03. 在GPU上优化SPAI0我们提出了一种基于优化SPAI预处理算法的0在CPU-GPU平台上优化SPAI的混合并行计算模型，称为HybridSPAI。混合并行计算模型如图1所示。HybridSPAI的并行框架如图2所示，包括以下三个阶段：Pre-HybridSPAI 阶段，Compute-HybridSPAI 阶段和Post-HybridSPAI 阶段。0图2. HybridSPAI的并行框架。03.1. 基于MPI和CUDA的混合并行编程0必须为混合并行编程模型设计0GPU和CPU的体系结构以提高计算性能，并具有扩展到更多设备的特性。在我们提出的模型中，作为CUDA设备，GPU可以由多核CPU的每个线程控制，也可以由每个单独的CPU控制。此外，数据从主机内存传输到GPU设备内存，然后CPU通过调用核函数在GPU上启动计算过程。0MPI提供了一个简单方便的并行计算能力，0在多核CPU上的多线程的并行计算模型如图1所示，�1，�2，...，��3，��4，是存储在主机内存中的子矩阵，Thread是分配给CPU核心的多线程。0请注意，使用此模型时，计算矩阵将会被分成多个子矩阵，这些子矩阵对应于CPU的调用线0对应于CPU的调用线程数，计算矩阵将被分成多个子矩阵，因此这些子矩阵将分别分配给每个GPU来执行。30Y. Wang, W. Li and J. Gao BenchCouncil Transactions on Benchmarks, Standards and Evaluations 1 (2021) 1000070表1HybridSPAI中的数组。0数组大小类型0AData非零 Double AIndex非零 Integer APtr n Integer I ns × n1max Integer atmoic n Integer Jns × n2max Integer jPTR ns Integer �� ns × n2max Double �� ns × n1max × n2max Double R ns ×n1max × n2max Double iPTR ns Integer03.2. 预-HybridSPAI阶段0在本文中，我们预先总结了M的稀疏性0[25]中的主要方法。如果 M(i,j) 被认为是非零的，那么0|�(�, �)| > (1 − �)max�|�(�, �)|，0 � � � 1 (3)0满足以下条件时，�是用户定义的容差参数（主对角线始终包括在内）。0接下来，A使用压缩稀疏列（CSC）存储在主机内存中0列(CSC)存储格式，M也以列存储。本地子矩阵（� 1 �，� 2 �）的维度通常对不同的k（k = 1, 2,...,n）是不同的。为了简化内存中数据的访问并增加合并度，所有本地子矩阵的维度都统一定义为（� 1 ���，� 2 ���），其中 � 1 ��� = max � {� 1 �}，� 2 ��� = max � {� 2 �}。0最后，提出了线程自适应分配策略。对于任何0矩阵，预处理器每列的线程数 � 由以下公式计算：0� = min (2�, ��) (4)0s.t. 2 � −1 < � 2 ��� � 2 � . (5)0在Eqs.(4)中，��是固定的线程块大小。�个线程被分组成一个线程组，分配给计算M的第k列。在本文中，Eqs.(4)中的小写‘L’需要计算合适的�线程。请注意，在本文中，我们使用了一个线程块的1D数组来组织计算网格，并且也使用了一个线程的1D数组来组织线程块。03.3. 计算-HybridSPAI阶段0在计算-HybridSPAI阶段，每个GPU的分配0全局内存在表1中显示。基于消息接口的特性，MPI非常方便在CPU的多个线程之间分发和收集数据。实现计算 M 的以下步骤。0查找�和�：在所有块中，每个线程组块大小都是0用于找到 J 和 I 相同，并且每个线程组（warpSize个线程）被分配来找到 J 和 I的一个子集，这样就可以同时获得许多 J 和 I的子集。此外，还在每个线程组内部利用并行性。对于找到 J的内核，warp（线程组）内的线程并行读取 M的一列，并使用原子操作将它们存储到共享内存中。对于找到 I的内核，建立了一个决策树，对于任何给定的 � 2 ��� 和 � 1���，这个优化的内核都是有效的。图3显示了找到 I 的决策树的一个片段。当4 256时，cuFindI内核在全局内存上执行，而不是在共享内存上执行，这类似于cuFindIBySharedMemory内核。0构造本地子矩阵：使用上述获得的J和I0本地矩阵集��是通过具有共享内存或具有全局内存的内核计算的，具体取决于建立的决策树。图5显示了构造��的决策树的一部分。当4 < n2max ≤8时，根据不同的�1���，将选择具有sharedSize大小的共享内存的cuComputeTildeABySharedMemory内核，或具有全局内存的cuComputeTildeA内核。图6显示了cuComputeTildeABySharedMemory内核的主要过程。对于每个GPU上计算���的线程组，线程组中的所有线程0首先并行将��中的值读入共享内存��中，并���0通过在共享内存中加载在��中索引的列并并行匹配它们到��来构造。当n1max >256时，cuComputeTildeA内核在全局内存上执行，而不是在共享内存上执行，这类似于cuComputeTildeABySharedMemory内核。0将本地子矩阵分解为QR：线程组大小0将本地子矩阵分解为QR的步骤在所有块中相同。与上述两个步骤类似，再次使用构建的决策树40Y. Wang, W. Li and J. Gao BenchCouncil Transactions on Benchmarks, Standards and Evaluations 1 (2021) 1000070图5. 构造��的决策树的一部分。0图6. cuComputeTildeABySharedMemory内核的主要过程。0图7. 将本地子矩阵分解为QR的决策树的一部分。0分解本地子矩阵。图7显示了将本地子矩阵分解为QR的决策树段。当4 < n2max≤8时，根据不同的n1max，将选择具有sharedSize大小和sharedQ大小的共享内存的cuQRByQRSharedMemory内核，或具有sharedR大小的共享内存的cuQRByRSharedMemory内核。图8显示了cuQRByQRSharedMemory内核的主要过程。此外，每个线程组负责一个QR分解。有关其详细实现的描述，请参阅文献[25]。在一个线程组中，本地子矩阵，例如���，是0通过以下四个步骤将每次迭代�中的本地子矩阵分解为QR。0图8. cuQRByQRSharedMemory内核的主要过程。0图9. 解决上三角线性系统的决策树的一部分。0第一步，将��的第i列并行读入共享内存sQ中。第二步，线程并行计算��的第i行上三角矩阵并放入共享内存sR中。第三步，��的第i列和sQ同时进行归一化，并计算投影因子��和sR。第四步，使用共享内存sQ和sR并行更新��的所有列的值。当n1max >128时，利用共享内存sR执行cuQRByRSharedMemory内核，而不是共享内存sQ，这类似于cuQRByQRSharedMemory内核。0解上三角线性系统：本节中，通过使用共享内存并行求解上三角线性系统0计算 �� � 的子集 = � −1 � � � � �� � 通过并行求解上三角0线性系统。图9 显示了解上三角线性系统决策树的一个片段。对于任何给定的n2max 值，选择具有256大小的共享内存和 warpSize 线程组大小的cuSolverBySharedMemory。例如，当 4 < n2max ≤ 8时，选择具有256大小的共享内存和线程组大小为8的cuSolverBySharedMemory 核心。图10 显示了 cuSolverBySharedMemory核心的主要过程。对于每个线程组T，计算 �� 的步骤，例如 � � �，包括：(1)并行计算 � � � �� �0保存到共享内存 ��，并且 (2) 获得 �� � 的值venkat0110000/351.312imagesensor10000/521.036cfd277685.16716133.518apache258138.06111063.032t2em16613.1227682.338thermal240959.77125849.748G3_circuit10000/4752.53venkat010.5060.5010.2620.1510.1020.8060.7710.7610.7360.7151.3121.2721.0230.8870.817imagesensor0.2280.2270.1790.1030.1040.8080.7850.7670.7450.7131.0361.0120.9460.8480.817cfd21.1871.1910.6310.3560.2242.3312.2312.2942.1782.1013.5183.4222.9252.5342.325apache20.2260.2190.1260.1010.1332.8062.7612.7462.7342.8383.0322.9802.8722.8352.971t2em0.0750.0700.0600.0640.1032.2632.2532.2412.2312.2682.3382.3232.3012.2952.371thermal20.3320.3290.2010.1650.1649.4169.4439.3679.3699.1879.7489.7729.5689.5349.351G3_circuit0.1670.1560.1130.0940.1152.3632.3022.3212.3292.2902.5302.4582.4342.4232.405ulaCSRILU) [32], a static sparse approximate inverse preconditioningalgorithm in ViennaCL (denoted by SSPAI-VCL) [33], and a latest paral-lel SPAI preconditioning algorithm(denoted by GSPAI-Adaptive) [29].Table 5 demonstrate the comparison results on GTX 1080 Ti GPUs.For each matrix and the preconditioner, the first row is the computingtime of these four preconditioning algorithms, and the second rowand the third row are respectively the execution time and the numberof iterations of GPUPBICGSTAB while the residual error is less than1𝑒−7. Note that ‘‘/’’ represents the number of iterations for HybridSPAIexceeds 10,000, and all the other rows for each matrix will be denotedexcept that the third row is denoted by ‘‘> 10000’’. The minimum valueof50Y. Wang, W. Li 和 J. Gao BenchCouncil Transactions on Benchmarks, Standards and Evaluations 1 (2021) 1000070图10. cuSolverBySharedMemory 核心的主要过程。0图11. 组装 M 。0表2 测试矩阵描述。0名称 种类 行 非零元 平均 最大 最小0venkat01 CFD 序列 62,424 1,717,792 27.52 44 16 imagesensor 半导体器件 118,758 1,446,39612.18 21 2 cfd2 CFD 问题 123,440 3,085,406 25.00 30 8 apache2 结构 715,176 4,817,870 6.748 4 t2em 电磁学 921,632 4,590,832 4.98 5 1 thermal2 热力学 1,228,045 8,580,313 6.99 11 1G3_circuit 电路仿真 1,585,478 7,660,826 4.83 6 20通过使用共享内存并行求解上三角线性系统 � � �� � = xE。03.4. HybridSPAI 后阶段0后-HybridSPAI 阶段是将 M 以 CSC 存储格式组装0从多个 GPU 中收集数据格式。图11 说明了在每个 GPU 上组装 MPtr、MIndex和 MData 数组的过程。首先，利用 jPTR 组装 MPtr。其次，利用 �� 和 J 组装MIndex 和 MData。最后，根据 GPU 的设备 ID，将每个 GPU 上的 MData数组传输到 CPU 的相应线程。在 CPU 上，每个线程利用 MPI 的 MPI_Gatherv()函数并行地将 MData 收集到一个完整的数组中。04. 评估和分析0本节评估了 HybridSPAI 的性能。测试0表2 中的矩阵用于评估从佛罗里达大学稀疏矩阵收集中选择的 NVIDIA GTX 1080Ti GPU 的性能。源代码使用 CUDA 工具包 10.1 编译和执行。04.1. 有效性分析0对于每个测试矩阵，调用 GPUPBICGSTAB 解决 Ax = b0GTX 1080 Ti，其中所有 b 的值为1，生成的 M 用作预处理器。当残差误差小于1 �−7，或迭代次数超过10,000时停止。表3 显示结果，时间单位为秒（s）。0此外，我们使用 GTX 1080 Ti 来研究单个0GPU 和增加线程数量对 HybridSPAI 和 GPUPBICGSTAB 执行时间的影响。表4示0表3 GTX 1080 Ti 上 GPUBICGSTAB 的迭代次数和执行时间。0矩阵 GPUBICGSTAB GPUPBICGSTAB0迭代次数 执行时间 迭代次数 执行时间0表4 HybridSPAI 和 GPUPBICGSTAB 的执行时间。0矩阵 GPU 1 线程 2 线程 4 线程 8 线程0执行时间。对于每个矩阵和给定的线程数，第一行和第二行分别是HybridSPAI和GPUPBICGSTAB的计算时间，第三行是前两行的时间总和。GPUPBICGSTAB在残差误差小于1�−7时停止。每个矩阵的第二行和第三行的最小值都用红色字体标记。此外，我们观察到当计算预处理器的时间保持在单GPU上小于228毫秒时，增加GPU数量不能提供显著的加速。04.2. 性能比较0我们通过将HybridSPAI的性能与流行的进行比较0从表5中，我们观察到在GTX 1080 Ti上，总时间0对于任何矩阵，HybridSPAI和GPUPBICGSTAB的差值都是所有算法中最小的。这表明HybridSPAI优于CSRILU和SSPAI-VCL，并且优于GSPAI-Adaptive。Y. Wang, W. Li and J. GaoBenchCouncil Transactions on Benchmarks, Standards and Evaluations 1 (2021) 100007venkat011.83538.8560.5060.1021.5740.0360.8060.715114835353.42738.8921.3120.817imagesensor//0.2280.104//0.8080.71310000100005252//1.0360.817cfd2//1.1870.224//2.3312.101100001000016131613//3.5182.325apache23.38643.5320.2260.1016.7762.9952.8062.73447525031106110610.16246.5273.0322.835t2em19.884/0.0750.0642998.63/2.2632.231427100007687683018.514/2.3382.295thermal25.502/0.3320.16445.008/9.4169.1871619100002584258450.510/9.7489.351G3_circuit5.245/0.1670.11512.475/2.3632.2902571000047547517.720/2.5302.40560表5 GTX 1080 Ti上所有预处理算法和GPUPBICGSTAB的执行时间。0矩阵CSRILU SSPAI-VCL GSPAI-Adaptive HybridSPAI05. 结论0我们提出了一种高效的并行稀疏逆预处理0本文中我们提出了一种基于MPI和CUDA的多GPU上的并行稀疏逆条件数算法，称为HybridSPAI。在我们提出的HybridSPAI中，采用了一种通用的并行框架来优化多GPU上的静态SPAI，并建立了一棵决策树来选择计算它的最佳核。实验结果表明了我们提出的HybridSPAI的显著性能和高效性。0参考文献0[1] Y. Saad, M.H. Schultz, GMRES: 一种广义最小残差算法，用于0解非对称线性系统，SIAM J. Sci. Stat. Comput. 7 (1986) 856–869.0[2] H.A: Bi-CGSTAB: 一种快速且平滑收敛的BiCG变体，用于解决0非对称线性系统，SIAM J. Stat. Comput. 13 (2) (1992) 631.0[3] Y. Saad, 稀疏线性系统的迭代方法, 第二版, SIAM,0费城，宾夕法尼亚州，2003年。0[4] J. Gao, R. Liang, J. Wang, 具有修改的不完全Cholesky预处理的共轭梯度算法研究0GPU上的修改不完全Cholesky预处理共轭梯度算法，J. Parallel Distr. Com. 74 (2) (2014)2088–2098。0[5] S.C. Rennich, D. Stosic, T.A. Davis, 加速稀疏Cholesky分解0GPU上的预处理Krylov求解器，Parallel Comput. 59 (2016)0[6] H. Anzt, M. Gates, J. Dongarra, M. Kreutzer, G. Wellein, M. Kohler, GPU上的预处理Krylov求解器0GPU上的预处理Krylov求解器，Parallel Comput. 68 (2017) 32–44。0[7] L.Y. Kolotilina, A.Y. Yeremin, 因式分解的稀疏近似逆预处理0因式分解的稀疏近似逆预处理I.理论，SIAM J. Matrix Anal. Appl. 14 (1) (1993) 45–58。0[8] M. Benzi, C.D. Meyer, M. Tuma, 用于稀疏线性系统的稀疏近似逆预处理器0共轭梯度法，SIAM J. Sci. Comput. 17 (5) (1996) 1135–1149。0[9] M. Ferronato, C. Janna, G. Pini, 用于非对称线性系统的广义块FSAI预处理器0非对称线性系统的多级方法，J. Comput. Appl. Math. 256 (2014) 230–241。0[10] V.A.P. Magri, A. Franceschini, M. Ferronato, C. Janna, 多级方法0FSAI预处理，Numer. Linear Algebr. (2018) e2183，http://dx.doi.org/10.1002/nla.2183。0[11] M.M. Dehnavi, D.M. Fern´andez, J.L. Gaudiot, D.D. Giannacopoulos, 并行稀疏近似逆预处理0GPU上的稀疏近似逆预处理，IEEE T. Parall. Distr. 24 (9) (2013) 1852–1861。0[12] Z. Jia, B. Zhu, 一种功率稀疏近似逆预处理程序0大型稀疏线性系统的稀疏近似逆预处理，Numer. Linear Algebr. 16 (4) (2009) 259–299。0[13] J.D.F. Cosgrove, J.C. Diaz, A. Griewank, 近似逆预处理0稀疏线性系统的近似逆预处理，Int. J. Comput. Math. 44 (1–2) (1992) 91–110。0[14] M. Grote, T. Huckle, 使用稀疏近似逆进行并行预处理0SIAM J. Sci. Comput. 18 (3) (1997) 838–853。0[15] E. Chow, Y. Saad, 通过稀疏-稀疏乘法得到近似逆预处理器0SIAM J. Sci. Comput. 19 (3) (1998) 995–1023。0[16] E. Chow, 并行稀疏近似逆的先验稀疏模式0SIAM J. Sci. Comput. 21 (5) (2000) 1804–1822。0[17] E. Chow, A. Patel, 细粒度并行不完全LU分解0Sci. Comput. 37 (2) (2015) C169–C193。0[18] J. Gao, Z. Li, R. Liang, G. He, 自适应优化l1-最小化求解器0GPU上的自适应优化l1-最小化求解器，Int. J. Parallel Program. 45 (3) (2017)0[19] K. Li, W. Yang, K. Li, 用于GPU的混合并行求解算法0对角线性方程组的GPU上的混合并行求解算法，IEEE Trans. Parallel Distrib. Syst. 27 (10) (2016)2795–2808。0[20] J. Gao, Y. Zhou, G. He, Y. Xia, 用于多GPU的并行优化模型0预处理共轭梯度算法，Para
我的内容管理
收起
我的收益 登录查看自己的收益
我的积分
登录查看自己的积分
我的C币
登录后查看C币余额
我的收藏 
我的下载 
下载帮助 
