用于大规模并行有限元模拟的可扩展混合全FETI方法

0用于大规模并行有限元方法仿真的可扩展混合全FETI方法0林科浩1，2，†，周春宝2，3，†，曾彦1，聂宁明2，3，王珏2，3，�，李世刚4，冯扬德2，3，王彦刚2，3，姚克汉1，2，姚铁垂2，3，张吉林1，万健101 杭州电子科技大学，中国杭州 2  中国科学院计算机网络信息中心，中国北京 3中国科学院大学，中国北京 4  北京邮电大学计算机科学学院，中国北京电子邮件：wangjue@sccas.cn0摘要混合全有限元撕裂和互连（HTFETI）方法在解决大规模和复杂工程问题中发挥着重要作用。该方法需要处理大量的矩阵-向量乘法。直接调用供应商优化的通用矩阵-向量乘法（gemv）库在GPU上性能较低，因为它未考虑HTFETI中不同矩阵大小（即不同的行和列大小）的优化。此外，最先进的图分割方法无法保证HTFETI的负载平衡，因为矩阵大小取决于子域边界的长度。为解决上述问题，我们首先将gemv移植到多流水线方案，并在GPU上开发了新的批处理内核函数，其吞吐量提高了15%~30%，平均GFLOP提高了37%。我们还提出了一种基于图重新分割和工作窃取的多粒度负载平衡方案，负载不平衡比例从1.5降至1.05~1.09。我们成功地将可扩展的HTFETI方法应用于模拟中国实验快堆（CEFR）的整芯装配进行稳态分析，弱可扩展性和强可扩展性的效率分别达到78%和72%，在12288个GPU上。据我们所知，这是HTFETI首次在大规模和高保真度的整芯装配模拟中使用。0CCS概念:  • 计算方法学  → 大规模并行算法。0关键词:  混合全FETI，GEMV，负载平衡0PPoPP ’23，2023年2月25日至3月1日，加拿大蒙特利尔©2023年计算机协会 ACM ISBN979-8-4007-0015-6/23/02。https://doi.org/10.1145/3572848.357751701 引言0随着机械工程、土木工程、生物力学和能源行业等复杂工程仿真对精度和规模的需求日益增加，大规模有限元方法（FEM）至关重要。直接解线性系统的方法更加稳健，但在大系统上可能存在内存使用问题[26]。有限元撕裂和互连（FETI）方法[13]是一种基于Krylov空间的迭代方法，将问题域划分为子域，并通过拉格朗日乘子（LM）实现子域之间的连续性。在全FETI（TFETI）方法[9]中，Dirichlet边界条件也通过LM来实施，这在数值解中具有优势。TFETI仍然受到粗问题规模的限制[12，32]，不能分解更多的子域。通过将少量相邻子域聚合成簇，从而产生更小的粗问题，混合全FETI（HTFETI）方法[23]使得原始问题的并行化能够扩展到成千上万个核心，因为需要更少的内存。HTFETI具有良好的可扩展性，因此更受欢迎。同时，GPU提供强大的计算能力，并且超级计算机往往使用GPU和CPU。HTFETI优化异构系统对于复杂工程仿真至关重要。计算开销和负载不平衡[5，14]使优化成为一项高回报的任务，这激发了我们的工作。在HTFETI中，所有的gemv操作占据了90%的运行时间。一组中等大小的方阵gemv占据了所有gemv时间的90%。行�较小且列�较大的大矩阵的操作占据了剩余的所有gemv时间。本文使用局部舒尔补（LSC）[30]方法将稀疏矩阵-向量乘法转换为密集矩阵-向量乘法。因此，主要的计算开销是密集矩阵-向量乘法，这是一个典型的内存访问密集型问题。对于异构计算系统，内存访问瓶颈特别突出。指令级并行ism（ILP）和线程级并行ism（TLP）可用于实现高0†  合作第一作者。�通信作者。0本作品采用知识共享署名国际4.0许可协议进行许可。=0PPoPP ’23，2023年2月25日至3月1日，加拿大蒙特利尔，林科浩，周春宝，曾彦等。0GPU上的内存吞吐量。对于中等大小的方阵gemv在AMDGNUGPU上，我们使用四个流，每个流独立完成一个稠密矩阵-向量乘法。gemv例程在操作大矩阵时性能较差，因为其线程块较少。许多大矩阵更适合批处理gemv，并且在矩阵之间分配线程块可以有效提高并行性。传统的有限元方法使用分割网格内部元素来组装矩阵，而FETI使用分割网格边界元素来组装矩阵。FEM的负载平衡基于图分割。虽然这种方法可以均匀地划分网格内部元素，但无法保证FETI中某些复杂模型中网格边界元素的平衡。在对中国实验快堆（CEFR）整芯装配进行稳态分析时，我们发现最大和最小运行时间差可以达到6倍，最大和平均运行时间差可以达到1.5倍。在进行CG迭代之前，需要调整计算量。我们的贡献如下：0• 我们将中等大小的方阵gemv移植到AMDGPU上，使用多流水线调度，并获得最高的内存吞吐量。对于许多大矩阵gemv，我们开发了可变大小的批处理内核函数，以提高GPU利用率。0•我们提出了一种多粒度的负载平衡方案。首先，根据上一次图分割后网格边界元素的统计信息，通过图重新分割来执行粗粒度的负载平衡。其次，通过使用gemv性能评估来进行繁忙进程和空闲进程之间的工作窃取来执行细粒度的负载平衡。0• 在AMD Radeon ProMI60集群上的评估结果显示，我们的优化方法在HTFETI的大规模仿真中实现了良好的可扩展性，并显著提高了负0•通过在3072个计算节点上使用12288个GPU应用优化的HTFETI，我们首次完成了100亿网格单元的CEFR整芯装配稳态模拟。在反应堆中，整芯装配的设计和寿命延长非常重要。我们使用真实环境的参数模拟整芯装配在热和压力下的变形，为整芯装配更换提供重要参考。0我们提出的多流水线调度和批处理内核函数，多粒度负载平衡方案为真实世界的大规模并行应用程序提供了参考0在配备有GPU加速器的超级计算机上进行。GNN[24]也可以应用所提出的负载平衡方法。02 背景在本节中，我们简要介绍HTFETI，预条件共轭投影梯度（PCPG）方法和负载平衡优化的动机。02.1HTFETI在结构力学中，FEM定义了全局刚度矩阵K，节点位移u和全局节点力f来描述系统Ku =f。传统的FEM通过直接方法解决线性系统。目的是找到原始未知数u的值。FETI方法是一种用于解决大规模问题的非重叠域分解方法[11, 13]。引入LM�以保持相邻子域之间的原始变量的相互连续性。B1定义了子域中边界节点与内部节点之间的映射关系。第i个子域的平衡方程式如式（1）所示。0K�u� = f� −BT1,��（1）TFETI将Dirichlet边界条件引入B1矩阵，通过�强制满足B1u =c。HTFETI将多个子域聚合成簇。选择子域的角节点来增加簇中解的连续性，并引入��。Bc定义了子域中角节点与内部节点之间的映射关系。我们假设系统中有四个子域，每两个子域属于一个簇。通过矩阵置换，系统可以从（2）变为（3）。0K1 O O O BTc,1 O BT1,1 O K2 O OBTc,2 O BT1,2 O O K3 O O BTc,3BT1,3 O O O K4 O BTc,4 BT1,4 Bc,1Bc,2 O O O O O O O Bc,3 Bc,4 O O OB1,1 B1,2 B1,3 B1,4 O O O0�����0������0u1u2u3u4�c,�c,�0���������0����������0f1f2f3f4ooc0��������������(2)0���0����0�������0����������0��������������(3)用符号表示，�K结合了K�和Bc,�，�u结合了u�和�c,�，�B结合了B1,�，�f结合了f�。我们得到全局方程式（4）。�R是�K的零空间的核，�K+是�K的Moore-Penrose逆矩阵，�是振幅。使用新符号�F = �B�K+�BT，�d = �B�K+�f − �c，�G = −�RT�BT，�e =−�RT�fT，我们得到要解决的新Schur补系统（5）。0用于大规模并行FEM模拟的可扩展混合总FETI方法PPoPP’23，2023年2月25日至3月1日，加拿大蒙特利尔0� �F �GT0�u = �K + � �f − �BT�� + �R�� (4)0� �� �0�GO0� (50�=��0Algorithm 1: PCPG02.2 PCPG系统（5）可以通过Krylov子空间迭代方法求解。PCPG [22]将在以下部分中描述，如算法1所示。01 r0 ← �d − �F�0;0Data: �0←�GT(�G��GT)−1e, � > 0, imax > 003 Project: w� ←05 Project: y� ← �Pz�;01w�;07 p1 ← y1;06 if � = 1 then09 � ← (y��, w�)/(y��−1,08 else011 if ||w�|| < � then010 p�←y�−1+�p�−1;013 � ← (y��, w�)/(p��, �Fp�);015 r� ← r�−1 − ��Fp�;0和预处理矩阵F−1给出。投影0在开始时，投影矩阵�P = I − �G��G��G��−1�G�0每次迭代的每个部分，簇LM �clu（绿色）由子域LM �sub（红色）组成，只有簇的边界LM0矩阵P包含粗问题G�G−1∈Rncp×ncp。在线性弹性问题中，其维度由ncp = d∙n给出。对于2D中的元素，d =3，它包括x，y方向的位移和x-y平面的旋转，3D中d =6。在TFETI中，n是子域的数量。在HTFETI中，n是簇的数量。HTFETI中的粗问题规模比TFETI中小得多。算法1由多个MPI进程并行实现。一次迭代包括以下部分：投影（第3行，第5行），预处理（第4行），点积（第9行，第13行）和残差范数（第11行）。在图1中，假设每个MPI进程对应一个簇，当完成每个迭代的每个部分时，簇LM �clu（绿色）由子域LM�sub（红色）组成，只有簇的边界LM�comm（蓝色）需要进行图案通信更新。其他通信在点积和残差范数的全局约简中发生，这意味着通信量较小，计算占据大部分运行时间。0为了计算F−1w�（第4行）和�Fp�（第13行），簇收集每个子域的矩阵-向量乘法结果。特别地，对于HTFETI，�Fp被划分为两部分[25]，如式（6）所示。�: �是属于第i个簇的子域的集合。0图1. 红色：子域LM �sub，绿色：簇LM �clu，蓝色：簇边界LM �comm0式（6）中的左部分是每个子域的中等大小的方阵B1K +B1。行数等于双重未知数�的数量，分布在1,000到8,000之间。式（6）中的右部分操作胖矩阵[31]。R的行数为6，Bc的行数是6的倍数且不超过48。K+的列数等于原始未知数u的数量，分布在4,000到48,000之间0Fp = B1,�: � K + �: � BT1,�: � p + B1,�: � R�: � � − K + �: � B�c,�: � � (6)02.3 动机FEM使用图分割库进行域分解[4]。图分割方法使得每个分区的内部元素大小均衡。它使得原始未知数u均匀分布。然而，这不能保证分区边界元素的大小平衡，对于HTFETI来说，这使得每个簇处理的双重未知数�的大小相当不同。子域矩阵B1K +B1的大小与�有关。每个簇中矩阵-向量乘法的工作负载是不同的。像飞机，汽车[17]和燃料棒这样的复杂模型，通常具有不均匀的密度，特别是某些区域是实体的，而其他区域是空心的。图分割后，分区边界大小不平衡。实体部分表示网格(a)(b)meshclusterssubdomainsmesh elementsdomains connected0PPoPP’23，2023年2月25日至3月1日，加拿大蒙特利尔，Kehao Lin，Chunbao Zhou，Yan Zeng等人0（a0网格簇子域0网格元素0连接的域0图2.a）两级分解：低空心杆网格被分成两个簇，簇网格被分成子域b）子域及其相连的子域0元素紧密相连。空心部分表示网格元素与空间空元素（如空气）接触，典型的空心部分例如门和螺旋桨。我们以CEFR整体堆芯装配的简化模型为例，该模型由一个具有实体顶部和空心底部的杆组成，以说明在FEM中图分割后边界大小不平衡的原因。图分割运行两次。它首先将模型的网格元素分为具有相同元素数量的簇，然后将簇网格元素内部分为多个子域，每个子域具有相同数量的元素。在图2（a）中，杆被分成两个簇。顶部簇的元素数量与底部簇相同。接下来，每个簇被分成子域。实体部分的元素与其他元素有更多的接触机会，因此子域的接触邻居数量更多且接触面积更大。空心部分的元素与其他元素的接触机会较少，因此子域的接触邻居数量较少且接触面积较小。如图2（b）所示，实体部分的子域有六个邻居，接触面积为完全接触表面。空心部分的子域有四个邻居，接触面积为窄接触表面。我们对整个堆芯装配运行了两级图分割。每个簇包含相同数量的子域。运行迭代时，每个簇的独立矩阵-向量乘法操作时间存在较大差异。这是一个重要的负载不平衡问题，可以进行优化。03 加速0对于HTFETI需要多次执行的矩阵-向量乘法，我们可以调用BLAS2操作gemv，它执行 �  =  ���  +  �� 。有许多广泛使用的BLAS库，例如cublas、rocblas[3]和magma。在我们的平台上，每个计算节点有4个AMD GPU，因此0我们关注如何调用rocblas例程或自定义新的内核函数。AMDGCN架构[1]的信息如下所示。一个GPU包含64个计算单元（CU），CU共享设备全局内存。每个CU包含4个SIMD单元，每个SIMD包含16个矢量算术逻辑单元（VALU）。工作负载管理器从命令队列中弹出内核并根据负载将其分配给CU运行。对于内核函数，线程组织在一个称为工作组的组中。一个内核可以包含多个工作组。工作组分布在许多CU上，每个工作组始终在单个CU上运行，直到结束。波前是CU上并行操作的64个宽度矢量线程。一个工作组可以包含多个波前，因此工作组中的线程数最好组织为64的倍数。SIMD单元需要4个周期来执行所有波前的线程操作，以锁定步进方式执行。SIMD执行多个锁步以完成一个波前。rocblas在方阵中的表现良好，但在胖矩阵中表现不佳。我们演示了如何为方阵安排rocblas例程调用（第3.1节），以及如何为胖矩阵开发批处理内核函数（第3.2节）。03.1多流水线调度作为参考方法，我们使用标准的gemv例程（rocblas_dgemv例程）。对于足够大的矩阵，rocblas_dgemv内核从16个工作组开始。内核可以放置在称为流的工作队列中。当资源足够时，GPU调度程序将自动调度流中的内核。我们逐渐增加流的数量，并测试rocblas_dgemv的GFLOP。对于AMD GNUGPU，从1到4个流的GFLOP显著增加，并且很难在4个流以上获得改进。流配置与内核工作组和CU的数量有关。在工作组太少时，GPU无法充分利用。最多有64个CU同时工作。将�乘法放入多个流中还有另一个好处[2]，即可以重叠设备到主机的复制（D2H）、主机到设备的复制（H2D）和流之间的内核。通过分配固定内存，可以使用DMA技术在GPU和CPU之间复制数据。这样可以避免从分页内存复制到临时固定内存的步骤（图3(a)），并且可以在图3(b)中异步发送CPU命令。多流水线调度在ILP级别上提高了性能。总内存吞吐量可以提高15%至30%。03.2 批处理内核函数0对于胖矩阵，rocblas_dgemv例程为一个矩阵使用  � /16个工作组，这意味着GPU利用率最低。在这种情况下，即使为并发配置了更多的流，也没有足够的工作量来覆盖启动开销。对于许多小矩阵，批处理线性代数例程针对更好的性能进行了优化，而rocblasDRAMPinned MemPaged MemHostDevicePCIeCPUStream1Stream2CPU H2D commandCPU kernel commandGPU H2D executionGPU kernel execution(a)(b)Wavefront 0Wavefront 1Wavefront 2Wavefront 3Wavefront 0Wavefront 1Wavefront 2Wavefront 3Workgroup 0Workgroup 1Wavefront………64 threads64 threads64 threads(a)(b)0用于大规模并行有限元模拟的可扩展混合总FETI方法 PPoPP’23，2023年2月25日至3月1日，加拿大蒙特利尔0主机0设备0CPU0流10流20CPU H2D命令 CPU内核命令0GPU H2D执行 GPU内核执行0(a) (b)0图3. a) 内存数据传输 b)流并发，CPU异步调用和GPU在复制和内核之间重叠执行0提供了两个例程：rocblas_dgemv_batched和rocblas_dgemv_strided_batched。批处理内核函数在不同矩阵之间分配线程，并在不同矩阵上使用更多的工作组。rocblas批处理例程要求每个矩阵具有相同的大小。如果我们对不同大小的矩阵使用rocblas批处理例程，则所有矩阵都需要调整为最大矩阵的大小，并填充零元素。此外，我们发现对于非常小的矩阵，即行数�为6的矩阵，rocblas批处理例程的性能表现不佳。我们开发了一个批处理内核函数，易于使用，并提高了应用程序的性能。0- 主机（CPU）端代码01. dim3  block (64, 4); 2. dim3  grid (m / 4,nmatrices); 3. gpukernel <<>>(...);0- 设备（GPU）端代码01. __global__  void  gpukernel02. ( double *  mat ,  double *  vec ,  double *  out , ...) 3. {       4. int  cur_mat  = blockIdx.y;5. int  cur_row  = blockIdx.x * blockDim.y + threadIdx.y; 6. if  (cur_row < mat_h[cur_mat]){ 7. double  res  = 0.; 8. int  laneId  = threadIdx.x % 64; 9. int  w  = mat_w[cur_mat]; 10. int kIteration  = (w - 1) / 64 / 2 + 1; 11. mat = &mat[mat_head[cur_mat] + cur_row * w]; 12.vec = &mat[mat_nstrtt[cur_mat]]; 13. for  ( int  i  = 0; i < kIteration; ++i) { 14. intcur_col_vec  = i * 64+ laneId; 15. double2  current_val  = reinterpret_cast(mat)[cur_col_vec]; 16. double2  current_x  = reinterpret_cast(vec)[cur_col_vec]; 17. res += current_val.x * current_x.x; 18. res += current_val.y *current_x.y; 19. } 20. for  ( int  offset  = 64 / 2; offset > 0; offset /= 2) { 21. res +=__shfl_down(res, offset); 22. } 23. if  (laneId == 0) out[mat_nstrtl[cur_mat] + cur_row] =res; 24. }0图4. GPU内核的伪代码0(a)0(b)0图5. a) 一个矩阵中的工作组 b) 一个波前中的线程组合0图4显示了新GPU内核的伪代码。注意，有 ���������个矩阵，我们为内核使用( � /4, ���������)个工作组。线程根据blockIdx.y判断当前操作的是哪个矩阵。对于不同的矩阵大小，根据偏移数组和大小数组确定工作数据和工作范围。这种安排是自然的，并提高了内存使用效率。图5支持一个矩阵的计算图像。一个工作组可以组织成 �个波前，rocblas_dgemv_batched将16个波前设置为一个工作组。但对于行数很小的矩阵，波前的数量是提前的，性能会下降。在我们的内核函数中，我们在图5(a)中为一个工作组放置了4个波前，为了更好的调度，工作组的数量也增加了。每个波前负责将矩阵  �  的一行与向量  �相乘，然后将其缩减为向量  �的一个列元素。在图5(b)中，波前中的所有线程连续访问64个double。这种安排满足了组合内存访问的条件，并确保从全局内存到寄存器的内存访问效率[6]。此外，我们执行矢量化内存访问，并尽可能使用像double2这样的128位内存访问指令，这带来了轻微的性能提升。04 负载均衡0负载平衡一直是有限元方法中的一个重要问题。矩阵-向量乘法是HTFETI中最耗时的部分。在CEFR整个堆芯装配体上进行模拟时，我们发现最大运行时间与最小运行时间之间的差异可以达到数倍之多。在第4.1节中，我们提供了粗粒度的负载平衡。考虑到分区的边界，我们再次运行图分割以实现负载平衡。在第4.2节中，我们提供了细粒度的负载0PPoPP ’23，2023年2月25日-3月1日，加拿大蒙特利尔，QC，加拿大 林科豪，周春宝，曾岩等0负载平衡。繁忙的进程将矩阵发送给空闲的进程进行工作窃取。04.1基于边界感知图重新分区的粗粒度负载平衡0图分区工具有METIS [19]、SCOTCH [28]和KaHIP[34]等，用于网格划分和静态负载平衡。现代图分区的目标是：减少跨分区的通信并生成平衡的分区[33]。它们被抽象为每个分区中顶点的数量平衡和最小边界切割。HTFETI过程之前的图分区基于顶点平衡，而由HTFETI组装的矩阵的大小可能不平衡。粗粒度负载平衡侧重于图分区结果，并找到实体部分和空心部分之间的差异。图分区的结果将告诉我们每个元素分配给哪个集群，以及分配给哪个子域。我们可以计算子域边界上的元素数量。边界元素的数量与每个集群的矩阵-向量乘法（即HTFETI的�）相关。图分区可以通过给图节点添加权重来影响输入。它使得每个分区的权重尽可能均衡，并完成元素数量的调整。图中有两个集群如图2（a）所示。如果我们给实体部分的元素分配高权重，并重新运行图分区，一些实体部分的元素可以被分配到空心部分，这反过来会影响图分区后的子域边界结果。粗粒度负载平衡的详细信息在算法2中描述。元素权重的确定是从图分区的统计边界结果中获得的。在开始时，进行两级图分区（第5行，第9行），我们计算每个子域分区上边界元素的数量，并得到每个集群中边界元素的总数量。我们为每个集群中的所有元素分配权重，该权重与边界元素的总量相关。如果集群中边界元素的总量很大，则元素位于实体部分，相应的权重较高。下一轮图分区将平衡划分的权重，并获得不同的划分结果。原先较重的集群减少了边界元素的数量和子域的数量，从而减少了边界元素的数量。该方案通过重新分区图来调整每个集群分区中的边界元素数量，并影响LM的数量。如果我们想要更多的负载平衡，需要运行时信息，因此这个方案被称为粗粒度负载平衡。但是，这个方案的优点如下：1）在组装方程之前可以调整每个集群的工作负载。0算法2：粗粒度负载平衡0输入：元素图：G，元素大小：esize，集群数：cnum，默认集群域数：dnum，平衡阈值：�，最大迭代次数：I0输出：集群元素集合：P，集群域数集合：dnums01 weighs ← (1, ...);02 dnums ← (dnum, ...);03 � ← 0;04 重复05 P ← graphPartition(G, weighs, cnum);06 wmax ← 0, wavg ← 0;07 for � = 1 to cnum do08 G� ← getSubGraph(G, P, c);09 P� ← graphPartition(G�, �, dnums(�));010 w ← getBoundaryStatistics(G�, P�);011 for � in P(�) do012 weighs(�) ← w / dnums(�);013 if wmax < w then014 wmax ← w;015 wavg ← wavg + w;016 wavg ← wavg / cnum;017 dnums(�) ← |P| × cnum / esize × dnum018 � ← � + 1019 直到 wmax / wavg < � 或 � < I020 返回 (P, dnums)0与解决阶段相比，图分区阶段运行非常快。2）我们可以提前使用低精度的网格而不是高精度的网格来提取模型密度特征。04.2基于工作窃取方法的细粒度负载平衡在矩阵组装完成后进行细粒度负载平衡。在这个阶段，确定每个集群中的子域数和对应的矩阵如图6（a）所示。在不同的硬件上，不同的数学库对不同大小的矩阵有不同的性能。我们使用rocblas库和在第3节中实现的批处理核函数，在当前集群处理的矩阵上进行矩阵-向量乘法的评估。然后我们将执行工作窃取。负载高于平均值的集群被视为被帮助者，负载低于平均值的集群被视为帮助者。高负载的集群将矩阵发送到低负载的其他集群中，如图6（b）所示。细粒度负载平衡的详细信息在算法3中描述。细粒度中的"performanceEvaluate"cluster 0cluster 1cluster 0 as helpee(a)0面向大规模并行FEM仿真的可扩展混合总FETI方法 PPoPP ’23，2023年2月25日-3月1日，加拿大蒙特利尔，QC，加拿大0算法3：细粒度负载平衡0输入：集群LM大小信息集合：D，集群数：cnum，硬件配置集合：T0输出：被帮助者集合：H0，帮助者集合：H1，被帮助者集合：P01 wavg ← 0;02 w ← (0, ...);03 for � = 1 to cnum do04 w(c) ← performanceEvaluate(D(c), T(c));05 wavg ← wavg + w(c);06 wavg ← wavg / cnum;07 perms ← sort_des(w);08 for � in perms do09 if w(�) > wavg then010 H0 ← H0 ∪ {�};011 else012 H1 ← H1 ∪ {�};013 h0 ← 1, h1 ←014 d0 ←015 while h0 < h1 do016 c0 ← perms(h0), c1 ← perms(h1);017 while w(c0) > wavg and w(c1) < wavg do018 D(c0) ← D(c0)\{d0};019 D(c1) ← D(c1) ∪ {d0};020 w(c0) ← performanceEvaluate(D(c0), T(c0));021 w(c1) ← performanceEvaluate(D(c1), T(c1));022 P(h0) ← P(h0) ∪ {d0};023 d0 ← d0 + 1;024 if w(c0) ≤ wavg then025 h0 ← h0 + 1;026 d0 ← 1;027 if w(c1) ≥ wavg then028 h1 ← h1 - 1;029 返回 (H0, H1, P)0集群0作为被帮助者 集群1作为帮助者0(a) (b)0图6. a) 每个子域对应一个矩阵和一个�的向量。 b)用于负载平衡的集群之间的工作窃取。0算法首先执行现有的dgemv任务，并记录时间。首先获取每个子域的矩阵大小，然后执行"performanceEvaluate"来获取dgemv的运行时间。我们将每个进程的dgemv运行时间作为负载平衡的度量标准。我们收集集群之间的所有运行时间，计算平均运行时间，然后根据运行时间对集群进行排序。我们优先考虑与平均值差异较大的帮助者和被帮助者之间的工作窃取。当帮助者（被帮助者）的运行时间下降（上升）到平均值时，我们切换到排序中的下一个帮助者（被帮助者）。直到所有集群达到平均值，工作窃取方法就会调整负载。这是负载平衡的细粒度调整。05 性能0实验是在通过Infiniband网络互连的先进计算平台上进行的。每个计算节点配备有一个32核心2.50GHz x86CPU，并由4个AMD Radeon Pro MI60GPU加速，通过PCIe进行互连。GPU基于GNU架构，由64个CU组成，共4096个核心。05.1 加速评估在实验中，我们首先在一个GPU上测试不同的rocblas调用方案的GFLOPs。我们使用多流dgemv（流程最多为5）和批处理dgemv，并测量内核函数的时间。GFLOPs的方程式表示为（7）。所有矩阵都是方阵，行数 �  等于列数 � ，矩阵数量 �大于16。对于gemv，乘法和加法操作数 f为2，目标向量维度 �  = 1。0GFLOPs  = �  ×  01 � 9 ×  � (7)0图7显示了各种方案的GFLOPs。结果显示，批处理的dgemv对于增长的矩阵尺寸具有最快的启动速度。例程rocblas_dgemv需要更多的流以实现全面加速，而4个流已足够。然而，在应用中，还需要考虑数据传输。多流可以利用流水线重叠。而矩阵  �  在GPU上提前固定，每次迭代从CPU输入向量  �到GPU，并输出向量  �从GPU到CPU。我们使用方程（8）评估内存吞吐量。对于双精度，字节大小  �  为8，单位为MB/s。0T  = �  ×  01024 2 ×  � (8)0对于HTFETI，考虑到内存能力，划分的子域数量既不能太多也不能太少。矩阵K的大小（子域自由度大小）从4000到48000分布。矩阵B1KBT1的大小从1000到8000分布。在实验中，矩阵大小从1024变化到8192。矩阵�的数量逐渐增加，可以分配的最大矩阵数量受到限制。方案的吞吐量如图8所示。与批处理方案相比，流程方案在矩阵尺寸为8192时表现更好。对于两种方案，随着矩阵尺寸的增加，吞吐量先增加后减少，矩阵尺寸约为2048时性能最佳。与批处理方案相比，流程方案带来了15% ~30%的吞吐量改进。大矩阵具有以下特点。矩阵尺寸的范围是 �  = [6, 48]，步长为6， �  = [4096,32768]，步长为4096。我们将矩阵数量� = 256设置为批处理内核函数实验。在图9中，当  �  < 24或  �>16348时，新的批处理内核函数在大矩阵上表现良好。在表1中，当  �  =6时，该内核可以实现高性能，平均GFLOPs改进为37%。0641282565121024204840968192Matrix Size1 stream2 stream3 stream4 stream5 streambatched010040060081632641282565121024 1536ThB/s]Batch Sizestreams 819240968192122881638420480245762867232768643.247.292.851.192.7 315.9 393.7 291.71213.86.921.82.821.4 113.0 123.3 98.41811.34.218.09.116.177.084.965.2240.2‐4.04.0‐6.31.441.144.932.5304.2‐0.14.5‐5.03.031.133.022.7364.60.03.6‐3.12.116.719.112.2427.95.58.11.78.116.717.511.5484.30.23.5‐2.11.917.218.914.6MN-60100percentage performance improvements0PPoPP ’23, 2023年2月25日至3月1日，加拿大蒙特利尔，Kehao Lin，Chunbao Zhou，Yan Zeng等人。0Gflops0矩阵尺寸01个流02个流程03个流04个流程05个流0批处理0图7.  不同方阵尺寸的DGEMV方案的GFLOPs0吞吐量[MB/s]0批处理 102401024个流0批处理 204802048个流0批处理 409604096个流0批处理 819208192个流0图8.  批处理和流程DGEMV方案的吞吐量比较024 0.2 ‐4.0 4.0 ‐6.3 1.4 41.1 44.9 32.5030 4.2 ‐0.1 4.5 ‐5.0 3.0 31.1 33.0 22.7036 4.6 0.0 3.6 ‐3.1 2.1 16.7 19.1 12.2048 4.3 0.2 3.5 ‐2.1 1.9 17.2 18.9 14.60M0性能改进百分比0图9.  批处理内核函数改进相对于rocBLAS的改进热图0表1.  DGEMV GFLOPs的格式为“kernel|rocBLAS”05.2 模拟评估0负载不平衡表现为集群的最大运行时间与平均运行时间之间的差异，我们定义最大运行时间差异与平均运行时间的比值为负载不平衡比例。负载不平衡现象在复杂模型中非常明显，我们对整个堆芯装配的4800万个网格单元进行负载平衡实验。在每个集群中有256个子域的情况下，我们使用16个计算节点并逐渐增加到64个计算节点。在图10中，粗粒度负载平衡可以将负载不平衡比例从1.5降低到1.15�1.20，而细粒度负载平衡可以进一步将负载不平衡比例降低到1.05�1.09。我们进行消融实验来观察GPUdgemv和负载平衡优化对模拟处理时间的影响。求解器运行时间阶段执行PCPG。PCPG迭代次数为200。我们补充了两组比较实验：首先，对在不同硬件平台上运行dgemv的实验比较 -在CPU上和在GPU上。其次，比较不同的负载平衡方法 -使用METIS图分区和使用多粒度负载平衡方法。在图11中，记录了模拟的处理时间。使用负载平衡时，加速比为1.09148.63%152.41%150.76%117.23%115.36%120.15%109.25%107.20%105.23%12163264126.969.943.0 115.9 63.1 38.0 51.3 36.8 28.9 47.8 33.9 26.1 0204060801001201401601632640一种可扩展的高性能并行有限元方法用于海量并行FEM模拟 PPoPP ’23, 2023年2月25日至3月1日，加拿大蒙特利尔，QC，加拿大0负载平衡比例0计算节点0理想的图分区0粗粒度方案 细粒度方案0图10. 整个堆芯装配模拟的负载平衡改进。0到1.13。使用GPUdgemv优化时，加速比为1.48到2.47。将上述两种优化结合起来，加速比为1.64到2.65。0处理时间（s）0计算节点0cpu+图分区cpu+多粒度负载平衡0gpu+图分区gpu+多粒度负载平衡0图11. 整个堆芯装配模拟的处理时间改进。05.3 应用0物体的变形是由于温度和压力而最常见的机械现象。在反应堆压力容器中，加热整个堆芯装配将导致堆芯组件的变形。不恰当的加热环境设置将导致堆芯破裂，导致更换的成本。为了延长反应堆中堆芯的寿命，有必要明确实际环境中堆芯的结构力学变形过程。一种有效的方法是逐渐增加反应堆压力，并观察堆芯的变形。我们进行多步结构力学0在整个堆芯装配的大规模网格单元上进行仿真。这个仿真可能需要几天时间，传统的有限元方法很难对大规模情况进行并行计算。我们使用HTFETI来模拟和跟踪整个堆芯装配的变形演变。我们模拟了100亿个CEFR整个堆芯装配的网格单元，其中包含712个反应堆堆芯组件，如图1