Symmetric positive-definite system solvers based on Cholesky factorization are fundamental to many scientific applications, such as climate modeling. We present a portable, nested recursive mixed-precision solver designed for Matrix Processing Units (MXUs), including NVIDIA Tensor Cores (H200) and AMD Matrix Cores (MI300X), that assigns low-precision FP16 arithmetic to large off-diagonal blocks, while preserving high precision on diagonal blocks to ensure numerical stability. The solver is implemented in Julia, providing a high-level, hardware-agnostic interface. We demonstrate up to a 5.07 × speedup relative to the diagonal-precision vendor baseline, with 100 × better accuracy than pure half precision on H200, providing higher accuracy than low-precision at higher speed than high-precision. Positive performance trends are also observed on MI300X, demonstrating broad applicability across GPUs.

We introduce a large language model (LLM)-driven approach for generating dynamic algorithmic dispatch heuristics in high-performance linear algebra. By combining prompt engineering with LLaMA 3 and a curated performance database, the model learns to synthesize selection heuristics that exploit structural patterns to identify fast algorithmic choices. A case study on LU factorization demonstrates the model's ability to replicate expert-designed strategies. This work, developed as part of the DARPA-MIT SmartSolve project, highlights the promise of LLMs for algorithmic discovery and the development of more adaptive, fast linear algebra software.

Agent-based epidemic models (ABMs) encode behavioral and policy heterogeneity but are too slow for nightly hospital planning. We develop county-ready surrogates that learn directly from exascale ABM trajectories using Universal Differential Equations (UDEs): mechanistic SEIR-family ODEs with a neural-parameterized contact rate κ_ϕ(u,t) (no additive residual). Our contributions are threefold: we adapt multiple shooting and an observer-based prediction-error method (PEM) to stabilize identification of neural-augmented epidemiological dynamics across intervention-driven regime shifts; we enforce positivity and mass conservation and show the learned contact-rate parameterization yields a well-posed vector field; and we quantify accuracy, calibration, and compute against ABM ensembles and UDE baselines. On a representative ExaEpi scenario, PEM-UDE reduces mean MSE by 77

Geothermal well arrays, which organize multiple geothermal wells into carefully planned geometric configurations, provide opportunities to enhance energy production capacity and increase fault tolerance. The development and adoption of these emerging geothermal technologies could be accelerated through the recent advances in large language models (LLMs) and high-level high-performance languages. A challenge in LLM-based applications is the reliability of the generated outputs, as they can be prone to subjective biases and hallucinations. This study assesses the potential of cutting-edge LLMs - such as ChatGPT, Gemini, Claude, Grok, and domain-specific models like AskGDR - as expert assistants that can synthesize insightful interpretations of complex geothermal data, as well as improve feature capabilities of geothermal models and numerical software. We developed a novel approach, leveraging Google's recently introduced AI assistant, NotebookLM, to accelerate the generation of unpublished quantitative geothermal benchmarks. The rapid generation of these evaluation instruments is essential for assessing the swiftly evolving capabilities of emerging language model technologies. In particular, we use these benchmarks and LLM-based interviews to analyze opportunities and limitations of two promising technologies: geothermal well arrays and closed-loop coaxial wells. Furthermore, we present a case study illustrating how LLMs can facilitate auto-parallelization of geothermal numerical models. Our analysis emphasizes their application in digital twins and underscores the importance of high-level, high-performance code generation. This line of research could play a transformative role in the geothermal sector by enabling the next-generation of decision-support applications, integrating data analysis, informed recommendations, and more dynamic numerical modeling workflows.

Positive-semidefinite matrices are most efficiently factored using the Cholesky decomposition. For indefinite matrices, the Cholesky factorization does not exist, and the alternatives face greater challenges in achieving numeric stability and preservation of banded structure. Here we pursue an analogy between the requirement for positive-semidefinite matrices and the solution of the quadratic equation x^2 = c for c <= 0. It is shown that a non-associative algebra, called the hemiplex numbers, allows the Cholesky factorization to be computed for arbitrary symmetric matrices. Crucially, the hemiplex Cholesky factorization does not require pivoting for its existence or stability, allowing it to preserve banded structure. For singular matrices it produces a parametrization of the null space, and provides opportunity for truncation of nearly-null directions in a manner similar to common usage of the singular value decomposition. The hemiplex Cholesky factorization may be a practically useful addition to the tools for solving symmetric linear equations.

