技术分享︱大型稀疏线性方程组求解技术——工业仿真的底层核心
一、背景
在工业仿真领域，对各种现实世界的问题进行数值模拟时，如流体动力学分析、电磁场仿真、结构力学应力应变分析等，其控制方程通常是偏微分方程组，在经过不同方法的隐式离散之后最终都可转化为大型稀疏线性方程组。随着人们对计算精度要求的不断提高，方程组的阶数也从上千阶、几十万阶提高到百万、千万阶甚至更高，所需的计算量以及存储需求也随之迅速膨胀。根据一般经验，方程组求解时间会占总计算时间的70%以上，往往是整个计算过程中的性能瓶颈。如果说求解器是工业CAE软件的核心模块，那么大型稀疏线性方程组的求解技术将毫无疑问是底层求解器的核心。
NASA翼型网格经过离散得到的稀疏矩阵（素材来源于网络）
二、方法
众所周知，稀疏线性方程组的求解方法可以分为直接法和迭代法 ，两类方法各有优劣，特点比较如下：
迭代法[1]：
对于不同类型稀疏矩阵表现差异较大，存在收敛性与收敛速度问题，催生了许多预处理技术（Preconditioners）；
对原矩阵的编辑很少，SpMV(Sparse matrix-vector multiplication)是其核心运算；
内存需求小，求解速度较快，算法复杂度低；
较易实现并行化。
直接法[2]：
通用、稳定；通过前后处理，能够保证计算的收敛性与精度；
对原矩阵的编辑多（分解、排序、缩放等）；
内存需求大，求解速度慢，算法复杂度更高；
并行度有限。
其中迭代法的种类很多，可以分为定常（Stationary）迭代法与非定常迭代法[3]。经典的定常迭代法有Jacobi迭代法、Gauss-Seidel迭代法、SOR迭代法等，这些方法均可基于矩阵分裂推导得到；而在数值模拟中，非定常迭代法则显得更加重要，常见的有共轭梯度法（Conjugate Gradient CG）、广义最小残量方法（Generalized Minimal Residual GMRES）、稳定双共轭梯度法（Biconjugate Gradient Stabilized Bi-CGSTAB）等，这些算法都属于Krylov子空间方法。 如今高性能共轭梯度（HPCG）测试包已成为国际上评测超级计算机性能的主要工具[4]，Krylov子空间方法也被评为20世纪最伟大的十大算法之一[5]。 还有一些迭代算法则比较特殊，它们基于特殊结构和性质，例如我们耳熟能详的多重网格（Multigrid MG）方法，它更多地被称为数值计算领域中一种加速迭代收敛的技术，而不仅仅是一种单纯的算法。
2014-2021年HPCG性能评测结果对比（素材来源于网络）
直接法的基础是矩阵的分解，常见的分解形式有LU分解、Cholesky分解、QR分解等。稀疏线性方程组的两类常见直接求解算法分别为超节点（Supernodal）方法和多波前（Multifrontal）法，其主要思想是将完整的稀疏矩阵的分解任务转化成许多个相对稠密的子矩阵的分解任务，任务间的依赖关系由消去树（Elimination tree）或其他类似的数据结构来确定。直接法的求解步骤通常分为矩阵重排、符号分解、数值分解与回代求解四个部分。
一个稀疏矩阵与其对应的消去树（来自文献6）
三、挑战
当前，国产超级计算机的峰值性能已达每秒十亿亿次量级，不久便将进入百亿亿次（E级）计算时代，我国的神威E级计算机和天河E级计算机已经蓄势待发。这些国际领先的超级计算机为我国科学与工程计算应用迈进超大规模计算时代、实现更高精细度的数值模拟提供了强力支撑。然而，超大规模计算也给高实用性与高性能的大型稀疏线性方程组求解的算法设计与优化带来了巨大挑战。
1.高实用
当前，国产超级计算机的峰值性能已达每秒十亿亿次量级，不久便将进入百亿亿次（E级）计算时代，我国的神威E级计算机和天河E级计算机已经蓄势待发。这些国际领先的超级计算机为我国科学与工程计算应用迈进超大规模计算时代、实现更高精细度的数值模拟提供了强力支撑。然而，超大规模计算也给高实用性与高性能的大型稀疏线性方程组求解的算法设计与优化带来了巨大挑战。
以SiP封装芯片的电磁-热-力耦合数值模拟为例，其稀疏矩阵具有明显的病态特征（来自文献7）
2.高性能
随着计算机硬件性能的提升，超级计算机呈现 “多级嵌套并行、异构众核加速” 的复杂体系结构特征，会导致大型稀疏线性方程组求解器的实现效率急剧下降。从下图也可以发现，随着高性能计算机系统变得更加复杂，特别是众核架构采用后，在每一代世界性能最为强大的超级计算机上，应用的实际求解能力变得更加低效，即解决问题时间（time-to-solution）变得越来越长。如何设计能 匹配机器体系结构特征 的算法与性能优化技术，是大型稀疏线性方程组求解技术以及其他科学计算核心算法中当前亟待解决的关键问题。
解决问题时间与超级计算机性能趋势对比
对于大规模稀疏线性方程组，原有串行和小规模并行模式下的数据结构和算法容易导致并行求解性能低下或失败。在分布式并行层面，需要解决以下几个问题：一是在数据和任务分解方面，如何设计良好的负载均衡策略、稀疏矩阵的高效存储格式以及计算通信重叠等优化策略；二是在负载均衡的前提下，如何设计以尽力避免节点间的通信；三是在内在串行特性导致并行化困难的算法方面，如何改进数据的分布方式以增加并行性。
在共享内存环境中，稀疏线性方程组求解算法的可扩展性问题也需要特别关注。因为现代多核/众核处理器上的核数在可预见的未来也将越多，在单个CPU上封装数十甚至上百个有较强处理能力的核心，或是在GPU上封装成千上万个轻量级处理单元将变得非常普遍。如何在这种共享内存节点上实现细粒度的并行仍然是很有挑战的研究内容。
左图：AMD霄龙CPU，64核128线程， 右图：英伟达Hopper架构GPU，1.8万核心
四、我们的探索——UNAP
面对来自“应用与机器”的双重挑战，神工坊团队与国产异构众核平台体系架构紧密耦合，发展了一套大型稀疏线性方程组求解库UNAP。该求解库是早期应基于非结构网格的仿真需求而开发的，其全称为“UNstructured Algebra Package”。
为了解决高实用性的挑战，UNAP结合异构众核处理器多级并行的特点和稀疏矩阵迭代解法的需求，初步探索了各种预处理方法在众核异构平台上的并行实现技术。UNAP已实现PCG、PBiCGStab、GMRES等Krylov子空间方法以及AMG代数多重网格求解器，预处理器包含Jacobi、DIC、DILU等，未来还将开发直接求解模块，以满足来自各领域的复杂应用需求。
为了达到高性能的目标，UNAP根据国产超算的异构特点，结合其处理器的多级缓冲区，实现了计算/通信混合的并行迭代算法；由于迭代算法的并发度天然较高，UNAP主要通过算法调整、浮点运算替换内存访问、通信同步等手段[9]，充分利用了申威主从众核浮点计算性能高而带宽受限的特点，且通过增加计算比例降低了全局集 合通信代价。在共享内存层面，UNAP还可以调用太湖之光超级计算机上的加速工具套件UNAT和向量计算加速库swArrays，充分发挥出从核阵列的计算能力，达到进一步的性能提升。
代数求解库UNAP的组织架构
UNAP的整体结构如上图所示， 底层并行环境主要包含对“神威·太湖之光”超级计算机申威众核异构芯片SW26010的从核并行计算支持和普通的进程间MPI层级并行支持。采用以容器模板为核心的架构体系， 基本容器层主要包含Vector分布式向量容器以及Matrix分布式矩阵容器。在容器层向下可调用向量计算加速库和非结构计算加速库，实现在国产申威众核芯片上的高效并行，同时，向上可扩展预条件子、迭代解法以及直接解法。这种架构可以将元素的内部处理通过模板泛化，从而降低求解算法实现的复杂度，也具有较好的 可移植性 ，易于兼容x86-CPU平台。核心的代数求解器模块包含了Krylov子空间迭代求解器、预条件子、代数多重网格求解器以及开发中的直接求解模块。
几个案例说明UNAP的应用情况：
发动机燃烧室大涡模拟： 在某航空发动机全环燃烧室的大涡模拟中，网格量达10亿，采用UNAP作为核心求解模块后，最终的并行规模达到1.6万进程，稀疏线性方程组求解部分加速达到20倍，收敛速度显著提升。
船舶水动力学应用： 在某水动力学软件中，通过UNAP的接入（主要使用了求解压力和压力修正方程的 代数多重网格算法 和求解速度等方程 的预条件稳定双共轭梯度法 ），获得了在代数求解过程中的自动多级并行能力，在神威·太湖之光超级计算机上完成了千万级网格对标算例计算， 并行规模达到万核级别，相对百进程 并行效率不低于50% ，经测试与商业软件FLUENT相当。
电机设备电磁场分析： 在某公司核心求解器向神威平台的移植部署中，采用UNAP代替原有直接求解库，进行了网格量为千万级的电机模型有限元算例并行计算测试，并对比了替换前后算例的节点磁通密度计算结果，UNAP表现良好。
一个典型算例展示UNAP的功能与性能：
算例 ： 方腔驱动流，不可压，上壁面滑移速度为1，其他为固壁边界，方腔的边长为1，雷诺数100。
并行计算结果： 以4进程为例，下图给出了方腔在各个进程的分割情况，流线图中可以清楚看到主涡和二次涡。
方腔驱动流的流线图
简单的并行计算效率测试：
强扩展性测试：
固定算例网格整体求解规模，通过比较不同并行核数下的计算时间可获得 强扩展性 并行计算效率。本测试采用的网格规模数量为2千万，网格类型为六面体，进程数量分别为64、256、1024，并行核数规模分别为4160、16640和66560，统计时间为前20步计算时间。
弱扩展性测试：
固定每个进程上的网格规模，通过比较不同并行核数下的计算时间来获得 弱扩展性 并行计算效率。本测试采用的网格规模数量分别为2千万、4千万、8千万，网格类型为六面体，进程数量分别为64进程、128进程和256进程，并行核数规模分别为4160、16640和66560，单个进程网格数量约为31.25万，统计时间为前20步计算时间。（本文作者：赵程鹏）
参考文献：
[1] Saad Y. Iterative methods for sparse linear systems[M]. Society for Industrial and Applied Mathematics, 2003.
[2] Davis T A, Rajamanickam S, Sid-Lakhdar W M. A survey of direct methods for sparse linear systems[J]. Acta Numerica, 2016, 25: 383-566.
[3] Barrett R, Berry M, Chan T F, et al. Templates for the solution of linear systems: building blocks for iterative methods[M]. Society for Industrial and Applied Mathematics, 1994.
[4] Marjanović V, Gracia J, Glass C W. Performance modeling of the HPCG benchmark[C]//International Workshop on Performance Modeling, Benchmarking and Simulation of High Performance Computer Systems. Springer, Cham, 2014: 172-192.
[5] Cipra B A. The best of the 20th century: Editors name top 10 algorithms[J]. SIAM news, 2000, 33(4): 1-2.
[6] Gupta A, Karypis G, Kumar V. Highly scalable parallel algorithms for sparse matrix factorization[J]. IEEE Transactions on Parallel and Distributed systems, 1997, 8(5): 502-520.
[7] Wang W, Liu Y, Zhao Z, et al. Parallel Multiphysics Simulation of Package Systems Using an Efficient Domain Decomposition Method[J]. Electronics, 2021, 10(2): 158.
[8] 刘伟峰. 高可扩展, 高性能和高实用的稀疏矩阵计算研究进展与挑战[J]. 数值计算与计算机应用, 2020, 41(4): 259.
[9] GU H, REN H U, LIU C, et al. An optimized Chebyshev smoother In GAMG solver of openfoam on sunway Taihulight supercomputer[C]//The 13th OpenFOAM Workshop. 2018.
