high performance computing engineers

Abstract: Historical experience can guide scientific coding at multiple levels, but task-level retrieval may miss implementation lessons that become relevant during execution. We introduce MUSE, a memory management framework combining multi-level memory with selective runtime injection for scientific coding agents. MUSE builds global, approach-level, and implementation-level memory from historical trajectories. At initialization, the agent receives global memory and problem-matched approach-level memory. During execution, current interactions are used to retrieve implementation-level memory. An LLM gate selects relevant experience and generates specific review advice, or skips injection. Across four benchmarks spanning numerical problem solving, kernel optimization, and scientific software repair, MUSE outperforms retrieval and memory baselines on all evaluation metrics. Ablations support the multi-level design and gated injection, while transfer experiments show that memory remains useful across kernel benchmarks. These results show that multi-level memory and selective runtime injection enable effective reuse of experience in scientific coding. Code is available at https://github.com/hr-jin/MUSE.

Abstract: We give an exposition of OpenAI's recent proof that the exponent $ω_\mathbb{C}$ of matrix multiplication over the complex numbers satisfies $ω_\mathbb{C}\leq9/4$. Our version of the proof gives $ω_\mathbb{F}\leq9/4$ over every field $\mathbb{F}$. We also treat the rectangular matrix multiplication exponent $ω_\mathbb{F}(a,b,c)$. Here we have $ω_\mathbb{F}(1,k,1)=2$ for $0\leq k\leq1/2$ and $ω_\mathbb{F}(1,k,1)\leq1+k+1/(4k)$ for $k\geq1/2$, and, for real numbers $a\geq b\geq c\geq0$, $ω_\mathbb{F}(a,b,c)=a+b$ if $c(a+b)\leq ab$ and $ω_\mathbb{F}(a,b,c)\leq(a+b+c)/2+(ab/c+ac/b+bc/a)/4$ otherwise. Our exposition is in the language of degeneration inequalities among polynomial multiplication tensors and Strassen's asymptotic spectrum duality. We moreover prove an operational barrier (which can be thought of as coming from an operational asymptotic spectrum) showing that the current analysis of all these upper bounds on $ω_\mathbb{F}$ and $ω_\mathbb{F}(a,b,c)$ is optimal and that in fact the current degeneration ingredients cannot improve the bounds.

Abstract: We design and optimize a quantum error correction (QEC) decoder for NVIDIA GPUs, combining Relay belief propagation, a state-of-the-art message-passing algorithm, with order-zero ordered statistics post-processing. Our implementation uses CUDA Graphs to control GPU execution, stops updates for converged samples, and collects unresolved samples in queues for post-processing. We optimize sparse traversal through coalesced message accesses, prefetching, and register reuse, and accelerate binary post-processing with warp-level collectives. In addition, mixed FP16-FP32 arithmetic and INT8 quantization reduce message storage and traffic while retaining accuracy using accumulators at higher precision. We evaluate surface, hypergraph product, and bivariate bicycle code workloads on NVIDIA A100 and GH200 data-center GPUs and a commodity RTX 5080. On the bivariate bicycle code at error probability 0.003, INT8 halves amortized decoding time relative to independently tuned FP32 with a negligible logical error rate increase on RTX 5080. With independent syndrome batches submitted back-to-back to concurrent decoder instances, RTX 5080 delivers 91% of GH200's highest measured throughput with an 8.9x shorter median submission-to-correction latency. The combination of performance optimization, mixed precision, INT8 quantization, and concurrent decoder instances opens a path toward QEC decoder clusters built from commodity GPUs.

Abstract: Stencil computations are a fundamental class of workloads in scientific and engineering applications. Although Tensor Cores offer new opportunities for stencil acceleration, existing tensorization methods suffer from substantial redundant computation and insufficient data reuse. We propose PrismStencil, a stencil tensorization algorithm that introduces a novel mapping scheme between matrices and stencil grid points. PrismStencil employs multidimensional vector mapping to accommodate the multidimensional data dependencies of stencil computations and a unified gather--scatter tensorization scheme to exploit the relative sizes of the matrix dimensions in Tensor Core operations, thereby reducing the redundant computation introduced by tensorization. Moreover, it introduces a cross-product data-reuse mechanism that exploits the overlap among neighboring stencil updates to further enhance data reuse. Across representative stencil kernels on an NVIDIA A100 GPU, PrismStencil achieves average speedups of 2.13x and 2.17x over state-of-the-art Tensor Core stencil methods with FP64 and TF32 precision, respectively.