The reduction of a banded matrix to bidiagonal form is a critical step in the calculation of Singular Values, a cornerstone of scientific computing and AI. Although inherently parallel, this step has traditionally been considered unsuitable for GPUs due to its memory-bound nature. However, recent advances in GPU architectures, such as increased L1 memory per Streaming Multiprocessor (SM)/Compute Unit and larger L2 caches, have shifted this paradigm. In this work, we present the first GPU-accelerated algorithm for reducing a banded matrix to bidiagonal form, integrated into an open-source software package. Our algorithm builds on prior multicore CPU cache-efficient bulge-chasing methods, adapted to modern GPU architectures to optimize throughput. Leveraging Julia's high-level array abstractions and KernelAbstractions.jl, we implement a single function that is both hardware-agnostic and data-precision-aware, running efficiently across NVIDIA, AMD, Intel, and Apple Metal GPUs. We develop a hardware-aware performance model to guide tuning and identify key hyperparameters that govern optimal GPU performance for memory-bound workloads. We show that such workloads, when carefully optimized, can achieve substantial speed-ups on modern GPUs: our implementation outperforms multithreaded CPU libraries (PLASMA, SLATE) starting from matrix sizes as small as $1024\times 1024$, and achieves over $\mathbf{100}\times$ speed-up on $32k\times 32k$ matrices. Moreover, the algorithm's performance scales linearly with the matrix bandwidth, enabling efficient reduction of matrices with larger bandwidths, previously considered impractical.

We present an exact sampling algorithm for Pfaffian point processes based on a skew-symmetric analogue of the Cholesky factorization. This algorithm enables efficient sampling of a wide range of statistics arising in random matrix theory and combinatorics. For instance, we can sample eigenvalues of the orthogonal and symplectic ensembles (β= 1,4). In addition, we introduce a symplectic Arnoldi method for computing skew-orthogonal polynomials associated with a general weight function. This method can be used to efficiently construct the 2 × 2 matrix valued skew-symmetric kernels that arise in β= 1,4 polynomial ensembles. We illustrate our approach with several numerical examples and experiments, including the symmetric corner growth model, the finite-N Gaussian (Hermite) orthogonal and symplectic ensembles, and the β= 1,4 Airy point processes and Tracy-Widom distributions.

Partial differential equations (PDEs) are foundational to modeling in science and engineering, but constructing reliable numerical solvers remains labor-intensive, demanding expert knowledge of discretization schemes, stability conditions, and boundary treatments. Recent work has begun to frame PDE solving as a code-generation task for large language models (LLMs), yet existing approaches operate primarily at inference time: relying on prompting, debugging, self-refinement, and test-time scaling rather than adapting the model itself. In parallel, reinforcement learning with verifiable rewards has emerged as a post-training paradigm for code and math reasoning, but its verifiers are typically binary: a compiler runs, or a test passes. Such signals discard the graded structure of scientific correctness, where two solvers may both execute and yet differ in solution accuracy by orders of magnitude. In this work, we introduce RLVP: Reinforcement Learning with Verifiable Physics, an RL post-training framework for multi-PDE solver code generation. RLVP addresses this verifiability gap with a hybrid verifier: hard program-validity checks ensure executability, while continuous physics rewards score function-space accuracy and PDE-residual consistency. A single policy is post-trained across diverse PDE families spanning hyperbolic, parabolic, elliptic, and incompressible-flow systems. RLVP improves over both pre-trained and supervised-only baselines on PDE benchmarks, and shows zero-shot improvement transfer to held-out PDEs. We show that a smaller LLM post-trained with RLVP can outperform prompting a frontier model on in-distribution PDE solver generation. The trained policy shows evidence of compositionality in numerical motifs: it recombines stencils, time-stepping schemes, and boundary-handling primitives learned from the PDEs used in training into generated solvers for unseen PDE problems.

Generative models have emerged as scalable surrogates for physical simulation, yet they offer no guarantee that their outputs respect the conservation laws, boundary conditions, and nonlinear invariants that govern the underlying physics. Constrained sampling closes this gap, enforcing such constraints exactly at inference time without retraining, but at a computational cost: projection, correction, and trajectory-optimization steps are repeated during sampling, with these steps becoming expensive for nonlinear constraints. Standard ML frameworks exacerbate this: their dense tensor algebra and limited sparse solver composability obscure the structure that physical constraints naturally induce, making efficient batched nonlinear optimization difficult to realize in practice. We address this bottleneck by exploiting the structure that sample-wise batching and local PDE couplings induce in the projection subproblems – namely, block-sparse Jacobian and KKT systems – exposing this structure using ExaModels.jl and solving the resulting sparse nonlinear programs with MadNLP.jl and GPU sparse factorization. Applied to Physics-Constrained Flow Matching (PCFM), on PDE benchmarks with linear, nonlinear, one-dimensional, and two-dimensional constraints, this approach accelerates nonlinear constraint projection while maintaining constraint satisfaction. These results show that sparse GPU nonlinear optimization is a practical foundation for constrained generative sampling in scientific machine learning.

