A Two-Level Direct Solver for the Hierarchical Poincar\'e-Steklov Method

We introduce a two-level direct solver for the Hierarchical Poincar\'e-Steklov (HPS) method for solving linear elliptic PDEs. HPS combines multidomain spectral collocation with a direct solver, enabling high-order discretizations for highly oscillatory solutions while preserving computational efficiency. Our method employs batched linear algebra routines with GPU acceleration to reduce the problem to subdomain interfaces, yielding a block-sparse linear system. This system is then factorized using a sparse direct solver that employs pivoting to achieve better numerical stability than the original HPS scheme. For a discretization of local order $p$ involving a total of $N$ degrees of freedom, the initial reduction step has asymptotic complexity $O(N p^6)$ in three dimensions. Nevertheless, the high efficiency of batched GPU routines makes the overall cost for practical purposes independent of polynomial order (for order $p=20$ or even higher). Additionally, the cost of the sparse direct solver is independent of the polynomial order. We present a description and justification of our method, along with numerical experiments on three-dimensional problems to evaluate its accuracy and performance.