⚛️ nuclear theory

HPRMAT: A high-performance R-matrix solver with GPU acceleration for coupled-channel problems in nuclear physics

This paper introduces HPRMAT, a high-performance, GPU-accelerated R-matrix solver that utilizes mixed-precision arithmetic and optimized linear algebra to achieve significant speedups over legacy codes, thereby enabling large-scale coupled-channel nuclear physics calculations on standard consumer workstations.

Original authors: Jin Lei

Published 2026-09-07

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Jin Lei

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

In the microscopic world of the atomic nucleus, particles do not simply bounce off one another like billiard balls; they interact through a complex web of forces that can bind them together or scatter them apart. To understand these interactions, physicists use a mathematical framework known as the R-matrix method. This approach divides the space around a nucleus into two regions: an inner zone where the particles are tightly packed and interacting, and an outer zone where they move freely. The challenge lies in calculating exactly how the particles behave in that inner zone, where the mathematics becomes incredibly dense and difficult to solve. For decades, researchers have relied on standard computer methods to crack these equations, but as the number of interacting particles and energy states increases, the calculations become so massive that they slow down to a crawl, often requiring access to expensive, specialized supercomputers that are out of reach for many university laboratories.

A researcher at Tongji University has now developed a new tool called HPRMAT that dramatically speeds up these calculations, making high-level nuclear physics accessible on standard desktop computers. The core of this breakthrough is a shift in how the computer solves the massive systems of equations required to describe nuclear scattering. Instead of using the traditional, slow method of inverting a giant matrix—a process akin to trying to solve a puzzle by working backward from the finished picture—the new software uses a more direct approach that cuts through the problem efficiently. By combining this direct method with modern graphics processing units, the kind of powerful chips found in high-end gaming computers, the software achieves a level of speed that was previously thought impossible for this type of work.

The software was tested on a consumer-grade graphics card, specifically an NVIDIA RTX 3090, which is widely available and costs a fraction of the price of the specialized data-center servers typically used for such heavy lifting. In tests involving a matrix size of 25,600, a common scale for complex nuclear reactions, the new GPU-powered solver completed the calculation in just 3.5 seconds. In contrast, the best traditional method running on a powerful multi-core processor took 52 seconds, and the older, legacy code that many physicists still use took over two minutes. This represents a speedup of roughly 15 times compared to the optimized modern CPU method and 41 times compared to the older code. The researchers demonstrated that this speed does not come at the cost of accuracy; the results matched the established reference data to a degree of precision that is well within the limits of experimental error, meaning the physical predictions remain trustworthy.

One of the most significant hurdles the software overcomes is the limitation of consumer graphics cards, which are designed to be incredibly fast at single-precision calculations but are intentionally slowed down for double-precision work, which is usually required for scientific accuracy. The new tool cleverly navigates this by performing the bulk of the heavy lifting in the fast single-precision mode and then using a refinement step to correct any small errors, ensuring the final result is as accurate as if it had been calculated entirely in the slower, more precise mode. This strategy allows a single desktop computer to handle problems that previously required a cluster of supercomputers. For the rare cases where a problem is too large for even the biggest consumer card, the software can distribute the work across multiple cards, pushing the limit of what can be solved on a desktop to matrix sizes of 76,800.

The impact of this work extends beyond raw speed; it democratizes access to high-performance computing in nuclear physics. By making these calculations feasible on standard workstations, the tool removes the barrier of needing to book time on expensive, centralized supercomputing facilities. Researchers can now run large-scale simulations of complex nuclear reactions, such as the breakup of weakly bound nuclei, directly in their own labs. The software is designed to be a seamless replacement for existing tools, allowing scientists to switch to this faster engine without rewriting their own research codes. It supports multiple programming languages, including Fortran, C, Python, and Julia, ensuring it can integrate into the diverse workflows of modern scientific teams.

The developers validated the software against five standard test cases used in the field, ranging from simple elastic scattering to complex interactions involving multiple channels and non-local potentials. In every case, the new solver reproduced the reference results with high fidelity. For the most demanding scenarios involving narrow resonances or deep sub-barrier energies, the software includes a safety mechanism that automatically switches to a fully precise calculation if the faster method encounters difficulties, ensuring that the results remain reliable. This combination of extreme speed, high accuracy, and robust safety features marks a significant step forward, transforming a computational bottleneck into a routine task and opening the door for more researchers to explore the intricate dynamics of the atomic nucleus.

Technical Summary of HPRMAT: A High-Performance R-Matrix Solver with GPU Acceleration

Problem Statement
The R-matrix method is a cornerstone of nuclear reaction theory, particularly for coupled-channel scattering calculations involving microscopic descriptions of continuum states (e.g., Continuum-Discretized Coupled-Channels or CDCC). In these calculations, the internal wave function is expanded in a Lagrange-Legendre basis, leading to the solution of large, dense, complex linear systems of dimension $N = n_{ch} \times n_{lag}$ , where $n_{ch}$ is the number of channels and $n_{lag}$ is the number of basis functions.

