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

Resumen Técnico de HPRMAT: Un Solucionador de R-Matrix de Alto Rendimiento con Aceleración por GPU

Planteamiento del Problema
El método de la matriz R es una piedra angular en la teoría de reacciones nucleares, particularmente para cálculos de dispersión de canales acoplados que involucran descripciones microscópicas de estados de continuo (por ejemplo, Canales Acoplados Discretizados por el Continuo o CDCC). En estos cálculos, la función de onda interna se expande en una base de Lagrange-Legendre, lo que conduce a la solución de sistemas lineales densos y complejos de dimensión $N = n_{ch} \times n_{lag}$, donde $n_{ch}$ es el número de canales y $n_{lag}$ es el número de funciones de base.

Para aplicaciones modernas que involucran fusión de iones pesados o núcleos de halo (por ejemplo, $^{11}$Be, $^{6}$He), las dimensiones de las matrices alcanzan rutinariamente $N \sim 10,000$ a $30,000$. El cuello de botella computacional es la solución de estos sistemas lineales. Las implementaciones tradicionales, como el ampliamente utilizado paquete de Descouvemont, dependen de la inversión explícita de la matriz ($x = A^{-1}b$). Este enfoque es numéricamente menos estable que los solucionadores directos, requiere almacenar la matriz inversa completa y no logra aprovechar las bibliotecas de álgebra lineal optimizadas o los aceleradores de hardware modernos. Además, los códigos de R-matrix de la física nuclear existente generalmente carecen de soporte para la aceleración por GPU o aritmética de precisión mixta, lo que limita la capacidad de realizar cálculos a gran escala en estaciones de trabajo estándar de escritorio.

Metodología
El artículo presenta HPRMAT, un marco de trabajo de solucionador autónomo diseñado para ser un reemplazo directo (drop-in replacement) para las rutinas de álgebra lineal del paquete de Descouvemont, manteniendo la compatía total de la API. HPRMAT aborda el cuello de botella computacional a través de cuatro backends de solución intercambiables:

1. Solucionador Denso (Tipo 1): Implementa la factorización LU directa ($M = P \cdot L \cdot U$) utilizando rutinas LAPACK optimizadas (específicamente ZGESV de OpenBLAS) en lugar de la inversión de la matriz. Esto utiliza multihilo y vectorización SIMD en CPUs.
2. Solucionador de Precisión Mixta (Tipo 2): Explota la disparidad en el rendimiento de punto flotante entre la precisión simple (FP32) y la doble (FP64), particularmente en GPUs de consumo, donde la relación FP32:FP64 puede ser de 64:1. El método realiza la costosa factorización LU de orden $O(N^3)$ en precisión simple y recupera la precisión completa de doble precisión mediante el refinamiento iterativo (calculando el residuo en FP64 y resolviendo para las correcciones).
3. Solucionador Woodbury-Kinetic (Tipo 3): Un algoritmo especializado para potenciales locales que aprovecha la estructura de bloques del Hamiltoniano. Descompone la matriz en términos de energía cinética de bloque diagonal y términos de acoplamiento diagonal fuera de la diagonal. Construye un sistema de complemento de Schur que se resuelve en precisión simple, obteniendo una aceleración de factor constante al utilizar aritmética real para la inversión del bloque cinético.
4. Solucionador Acelerado por GPU (Tipo 4): Utiliza la biblioteca cuSOLVER de NVIDIA. Emplea una estrategia de precisión mixta donde la factorización y las soluciones triangulares ocurren en la GPU en precisión simple. Para mantener la exactitud, ofrece un modo "híbrido de refinamiento en el host": la factorización de precisión simple se ejecuta en la GPU, mientras que el cálculo del residuo en doble precisión y los pasos de refinamiento se realizan en el host de la CPU, evitando el rendimiento limitado de FP64 en las GPUs de consumo.

