基于CUDA并行技术加速2维矩阵MOC方法的研究

郑 勇 彭敏俊 安 萍 强胜龙 芦 韡

（1.中国核动力研究设计院 核反应堆系统设计技术重点实验室，四川 成都610213；2.哈尔滨工程大学，黑龙江 哈尔滨150001）

0 引言

随着反应堆物理计算对几何适应性和计算性能要求的提高，特征线方法在中子输运计算领域应用越来越广泛，而且由传统特征线方法出发，可以推得矩阵形式的特征线方程（MMOC）[1]，该方法仅需在构造系数矩阵阶段扫描特征线信息，而在随后的内迭代和源迭代过程中，均是对矩阵方程组的数值求解。由于减少了繁复的特征线扫描次数，因而，MMOC具有减少输运计算时间的潜力，但是否具备效率优势还依赖于矩阵方程的求解过程。以往的研究结果表明[2]，该矩阵方程具有较强的稀疏性和不规则的非零元分布，这些特点使得直接求解和常规迭代算法束手无策，而隶属于Krylov子空间方法的极小残余算法（GMRES）具有较好的收敛性能和计算效率[3]。

重启型的GMRES算法具有良好的收敛特性和数值稳定性，目前国际上已有大量的研究以期在算法层次上进一步减少GMRES算法的迭代次数。针对矩阵特征线方法中系数矩阵的特点，本文基于GPU并行技术和CUDA程序开发规范，从代码层面上对重启型GMRES算法进行并行化，以提升该算法的计算效率，从而加速矩阵MOC程序的求解。

1 矩阵特征线方法概述

特征线方法是将偏微分方程简化为常微分方程，然后使用成熟的数学方法求解常微分方程的一种方法。在反应堆中子输运计算领域，特征线方法在每个离散方向上生成一族射线覆盖待求解区域，并跟踪每条射线来研究中子通量的分布规律。沿着每条射线，中子输运方程可以写成如下的特征线形式：

式中，ψ（s）为中子角通量；s为沿射线的当地坐标；∑t和q（s）分别为宏观总截面和中子源项。假设宏观截面和源项在计算网格内呈均匀分布时，可以解析求解上式，得到中子出射角通量和入射角通量之间的关系：

由方程（2）不难看出，根据特征线在相邻网格之间的角通量连续条件，可以从外边界入射角通量依次求得任意网格的入射和出射角通量，从而获得中子角通量的传播方程：

结合中子角通量和通量密度之间的加权求和关系，可以建立起任意计算网格上中子通量密度同外边界入射角通量和网格源项的方程；同时利用外边界条件，还可以建立各条特征线外边界处入射角通量之间的关系。从而，构成了待求解的矩阵方程：

式中，S1、S2、S′1、S′2、S″1、S″2均为扫描特征线得到的系数矩阵，需要在迭代计算开始前计算完毕；Φ和Ψ分别代表由中子通量密度和外边界入射角通量组成的向量；E为单位矩阵；Q表示无自散射源项的其他源项之和组成的向量。

2 极小残余算法的分析

由特征线方法得到的矩阵方程通常规模较大，且系数矩阵具有较强的稀疏性，难以应用直接求解算法，同时非零元的非对称分布也限制了迭代算法的可选择范围。以往的研究表明，重启型极小残余算法（GMRES）具有比较稳定的收敛特性，若能选用恰当的预条件技术，还能够进一步加速其收敛过程。下面给出带左预条件的重启型GMRES算法的标准过程（其中，M为预条件矩阵）：

在矩阵特征线方法中，当子空间维度m不大于50时，算法即可收敛。因此，该算法中最耗时的部分是用于构造正交基的Arnold循环，而最小二乘问题的求解时间几乎可以忽略不计，表1分析了Arnold过程中涉及向量、矩阵的各种基本操作的计算量。

表1 Arnold过程基本操作计算量

假设矩阵方程的阶数为n，而矩阵A及M-1的非零元个数分别为nza、nzm，那么向量-向量与向量-标量运算的总浮点计算量为（2m2+7m+4）n，矩阵-向量乘操作需要的总计算量为4（m+1）（nza+nzm）。不妨定义稀疏度sp=（nza+nzm）n，则可以得到两种运算量的比值为：

如前所述，当子空间维度m∈（0，50）时，f∈（1，26），因此当稀疏度sp小于26时，向量运算的浮点计算量可能会更大，反之则必然是SpMV的浮点计算量更大。通常情况下，特征线所穿过的计算网格数量以及覆盖网格的特征线条数都大于26，因此矩阵的稀疏度sp通常会远大于26。

