[论文解读] Matrix-free GPU implementation of a preconditioned conjugate gradient solver for anisotropic elliptic PDEs

Eike Mueller, Xu Guo+2 more|arXiv (Cornell University)|Feb 28, 2013

Advanced Numerical Methods in Computational Mathematics参考文献 9被引用 10

一句话总结

本论文提出了一种无矩阵、基于GPU加速的预条件共轭梯度（PCG）求解器，用于数值天气预报（NWP）中出现的各向异性椭圆型偏微分方程（PDE）。通过利用垂直方向上的各向异性特性，结合无矩阵模板计算、循环融合以及优化的内存访问方式，该方法实现了高吞吐量并降低了全局内存带宽，相较于大规模问题下的串行CPU实现和矩阵显式GPU实现，性能表现更优。

ABSTRACT

Many problems in geophysical and atmospheric modelling require the fast solution of elliptic partial differential equations (PDEs) in "flat" three dimensional geometries. In particular, an anisotropic elliptic PDE for the pressure correction has to be solved at every time step in the dynamical core of many numerical weather prediction models, and equations of a very similar structure arise in global ocean models, subsurface flow simulations and gas and oil reservoir modelling. The elliptic solve is often the bottleneck of the forecast, and an algorithmically optimal method has to be used and implemented efficiently. Graphics Processing Units have been shown to be highly efficient for a wide range of applications in scientific computing, and recently iterative solvers have been parallelised on these architectures. We describe the GPU implementation and optimisation of a Preconditioned Conjugate Gradient (PCG) algorithm for the solution of a three dimensional anisotropic elliptic PDE for the pressure correction in NWP. Our implementation exploits the strong vertical anisotropy of the elliptic operator in the construction of a suitable preconditioner. As the algorithm is memory bound, performance can be improved significantly by reducing the amount of global memory access. We achieve this by using a matrix-free implementation which does not require explicit storage of the matrix and instead recalculates the local stencil. Global memory access can also be reduced by rewriting the algorithm using loop fusion and we show that this further reduces the runtime on the GPU. We demonstrate the performance of our matrix-free GPU code by comparing it to a sequential CPU implementation and to a matrix-explicit GPU code which uses existing libraries. The absolute performance of the algorithm for different problem sizes is quantified in terms of floating point throughput and global memory bandwidth.

研究动机与目标

• 为解决数值天气预报（NWP）模型中椭圆型求解器的性能瓶颈问题，该类问题因具有各向异性、扁平的三维几何结构而计算成本高昂。
• 为GPU上的各向异性椭圆型PDE开发一种高性能、内存高效的PCG求解器，目标是面向实际NWP与地球物理建模工作负载。
• 通过无矩阵实现与PCG算法中的循环融合，降低全局内存带宽与访问开销。
• 通过标准CUDA库对比，展示该实现相较于串行CPU与矩阵显式GPU实现的优越性能。

提出的方法

• 采用无矩阵方法避免显式存储矩阵，在矩阵-向量乘积过程中动态计算局部模板。
• 基于椭圆算子的强垂直各向异性特性，构建物理启发式预条件器，以提升收敛速度。
• 应用循环融合技术，将多个PCG操作（如向量更新、点积、矩阵-向量乘积）合并为单个核函数，以减少内存流量。
• 实现采用CUDA核函数，线程块映射到水平网格平面，实现内存访问的合并与高效并行化。
• 通过BLAS调用对二维水平向量执行全局约简操作（如点积与范数计算），最大限度降低通信开销。
• 求解器设计具备可扩展性，可通过交换边界层（halo exchange）扩展至多GPU系统，最小化设备间数据传输开销。

实验结果

研究问题

• RQ1 无矩阵与循环融合的GPU内核在求解各向异性椭圆型PDE时，如何减少全局内存访问并提升性能？
• RQ2 在大气与海洋PDE中，通过预条件器利用垂直各向异性可获得多大的性能提升？
• RQ3 与使用CUDA库的矩阵显式GPU实现相比，无矩阵GPU PCG求解器在吞吐量与内存带宽方面表现如何？
• RQ4 在GPU架构上，循环融合在多大程度上减少了PCG算法的运行时间？
• RQ5 所提出的实现能否在多GPU集群上实现高效扩展，同时保持极低的设备间通信？

主要发现

• 无矩阵GPU PCG求解器在浮点运算吞吐量方面显著优于串行CPU与矩阵显式GPU实现，且全局内存带宽使用更低。
• 循环融合通过减少冗余内存访问并提升PCG内核中的内存合并程度，有效降低了运行时间。
• 无矩阵方法消除了存储大型稀疏矩阵的需求，降低了全局内存占用，提升了内存访问效率。
• 该求解器在典型NWP与地下水流模拟中的大规模三维各向异性问题上表现出高性能，且在aquila超级计算机上实现了良好的强可扩展性。
• 预条件器有效利用了垂直各向异性，使收敛迭代次数稳定在约100次，性能与多重网格方法相当。
• 该实现展现出极高的能效与绝对性能，适用于实际天气与气候模型。