Abstract: The Polynomial Modular Number Systems (PMNS) aim to represent elements of fields or rings of large characteristics using polynomials satisfying bounds on some parameters (degree, absolute values of the coefficients). Those PMNS, while some conditions are fulfilled, using convenient parameters and implementation features, allow some speed-ups in cryptographic computations. Recent works (Meloni \emph{et al.} \cite{MeloniPV25}) propose better speed-ups by improving the parameter generation of the system, and by using multi-precision polynomial coefficients, i.e. coefficients stored using several machine words. In this work, we first present PMNS schoolbook software implementation better in some context than the Toeplitz state-of-the-art counterpart of \cite{MeloniPV25}, taking advantage of smaller memory cost by being free to choose a smaller $n$ parameter, minimizing at the same time the complexity of the multiplication. This implementation for 4096 bit modulo size shows element size smaller by 11\% and 20 \% speed-up. We thus present a new improvement on the PMNS, applying to the \textsl{internal reduction} an approach similar to the truncated Montgomery reduction technique presented by Didier \emph{et al.} in \cite{DidierEGR24}. In the context of software implementations using \texttt{AVX512} instruction set extension, and modulo size up to 8192 bits, this new approach allows speed-ups up to 15\% in modular multiplication computation, in comparison with conventional approaches.

Abstract: The performance of stationary iterative methods is memory-bound, due to the low arithmetic intensity. So, we propose a mixed-precision technique to reduce the amount of data movement in matrix-splitting iterative methods. We demonstrate both theoretically and experimentally that this technique will generally converge at a similar rate as an exact-precision iteration, depending on the precision used and the conditioning of the $M$ matrix. Furthermore, we demonstrate that iterative refinement can be used to refine the solution to full double-precision accuracy. Using this approach, we mixed half and double precision to achieve an average speedup over uniform double precision of $1.60\times$ for scalar and block Jacobi and $1.25\times$ for scalar and block Gauss-Seidel. Furthermore, we reduced the precision in the SSOR iteration used by OVERFLOW, a NASA code for compressible fluid dynamics simulations, and achieved a $1.14\times$ speedup to the overall application.

Abstract: Strassen showed that two 2x2 matrices can be multiplied with 7 multiplications instead of 8. Applied recursively, his algorithm multiplies two nxn matrices with O(n^2.807) multiplications, beating the naive O(n^3). The best known 3x3 recursive matrix multiplication algorithm uses 23 multiplications O(n^2.854). The best published lower bound of 21 (on algorithms with integer constants) leaves room for an algorithm with O(n^2.771) multiplications, and thus does not rule out the possibility of an algorithm that would beat Strassen's. We prove a lower bound of 22 multiplications for any 3x3 recursive algorithm with integer constants, proving that no such algorithm can do better than O(n^2.814) multiplications, and eliminating the possibility of a 3x3 algorithm that beats Strassen's 2x2 method. The proof builds on a recent decomposition method from Wang, who approached the problem by turning it into 496 subproblems. We provide exact solutions for 359 of them. The proof is in Lean; verification requires auditing only a few short files. The Lean formalization directly encodes statements about the limitations of recursive algorithms for matrix multiplication, as opposed to just a statement about the rank of the problem.

Abstract: We validate memory-optimal cost functions for transformer kernels derived via the Mathematics of Arrays (MoA). Companion Papers I-IV formally derive kernels for attention forward, backward, fused forward+backward, decode, and the complete block (RMSNorm, gated MLP) as a hardware-independent specification (DNF) transformed to a machine-specific realization (ONF) via gamma, with verification to machine precision against PyTorch. This paper checks those predictions against measured performance on two HPC clusters (Purdue Anvil, NCSA Delta) across CPU and GPU. Three results stand out. (1) We identify and fix a GPU regression: fusing forward+backward, proven to avoid materializing an O(n^2) intermediate, initially ran slower than naive on GPU due to atomic contention. Profiling confirmed 2.00x more atomic instructions; a targeted ONF rewrite reversed it, yielding up to 2.5x speedup. (2) Identical derivations produce markedly different real costs by topology: 535x NUMA-locality penalty on one cluster vs <3x oversubscription on another, showing optimal deployment is a function of the machine's array structure. (3) We report a partially resolved anomaly: identical denotational computations run faster in C than Fortran on CPU but faster in Fortran than C on GPU, narrowed to one dominant kernel and one memory-latency stall mechanism (3.17x time gap matches 3.35x stall gap). We treat hardware-specific optimization as a routine ONF rewrite with fixed, verified DNF, a candidate methodology for scaling AI onto evolving hardware without re-deriving correctness.