Traditionally, GPU kernels have been developed and optimized within vendor-specific programming models to achieve high performance, resulting in software that is difficult to optimize and adapt across increasingly heterogeneous computing systems. Hardware-agnostic programming models offer a more sustainable approach to GPU software development by improving portability and maintainability, but achieving efficient execution across diverse architectures remains challenging. We address this challenge by integrating auto-tuning into hardware-agnostic GPU kernels written in Julia. We rebuild the established Kernel Tuner auto-tuning framework with Julia support, enabling systematic exploration of kernel configurations for hardware-agnostic GPU kernels targeting NVIDIA, AMD, Intel, and Apple GPUs. We demonstrate this approach on hardware-agnostic singular value decomposition (SVD) as implemented in the NextLA.jl linear algebra library. The results show that auto-tuning is essential for creating resource-efficient hardware-agnostic GPU kernels across a variety of hardware. Optimal configurations improve kernel performance by a factor of 3x to 7x compared to median parameter configurations, demonstrating the substantial impact of tuning on efficient hardware utilization.

AI coding agents have quickly become omnipresent in software engineering. Their serial performance, both in terms of accuracy and speed, has been extensively covered. However, recent initial results suggest their parallel programming capabilities lag behind serial programming capabilities. This paper presents a cross-language evaluation of three coding agents – Cursor's Composer 2.0, GPT 5.4, and Claude Sonnet 4.6 – on parallel code generation across three algorithm categories – sorting, graph traversal, and search – in C++, Python, and Julia. For each algorithm and language pair, we prompt a coding agent to produce a parallel implementation from a serial baseline, track the prompting effort required to achieve both functional correctness and performance improvements, and measure speedup against both custom serial baselines and third-party library implementations. We find that coding agents can produce correct parallel implementations with modest prompting effort, but that achieving meaningful speedup is heavily algorithm- and language-dependent. Sonnet 4.6 delivers the strongest overall performance gains, whereas GPT 5.4 produces no measurable speedups despite consistent correctness. C++ is most consistently parallelizable for graph algorithms, while Python and Julia achieve the largest speedups on search algorithms: no single language dominates across all categories. Python and Julia each achieve speedup on some graph algorithms but regress on others. These findings underscore the impact of including runtime performance efficiency as a main LLM performance metric, in addition to accuracy, particularly for parallel implementations.

Efficiently solving nonlinear equations underpins numerous scientific and engineering disciplines, yet scaling these solutions for challenging system models remains a challenge. This article presents NonlinearSolve. jl-a suite of high-performance open source nonlinear equation solvers implemented natively in the Julia programming language. NonlinearSolve.jl distinguishes itself by offering a unified API that accommo-dates a diverse range of solver specifications alongside features such as automatic algorithm selection based on runtime analysis, support for static array kernels for improved GPU computation on smaller problems, and the utilization of sparse automatic differentiation and Jacobian-free Krylov methods for large-scale problem-solving. Through rigorous comparison with established tools such as PETSc SNES, Sundials KINSOL, and MINPACK, NonlinearSolve.jl demonstrates robustness and efficiency, achieving significant ad-vancements in solving nonlinear equations while being implemented in a high-level programming language. The capabilities of NonlinearSolve.jl unlock new potentials in modeling and simulation across various domains, making it a valuable addition to the computational toolkit of researchers and practitioners alike.

The Singular Value Decomposition (SVD) is a foundational building block in many applications, including low-rank adaptation (LoRA) for large language models (LLMs). Historically, separate SVD implementations have been designed for each data precision, for each hardware vendor, and for each hardware type (personal computer and HPC). This divergence leads to increased development time, the need to redevelop entire libraries when new architectures or data types emerge, and significant complexity for the end user. In this abstract, we discuss a work in progress to develop an alternative: a unified SVD, enabled by abstraction layers. We demonstrate that state-of-the-art performance across the board can be reached using abstraction frameworks, and investigate the performance engineering process and the characteristics that enable adaptable performance.

This paper demonstrates new methods and implementations of nonlinear solvers with higher-order of convergence, which is achieved by efficiently computing higher-order derivatives. Instead of computing full derivatives, which could be expensive, we compute directional derivatives with Taylor-mode automatic differentiation. We first implement Householder's method with arbitrary order for one variable, and investigate the trade-off between computational cost and convergence order. We find that the second-order variant, i.e., Halley's method, to be the most valuable, and further generalize Halley's method to systems of nonlinear equations and demonstrate that it can scale efficiently to large-scale problems. We further apply Halley's method on solving large-scale ill-conditioned nonlinear problems, as well as solving nonlinear equations inside stiff ODE solvers, and demonstrate that it could outperform Newton's method.

