A GEMM-based direct solver for finite-difference Poisson problems in non-uniform grids

TL;DR

This paper introduces a GPU-accelerated direct Poisson solver using GEMMs for non-uniform grids, offering improved performance and scalability for high-resolution simulations in computational physics.

Contribution

A novel tensor-based, GEMM-optimized direct solver for Poisson problems on non-uniform grids, integrating seamlessly with existing methods and hardware architectures.

Findings

Achieves superior time-to-solution compared to multigrid and FFT methods.

Demonstrates high parallel efficiency in strong scaling tests.

Enables efficient high-resolution simulations on modern heterogeneous systems.

Abstract

We present a direct Poisson solver for massively parallel simulations on three-dimensional Cartesian grids with non-uniform spacing. The method uses a tensor-based formulation in which the operator is diagonalized numerically along two directions through one-dimensional eigendecompositions, while the third direction is solved directly. The resulting dense transforms are evaluated efficiently as GEMMs (General Matrix--Matrix Multiplications), allowing many independent one-dimensional operations to be combined into matrix-matrix products that map well to modern CPU and GPU hardware. For uniform grids, the method reduces to the classical eigenfunction-expansion approach, and it naturally supports hybrid combinations of FFT-based and GEMM-based transforms depending on grid uniformity. After coupling the solver to an incompressible Navier-Stokes code, we assess its accuracy and performance against geometric multigrid and block cyclic reduction with FFT diagonalization. The results show that the proposed method is robust and consistently achieves the best time-to-solution. In strong scaling, the more compute-intensive GEMM-based variants attain higher parallel efficiency by better amortizing communication costs, while weak scaling highlights the expected trade-off between FFT-based and dense-transform formulations. Overall, the method enables efficient high-resolution stretched-mesh simulations on modern heterogeneous systems.