Abstract: Verifying a numerical algorithm in a large scientific simulation framework is challenging: the framework is too big to model-check, and the unit tests exercise only sampled inputs. We report a case study in which we developed a new cloud-in-cell (CIC) deposition algorithm for Flash-X, a large-scale multiphysics simulation framework, keeping the CIVL model checker in the development loop. Rather than verify the algorithm within Flash-X's hefty infrastructure, we extract only the interfaces that the algorithm needs into a small, self-contained C model, which abstracts away implementation details of the Flash-X infrastructure that are unrelated to the new algorithm. The new algorithm is then built and checked within this C model. The CIVL model checker enables verification of the required physical properties of the CIC deposition algorithm using symbolic values for the particle positions. It proves two physical properties---mass conservation and the deposition location---for the continuum of admissible positions in the simulation domain. CIVL also verifies the algorithm's memory safety and freedom from MPI deadlocks and data race conditions over all rank distributions within specified bounds. Writing the verifying properties first and continuously checking them at each stage of the bottom-up prototyping workflow turned CIVL into a design guardrail that greatly increased confidence in the extended algorithm. During our case study, CIVL surfaced a concurrency defect in the algorithm that our random-seed-based tests failed to exercise. This paper shows our workflow, with the goal of helping readers understand its benefits, cost, and tradeoffs compared with testing.