This paper presents a portable, GPU-accelerated implementation of a QR-based singular value computation algorithm in Julia. The singular value decomposition (SVD) is a fundamental numerical tool in scientific computing and machine learning, providing optimal low-rank matrix approximations. Its importance has increased even more in large-scale machine learning pipelines, including large language models (LLMs), where it enables low-rank adaptation (LoRA). The implemented algorithm is based on the classic two-stage QR reduction, consisting of successive matrix reduction to band form and bidiagonal form. Our implementation leverages Julia's multiple dispatch and metaprogramming capabilities, integrating with the GPUArrays and KernelAbstractions frameworks to provide a unified type and hardware-agnostic function. It supports diverse GPU architectures and data types, and is, to our knowledge, the first GPU-accelerated singular value implementation to support Apple Metal GPUs and half precision. Performance results on multiple GPU backends and data types demonstrate that portability does not require sacrificing performance: the unified function outperforms most linear algebra libraries (MAGMA, SLATE, rocSOLVER, oneMKL) for matrix sizes larger than 1024 x 1024, and achieves 80%-90% of the performance of cuSOLVER for large matrices.

Symbolic computation systems suffer from memory inefficiencies due to redundant storage of structurally identical subexpressions, commonly known as expression swell, which degrades performance in both classical computer algebra and emerging AI-driven mathematical reasoning tools. In this paper, we present the first integration of hash consing into JuliaSymbolics, a high-performance symbolic toolkit in Julia, by employing a global weak-reference hash table that canonicalizes expressions and eliminates duplication. This approach reduces memory consumption and accelerates key operations such as differentiation, simplification, and code generation, while seamlessly integrating with Julia's metaprogramming and just-in-time compilation infrastructure. Benchmark evaluations across different computational domains reveal substantial improvements: symbolic computations are accelerated by up to 3.2 times, memory usage is reduced by up to 2 times, code generation is up to 5 times faster, function compilation up to 10 times faster, and numerical evaluation up to 100 times faster for larger models. While certain workloads with fewer duplicate unknown-variable expressions show more modest gains or even slight overhead in initial computation stages, downstream processing consistently benefits significantly. These findings underscore the importance of hash consing in scaling symbolic computation and pave the way for future work integrating hash consing with e-graphs for enhanced equivalence-aware expression sharing in AI-driven pipelines.

Deep generative models have recently been applied to physical systems governed by partial differential equations (PDEs), offering scalable simulation and uncertainty-aware inference. However, enforcing physical constraints, such as conservation laws (linear and nonlinear) and physical consistencies, remains challenging. Existing methods often rely on soft penalties or architectural biases that fail to guarantee hard constraints. In this work, we propose Physics-Constrained Flow Matching (PCFM), a zero-shot inference framework that enforces arbitrary nonlinear constraints in pretrained flow-based generative models. PCFM continuously guides the sampling process through physics-based corrections applied to intermediate solution states, while remaining aligned with the learned flow and satisfying physical constraints. Empirically, PCFM outperforms both unconstrained and constrained baselines on a range of PDEs, including those with shocks, discontinuities, and sharp features, while ensuring exact constraint satisfaction at the final solution. Our method provides a flexible framework for enforcing hard constraints in both scientific and general-purpose generative models, especially in applications where constraint satisfaction is essential.

The growth factor in Gaussian elimination measures how large the entries of an LU factorization can be relative to the entries of the original matrix. It is a key parameter in error estimates, and one of the most fundamental topics in numerical analysis. We produce an upper bound of n 0.2079 ln n + 0.91 $n^{0.2079 \ln n +0.91}$ for the growth factor in Gaussian elimination with complete pivoting — the first improvement upon Wilkinson's original 1961 bound of 2 n 0.25 ln n + 0.5 $2 \, n ^{0.25\ln n +0.5}$ .

This paper presents a performant and portable recursive implementation of triangular matrix-matrix multiplication (TRMM) and triangular solve (TRSM) operations in Julia for GPUs, which form the backbone of many other linear algebra algorithms. This work is based on an existing recursive implementation for TRMM and TRSM, which restructures the operations to include general matrix-matrix multiplication (GEMM) calls, facilitating better utilization of the GPU memory hierarchy, and reducing latency overhead. The unified implementation in Julia harnesses the language’s multiple-dispatch and metaprogramming capabilities through the existing GPUArrays and KernelAbstractions frameworks, enabling performant hardware-agnostic execution across different GPU architectures. By supporting a consistent API, this implementation allows users to seamlessly switch between different GPU backends. The recursive hardware-agnostic implementation we present achieves performance comparable to vendor-optimized (cuBLAS/rocBLAS) libraries for larger matrix sizes and provides such methods for the first time to Apple Silicion hardware with only a few hundred lines of code, demonstrating the power of unified implementations.