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

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.