3 基于CUDA技术的GMRES算法并行优化

SpMV运算可以分成如算法1所示的两步，首先将矩阵非零元素与向量对应元素相乘，并把乘积保存在临时矩阵temp中；然后对每个行向量进行累加求和，得到目标向量对应的元素。

表2 CSR格式的稀疏矩阵-向量乘算法

在第1步中，由于对向量x的访问模式是随机的，因此，可以将它绑定到访存效率较高且适合随机访问的纹理内存中，其访存延迟可以通过激活海量线程而有效隐藏。该算法的并行优化关键在于第2步，由于稀疏矩阵每一行的非零元个数都是不相等的（如图1所示），如何同时高效地求得临时矩阵temp每个行向量的累加和是优化的难点。

图1 系数矩阵中各行非零元个数

针对稀疏矩阵不同部分稀疏程度的不同，本文结合文献[4]提出合并访存的方法，采取分部并行优化策略：第2部分由于每行平均非零元个数较少，因此，可以直接采用合并访存和共享内存技术实现CUDA并行；对第1部分提出了两种CUDA并行计算优化方案，在合并访存的基础上将每行的求和操作分别布局到线程和线程块上。以上三种并行方案，在后文中分别简写为：gpu_cv/sm、gpu_thread及gpu_block。

4 数值验证与讨论

为了评估并行程序的加速性能，本节对2维C5G7问题进行了计算。计算环境参数配置如表3所示，虽然两块GPU均支持CUDA编程，但本文仅采用Tesla专业计算显卡作分析比较。

表3 计算环境参数配置表

国际上，C5G7基准题被广泛用于输运程序计算精度的验证，本文针对C5G7基准题的1/4堆芯进行计算，以比较上述几种方法在处理较大规模问题时的加速效果。计算区域共包含4个17×17的燃料组件和外围的5个反射层组件，采用文献[5]给出的7群截面参数。计算条件为16个方位角，2个优化的Leonard极角[6]，特征线间距为0.05 cm，采用粗网有限差分求解器加速；求解器参数设置为子空间维度为10，允许的最大迭代次数为10。

计算结果表明，无论采取串行还是并行计算，2维C5G7基准题的计算结果都是完全相同的，且并行方案对计算精度几乎没有影响，部分结果如表3所示，可以发现与国际上其他输运计算结果相比[5]，本文给出的计算精度完全在可以接受的范围内。

表5给出了采取不同并行方案的耗时，同时CMFD加速方法和Krylov子空间维度对计算效率的影响也在表中给出。从表中可以看出，无论是cpu串行计算还是gpu并行计算，CMFD都能带来超过6倍的加速比。当子空间维度越小，GMRES（m）求解器需要的重启动次数越多，自然也就需要更多的SpMV操作；而子空间维度增大后，对向量的操作也会呈平方关系增加，这种平方增加关系不允许选择太大的子空间维度。与gpu并行程序比较，在cpu串行计算中，子空间维度对计算效率的影响更为显著，这是由于gpu并行程序对向量的操作耗时几乎可以忽略不计。

表5 2维C5G7问题运行时间比较/单位：秒

图2给出了以cpu串行计算为参考，不同计算条件下获得的加速比，图中的Case编号对应于表4中的第2列。由图中可以看到，不同的并行策略对加速比有显著影响，优化后的并行方案获得的加速比成倍增加，且gpu_block具有最高的计算效率，能够达到6倍以上的加速比。在每种并行策略中，CMFD加速方案对加速比均无明显影响。随子空间维度增加，加速比有轻微的变化，这种变化主要是因为cpu串行计算效率对子空间维度更为敏感。

图2 2维C5G7问题不同计算条件下的加速比

表4 C5G7问题计算精度比较

5 结论

本文首先简要介绍了矩阵特征线方法的基本模型，引出了求解矩阵方程的极小残余算法，继而对该算法中涉及矩阵和向量的基本运算量进行了分析。理论分析结果表明，稀疏矩阵-向量乘操作是算法中的热点，开发并行求解器时时需要重点优化。基于CUDA并行技术，采用访存合并和共享内存并行方案，编制了GMRES算法的GPU并行版本，并在此基础上提出了两种并行优化策略。为了验证并行策略的有效性，评估并行求解器的加速效果，对2维C5G7基准题进行了对比计算。计算结果表明，并行方案的不同对计算精度没有显著影响，而优化后的并行程序较原始并行程序有更高的加速比。