El marco está implementado en Fortran 90/95 con extensiones de CUDA C. Incluye enlaces (bindings) para C, Python y Julia. El código gestiona automáticamente los mecanismos de retroceso (fallback): si la inicialización de la GPU falla o si una factorización de precisión simple encuentra una matriz singular, revierte a los solucionadores de doble precisión basados en CPU.

Resultados Clave
Se realizaron pruebas de rendimiento (benchmarks) en tres plataformas: un Apple M3 Ultra (solo CPU), una estación de trabajo Intel Xeon con una NVIDIA RTX 3090 (CPU + GPU), y un Intel i9-12900 (solo CPU).

• Rendimiento:

  • GPU vs. Legado: Con $N=25,600$, el solucionador de GPU (Tipo 4) logra una aceleración de 41× respecto al código basado en inversión legado (la referencia de Descouvemont).
  • GPU vs. CPU Optimizado: El solucionador de GPU logra aproximadamente una aceleración de 15× respecto al solucionador directo de CPU optimizado (Tipo 1) en la RTX 3090.
  • Mejoras de CPU: Incluso sin GPUs, los solucionadores directos (Tipos 1–3) proporcionan aceleraciones de 5–7× sobre el código de inversión legado en CPUs multinúcleo modernas debido al cambio de inversión a resolución directa y al uso de OpenBLAS.
  • Escalabilidad: El solucionador de GPU escala como $O(N^3)$. Una sola tarjeta de consumo (RTX 3090, 24 GB) puede manejar matrices de hasta $N \approx 25,600$. Un backend multi-GPU (probado con cuatro RTX 3090) extiende esto a $N=76,800$ al distribuir la factorización de precisión simple mientras se mantiene la matriz de doble precisión en el host.

• Exactitud:

  • Los solucionadores de CPU de doble precisión reproducen los resultados de referencia con precisión de máquina ($\sim 10^{-18}$).
  • El solucionador de precisión simple de GPU (Tipo 4) produce errores en la sección eficaz de $\sim 10^{-3}$, lo cual está bien dentro de las incertidumbres experimentales típicas en física nuclear ($\gtrsim 1\%$).
  • El modo híbrido de refinamiento en el host recupera la exactitud de precisión de máquina ($\sim 10^{-16}$) en el mismo hardware de consumo con solo un aumento modesto en el tiempo de ejecución (2.5–3× el tiempo de precisión simple).
  • La validación contra cinco casos de prueba estándar (incluyendo la dispersión $\alpha + ^{208}$Pb y la dispersión nucleón-nucleón) confirma que todos los solucionadores reproducen los observables de referencia dentro de las incertidumbres experimentales.

Significancia y Reivindicaciones
El artículo afirma que HPRMAT representa el estado del arte actual en el rendimiento de los solucionadores de R-matrix para la física nuclear. Su principal significancia radica en la democratización de la computación de alto rendimiento para la teoría de reacciones nucleares. Al aprovechar la aritmética de precisión mixta y las GPUs de grado de consumo (por ejemplo, RTX 3090/4090), HPRMAT permite a los investigadores realizar cálculos de CDCC y de canales acoplados a gran escala en estaciones de trabajo de escritorio estándar, eliminando la necesidad de las costosas GPUs de centros de datos (por ejemplo, A100/H100) que a menudo son inaccesibles para los grupos universitarios.

El trabajo enfatiza que la transición de la inversión de matrices a la resolución directa, combinada con la aceleración por GPU, hace que problemas previamente intratables (como los cálculos XCDCC extendidos para núcleos de halo) sean accesibles. La compatibilidad de tipo "drop-in" asegura que los códigos de usuario existentes puedan beneficiarse inmediatamente de estas ganancias de rendimiento sin requerir modificaciones en el código fuente, solo el relincado. El artículo posiciona a HPRMAT no solo como un solucionador más rápido, sino como una herramienta práctica que cierra la brecha entre las capacidades de la computación de alto rendimiento y las necesidades rutinarias de la comunidad de física nuclear.