Abstract: Quantum-inspired evolutionary optimization (QIEO) is a new class of population-based metaheuristic optimization algorithms which represents design variables as a set of qubits and searches a continuous, multi-dimensional landscape through rotation of the qubit's amplitude pair. Every generation rotates those amplitudes toward a single elite, which corresponds to that generation's best. The per-generation cost scales as $O(N_p N_g)$ for $N_p$ chromosomes and $N_g$ genes (decision variables). Production use of such solvers is rarely confined to a single machine class. Prototypes are run on laboratory servers- before moving to rented cloud workstations for more involved campaigns. The largest problems are reserved for leadership-class accelerators. This paper asks whether a \emph{single} OpenMP~5 source of QIEO, offloaded with \texttt{\#pragma omp target}, is a viable production path in each of those settings. We report three independent, campaigns of the 0/1 knapsack problem against a same-source multi-core Intel CPU baseline. The study comprises approximately 3,000 runs spanning varying chromosome and gene counts, evaluated using both chromosome-level and gene-level offload strategies on the NVIDIA Tesla V100 SXM2, NVIDIA A100 80GB, and AMD Instinct MI300X GPUs. Deployment-specific nuances such as Volta's constant-memory cliffs, Ampere's L2 persistence and \texttt{cp.async}, CDNA~3's Infinity Cache and XCD occupancy are addressed to ensure high performance of these platforms. Results reveal gene-parallel offload achieved geometric-mean speedups of 90$\times$, 136$\times$, and 155$\times$ over a single CPU core on the V100, A100, and MI300X, respectively, and 12$\times$, 17$\times$, and 16.6$\times$ over 72 host threads. Furthermore DetermineElite, the $O(N_p)$ selection of the generation-best chromosome, is found to be better suited to the host than to the device.

Abstract: We study performance portability for low-order, matrix-free finite element kernels, using the example of a vectorial, variable-coefficient PDE operator originating in geophysical models. Written in Kokkos, the kernel is compared on NVIDIA H100, AMD MI250X, AMD MI300A and Intel PVC Max 1550 GPUs. Owing to its low order and to optimizations that reduce the arithmetic, the kernel has a low arithmetic intensity, so that its performance is determined by how the finite element assembly is mapped onto the memory hierarchy. This is a dimension in which the architectures differ even within one vendor family, causing different performance characteristics. We examine how Kokkos' hierarchical parallelism and shared scratch memory, which are used for the shared degrees of freedom of the conforming discretization, behave on each device. Finally, we show how portability gaps can be narrowed with tuning levers such as the size of the thread groups, the balance between occupancy and register use, and the atomic accumulation strategy at the end of the kernel.

Abstract: We present TERRA-NG, a portable, GPU-accelerated, matrix-free mantle-convection code. A single Kokkos C++ implementation runs at scale on NVIDIA, AMD, and Intel GPU supercomputers. TERRA-NG has a deliberately narrow design: built on a radially extruded mesh of spherical wedges, tailored to the spherical shell geometry, which enables domain-specific optimizations like single quadrature-point integral-evaluations, radial coordinate storage compression and radial shared-memory tiling. The corresponding low-order $W_1$-iso-$W_2/W_1$ wedge-based Stokes--energy discretisation is verified against the Zhong et al.(2008) spherical-shell convection benchmark suite. We showcase TERRA-NG through strong- and weak-scaling on the JUWELS Booster (NVIDIA A100), MareNostrum 5 (NVIDIA H100), LUMI-G (AMD MI250X), Hunter (AMD MI300A APU), and SuperMUC-NG Phase 2 (Intel PVC) supercomputers. Coupled mantle convection simulations at $\sim\!11$ km and $\sim\!5.6$ km radial spacing ($\sim 2.8$ B and $\sim 22$ B DoFs) can be run routinely on standard node partitions of all considered systems. Global $\sim\!1$ km-per-gridpoint mantle convection ($\sim 1.4$ T DoFs) is feasible on an extreme-scale allocation, and a sub-km hero-run at $\sim\!0.7$ km grid spacing scaling up to $\sim 11,000$ GPUs of LUMI-G ($\sim 11$ T DoFs) shows the potential of the code on future, larger machines.

Abstract: We introduce TokaGLINT, a GPU-accelerated implicit solver for electromagnetic field computations in full 3D tokamak simulations, aimed at efficient large-scale parallel GPU computing. Its central innovation lies in the co-design of hierarchical domain decomposition and a fast exact local solver, where hierarchical partitioning is tailored to match fine-grained intra-card subdomains and exploit the tensor-based solver dedicated to curvilinear-coordinate symplectic CN-FDTD-discretized 3D Maxwell equations. Backed by automated operator fusion and batching customized for the intra-card multi-subdomain structure, the solver decouples unknowns through discrete transformations and leverages tensor-structured computations to achieve high hardware utilization, while preserving the long-time stability characteristic of symplectic discretizations. TokaGLINT scales the electromagnetic field solve beyond 10,000 GPUs, achieving 90.1% weak and 53.9% strong scaling efficiency, while delivering a 2.67X single-node speedup over an unpreconditioned BiCGStab baseline (HIP-enabled HYPRE). It is validated in EAST tokamak simulations within the SymPIC plasma simulation code, enabling high-fidelity long-duration modeling.

Abstract: The Max-Host Dragonfly+ topology's original paper proves that there is a 2:1 worst-case oversubscription ratio in expectation for the permutation traffic pattern. We show that the proof only covers a subset of permutation patterns, specifically those in which all host pairs are in different groups, and for any receiver group there are at least two sending hosts in distinct groups. The proof contains two mistakes that cancel out to produce the correct result. Furthermore, the proof implicitly uses an important observation without a backing argument: traffic leaving an indirect group, now forced to follow min-cost paths to receiver groups, can still be split across almost all outgoing global links, in contrast with splitting traffic from a sender group that can naturally select any outgoing group. We prove that the 2:1 ratio holds in expectation for a larger set of patterns, including any permutation. We only need to know that any host sends and receives at most line rate traffic. This constraint can be altered to obtain approximations for expected FCT bounds on any pattern. We also attack the oversubscription problem without the expectation assumption, and find bounds that hold with high probability on permutations for small switch radixes. While in expectation the global layer wiring doesn't affect the oversubscription rate, we experimentally find that it matters in the general case. We find topologies that obtain visible speedups against the default global wiring from the original paper, under $7\%$ for radixes at most $8$.