For modern applications involving heavy-ion fusion or halo nuclei (e.g., $^{11}$ Be, $^{6}$ He), matrix dimensions routinely reach $N \sim 10,000$ to $30,000$. The computational bottleneck is the solution of these linear systems. Traditional implementations, such as the widely used package by Descouvemont, rely on explicit matrix inversion ( $x = A^{-1}b$ ). This approach is numerically less stable than direct solvers, requires storing the full inverse matrix, and fails to leverage modern optimized linear algebra libraries or hardware accelerators. Furthermore, existing nuclear physics R-matrix codes generally lack support for GPU acceleration or mixed-precision arithmetic, limiting the ability to perform large-scale calculations on standard desktop workstations.

Methodology
The paper presents HPRMAT, a self-contained solver framework designed to be a drop-in replacement for the linear algebra routines in Descouvemont's package while maintaining full API compatibility. HPRMAT addresses the computational bottleneck through four interchangeable solver backends:

Dense Solver (Type 1): Implements direct LU factorization ( $M = P \cdot L \cdot U$ ) using optimized LAPACK routines (specifically ZGESV from OpenBLAS) instead of matrix inversion. This utilizes multi-threading and SIMD vectorization on CPUs.
Mixed-Precision Solver (Type 2): Exploits the disparity in floating-point performance between single (FP32) and double (FP64) precision, particularly on consumer GPUs where the FP32:FP64 ratio can be 64:1. The method performs the expensive $O(N^3)$ LU factorization in single precision and recovers full double-precision accuracy via iterative refinement (computing the residual in FP64 and solving for corrections).
Woodbury-Kinetic Solver (Type 3): A specialized algorithm for local potentials that exploits the block structure of the Hamiltonian. It decomposes the matrix into block-diagonal kinetic energy terms and off-diagonal diagonal coupling terms. It constructs a Schur complement system solved in single precision, gaining a constant-factor speedup by using real arithmetic for the kinetic block inversion.
GPU Accelerated Solver (Type 4): Utilizes NVIDIA's cuSOLVER library. It employs a mixed-precision strategy where the factorization and triangular solves occur on the GPU in single precision. To maintain accuracy, it offers a "host-refined hybrid" mode: the single-precision factorization runs on the GPU, while the double-precision residual calculation and refinement steps are performed on the CPU host, avoiding the throttled FP64 performance of consumer GPUs.

The framework is implemented in Fortran 90/95 with CUDA C extensions. It includes bindings for C, Python, and Julia. The code automatically handles fallback mechanisms: if GPU initialization fails or if a single-precision factorization encounters a singular matrix, it reverts to CPU-based double-precision solvers.

Key Results
Benchmarks were conducted on three platforms: an Apple M3 Ultra (CPU only), an Intel Xeon workstation with an NVIDIA RTX 3090 (CPU + GPU), and an Intel i9-12900 (CPU only).

Performance:
- GPU vs. Legacy: At $N=25,600$ , the GPU solver (Type 4) achieves a 41× speedup over the legacy inversion-based code (Descouvemont's reference).
- GPU vs. Optimized CPU: The GPU solver achieves approximately 15× speedup over the optimized CPU direct solver (Type 1) on the RTX 3090.
- CPU Improvements: Even without GPUs, the direct solvers (Types 1–3) provide 5–7× speedups over the legacy inversion code on modern multi-core CPUs due to the switch from inversion to direct solving and the use of OpenBLAS.
- Scalability: The GPU solver scales as $O(N^3)$ . A single consumer card (RTX 3090, 24 GB) can handle matrices up to $N \approx 25,600$ . A multi-GPU backend (tested on four RTX 3090s) extends this to $N=76,800$ by distributing the single-precision factorization while keeping the double-precision matrix on the host.
Accuracy:
- The double-precision CPU solvers reproduce reference results to machine precision ( $\sim 10^{-18}$ ).
- The GPU single-precision solver (Type 4) yields cross-section errors of $\sim 10^{-3}$ , which is well within typical experimental uncertainties in nuclear physics ( $\gtrsim 1\%$ ).
- The host-refined hybrid mode recovers machine-precision accuracy ( $\sim 10^{-16}$ ) on the same consumer hardware with only a modest increase in runtime (2.5–3× the single-precision time).
- Validation against five standard test cases (including $\alpha + ^{208}$ Pb scattering and nucleon-nucleon scattering) confirms that all solvers reproduce reference observables within experimental uncertainties.

Significance and Claims
The paper claims that HPRMAT represents the current state-of-the-art in R-matrix solver performance for nuclear physics. Its primary significance lies in democratizing high-performance computing for nuclear reaction theory. By leveraging mixed-precision arithmetic and consumer-grade GPUs (e.g., RTX 3090/4090), HPRMAT enables researchers to perform large-scale CDCC and coupled-channel calculations on standard desktop workstations, removing the necessity for expensive data-center GPUs (e.g., A100/H100) which are often inaccessible to university groups.

The work emphasizes that the transition from matrix inversion to direct solving, combined with GPU acceleration, makes previously intractable problems (such as extended XCDCC calculations for halo nuclei) accessible. The "drop-in" compatibility ensures that existing user codes can immediately benefit from these performance gains without requiring source code modifications, only relinking. The paper positions HPRMAT not just as a faster solver, but as a practical tool that bridges the gap between high-performance computing capabilities and the routine needs of the nuclear physics community.

More like this