DEPARTMENT OF COMPUTATIONAL AND DATA SCIENCES
Ph.D. Thesis Colloquium
Speaker : Mr. Gourab Panigrahi
S.R. Number : 06-18-01-10-12-21-1-19282
Title : “Hardware-Aware Matrix-Free Algorithms for Accelerating Large-Scale
Finite-Element Simulations on Distributed CPU and GPU Architectures”
Research Supervisor : Dr. Phani Motamarri
Date & Time : July 30, 2026 (Thursday), 15:30 PM
Venue : #102, CDS Seminar Hall
ABSTRACT
Most large-scale simulations governed by partial differential equations (PDEs) reduce, after numerical discretization, to two recurring computational problems: eigenproblems solved by subspace iteration approaches, and linear systems by preconditioned Krylov methods. On present-day supercomputers, the performance of the underlying solvers is limited by data movement through deep memory hierarchies and across interconnects rather than by floating-point throughput. This thesis addresses these issues within the framework of higher-order finite-element (FE) discretizations, where their dominant kernels are the action of the FE-discretized operator on large blocks of vectors within eigensolvers and the solution of FE-discretized generalized Poisson or Helmholtz systems within linear solvers. These kernels have traditionally relied on assembled global sparse-matrix or cell-matrix strategies, both of which scale poorly at high FE polynomial orders: the former limited by data movement, the latter by an operation count that grows as the sixth power of the polynomial order. Addressing these limitations, we develop hardware-aware, matrix-free algorithms that evaluate operator actions on the fly through sum factorization, exploiting the tensor-product structure of the FE basis to reduce both operation count and memory footprint. Combined with communication-efficient approximation-tolerant eigensolvers, tensor-structured direct linear solvers, and device-native multigrid preconditioners, the performance of these methods has been demonstrated on the world’s foremost supercomputing platforms. As the driving application, we adopt the FE discretization of Kohn–Sham density functional theory (DFT), one of the most demanding HPC-intensive workloads in computational materials science, on which these algorithms are shown to substantially accelerate the time-to-solution.
We first develop matrix-free algorithms for FE-discretized Helmholtz matrix–multivector products on multi-node CPU and GPU architectures, where existing implementations are not designed for operator actions on large blocks of vectors. To this end, we design batched evaluations with hardware-tailored batch sizes to improve data locality, together with even–odd decomposition and SIMD vectorization on CPUs and kernel fusion, constant-memory strategies and register blocking on GPUs, overlapping compute with data movement. For the operator action on 1024 vectors, the implementation is up to 2.8× faster than the cell-matrix baseline on a single V100 node of OLCF Summit and 2.6× on a single Intel Xeon node of PARAM Pravega, and up to 17× and 6.5× faster than the single-vector deal.II matrix-free implementation and libCEED–MFEM, respectively. Embedded in a Chebyshev-filtered subspace iteration (ChFSI) eigensolver for the Helmholtz operator, the speedups over the cell-matrix approach reach 4.0× and 2.2× on a single CPU and GPU node, and 3.7× on ~3000 Xeon cores of Pravega and 1.4× on ~25 V100 GPUs of Summit. This framework also underpinned the DFT electrostatics Poisson solves in the ACM Gordon Bell Prize-winning SC’23 study of quasicrystals and interacting extended defects, which scaled to 64,000 AMD MI250X GCDs of the OLCF Frontier exascale system.
Building on these ideas, we develop a matrix-free algorithm for the FE-discretized Kohn–Sham DFT eigenproblem on distributed CPU architectures. The underlying local and nonlocal operator actions on multivectors are recast as structured tensor contractions over one-dimensional FE basis functions and quadrature data. A unified batched layout is proposed that handles the action of both real- and complex-valued operators on multivectors while maximizing cache reuse and SIMD utilization across AVX2, AVX-512, and SVE, complemented by even–odd decomposition and mixed-precision intrinsics. Matrix–multivector product benchmarks on OLCF Frontier (AMD EPYC), PARAM Pravega (Intel Xeon), and RIKEN Fugaku (Fujitsu A64FX) show speedups of up to 4–6× on representative material systems over the cell-matrix baseline of the DFT-FE library, with further end-to-end gains from a residual-based ChFSI (R-ChFSI) eigensolver that tolerates inexact matrix–vector products without sacrificing robustness.
Taking this further to GPU architectures, we fuse the entire discretized Kohn–Sham Hamiltonian operator action into a single kernel that combines a register–shared-memory hybrid sum factorization, constant-memory shape-function matrices, even–odd decomposition, and a precision-aware batched layout, distributed across ranks by the Morton (Z-order) numbering for locality. FE constraints arising from boundary conditions or from hanging nodes on non-conforming meshes are handled inline to avoid atomic operations, and a reduced-precision (BF16) halo exchange is overlapped with computation without affecting the self-consistent field (SCF) convergence of the underlying DFT nonlinear eigenproblem. TF32 tensor cores and A100 asynchronous copies (cp.async) further improve compute throughput and memory bandwidth. On the NVIDIA A100 GPUs of ALCF Polaris, the filtering step in R-ChFSI eigensolver achieves speedups up to 3.9× on representative benchmarks, sustained in strong scaling to 128 A100 GPUs.
We now turn to the linear systems arising in FE discretizations. When the mesh retains a global tensor-product structure, iterative solution can be avoided altogether. We develop a distributed-memory GPU fast-diagonalization (FD) solver that inverts the assembled FE Laplacian exactly in a single, non-iterative solve at O(N^(4/3)) FLOPs, where N is the number of degrees of freedom, supporting Dirichlet, periodic, and semi-periodic boundary conditions on rectilinear meshes. For the singular fully periodic case, we deflate the constant mode rather than pinning a node, keeping the tensor structure intact and avoiding the ill-conditioning that pinning introduces. With O(N/P) compute and memory on each of P ranks, the solver scales efficiently on distributed GPUs. Against a conjugate-gradient (CG) solver with Jacobi preconditioner, the proposed FD solver is 125–140× faster on uniform meshes. On an atomistic-scale battery electrode–electrolyte interface (LLZO) discretized on an FE mesh (12.3 million DoF, periodic in two directions, Dirichlet in the third) and solved to machine precision on 64 GPUs of Polaris, the speedup is about 500×, rising to 4600× on strongly graded meshes. This advantage widens as the mesh becomes more graded, because the cost and accuracy of the solver are independent of the mesh-induced condition number.
On general meshes, where the tensor-product structure is lost and iterative solution is unavoidable, we develop a GPU-accelerated matrix-free preconditioned CG solver. It fuses the Laplace kernel with the CG iteration, uses tree-reduction inner products and a Chebyshev–Jacobi preconditioner free of global reductions. On 32 nodes (384 Intel Data Center Max GPUs) of the ALCF Aurora exascale system, it attains 5.86× higher throughput than the deal.II matrix-free implementation. Since the iteration count of this single-level preconditioner grows with FE polynomial order, we extend it to adaptively refined meshes by developing a fully device-native p-multigrid preconditioner whose iteration count is nearly independent of the polynomial order. We further test this preconditioner on the variable-coefficient generalized Poisson operators arising in the continuum solvation models used for electrochemical systems in DFT, for which no fast direct solver exists. On a mesh that models the Trp-cage protein (304 atoms, 21.4 million DoF) in implicit solvent, it reduces CG iterations from 2583 to 147, a 17.6× reduction.
Overall, this thesis develops a set of hardware-aware algorithms for accelerating the two classes of computational problems that arise in FE discretizations of PDEs, namely eigenproblems and linear systems. Implemented within the DFT-FE library and validated against production baselines at desired accuracy, these methods reduce the time-to-solution of large-scale simulations from single-node CPUs to exascale GPU machines, in several cases by orders of magnitude. The kernels, data layouts, and solver designs make no assumptions specific to DFT, and therefore carry over directly to higher-order FE discretizations in other scientific areas on modern HPC architectures.
ALL ARE WELCOME
