arXiv:2611.14041v1 Announce Type: new Abstract: The time to solution of an iterative solver depends on the number of iterations and the cost per iteration. This paper develops algorithms that reduce both for elliptic problems discretized with high-order continuous finite elements on distributed GPUs, implemented in a unified matrix-free code supporting CUDA, HIP and SYCL. A distributed direct fast-diagonalization (FD) Poisson solver provides exact solutions on tensor-product meshes with direction-wise grading and coefficients varying along one coordinate direction. Its tensor data stay distributed, moved between slab layouts by collective transposes or kept in place by a Schur-complement approach. It is also an exact coarse solver. A two-level hybrid Schwarz preconditioner with FD cell solves (HS-FD) is extended to 2:1 hanging-node octrees with exact local cell inverses, while a \(p\)-multigrid preconditioner with Chebyshev--Jacobi smoothing (CJ-\(p\)MG) resolves hanging-node constraints at every level within fused kernels. Both approaches exploit FP32, where cell shapes and coefficients allow, within an FP64 conjugate-gradient iteration, adaptively select replicated or distributed coarse levels, support block solves for multiple right-hand sides and determine all parameters from the operator without case-specific tuning. Comparisons with Jacobi-preconditioned conjugate gradients (Jacobi PCG), MFEM, nekRS and deal.II span uniform, graded and octree meshes, layered coefficients, and four applications across NVIDIA, AMD and Intel GPUs. The direct solver is two to three orders of magnitude faster than Jacobi PCG on uniform meshes. HS-FD converges in 11 iterations for an octree generalized Poisson problem across 8 to 1024 MI250X graphics compute dies. The results guide solver selection by mesh structure and coefficients