Randomized sketching is a core primitive in randomized numerical linear algebra. On modern hardware architectures, in particular on GPUs, the performance of sparse sketches is limited by memory traffic and atomic accumulation rather than floating-point throughput. This makes sketching a natural target for mixed precision, provided that low-precision accumulation does not degrade the embedding quality. We study mixed-precision GPU implementations of sparse oblivious subspace embeddings, focusing on a SparseStack generalization of the GPU CountSketch kernel of Higgins et al. SparseStack improves embedding quality relative to CountSketch on coherent inputs, but its additional nonzeros per column increase atomic-update contention and reduce throughput. We therefore implement FP16 SparseStack variants using deterministic round-to-nearest, exact stochastic rounding, and dithered rounding, and compare them with FP32 SparseStack, CountSketch, mixed-precision CountSketch, and FlashSketch. Our main empirical finding is that, for the tested regimes, SparseStack embedding quality is insensitive to the FP16 rounding rule. Deterministic, stochastic, and dithered rounding FP16 SparseStack produce nearly identical subspace distortion and sketch-and-solve least-squares accuracy across incoherent, coherent, and adversarial test problems. The dominant accuracy factor is the sketch distribution rather than the quantization rule: SparseStack variants substantially improve distortion on coherent inputs, while all methods behave similarly on incoherent inputs. Since deterministic rounding has the lowest overhead, it provides the best performance–accuracy tradeoff among the FP16 SparseStack variants.

We analyze a Balancing Domain Decomposition by Constraints (BDDC) preconditioner for the solution of three dimensional composite Discontinuous Galerkin discretizations of reaction-diffusion systems of ordinary and partial differential equations arising in cardiac cell-by-cell models like the Extracellular space, Membrane and Intracellular space (EMI) Model. These microscopic models are essential for the understanding of events in aging and structurally diseased hearts which macroscopic models relying on homogenized descriptions of the cardiac tissue, like Monodomain and Bidomain models, fail to adequately represent. The modeling of each individual cardiac cell results in discontinuous global solutions across cell boundaries, requiring the careful construction of dual and primal spaces for the BDDC preconditioner. We provide a scalable condition number bound for the precondition operator and validate the theoretical results with extensive numerical experiments.

As systems and simulations grow in size and complexity, it is challenging to maintain efficient use of resources and avoid failures. In this scenario, monitoring becomes even more important and mandatory. This paper describes and discusses the benefits of the advanced monitoring and control tool JuMonC, which runs under user control alongside HPC simulations and provides valuable metrics via REST-API. In addition, plugin extensibility allows JuMonC to go a step further and provide computational steering of the simulation itself. To demonstrate the benefits and usability of JuMonC for large-scale simulations, two use cases are described employing nekRS and ICON on JURECA-DC, a supercomputer located at the J & uuml;lich Supercomputing Centre (JSC). Furthermore, a large-scale use case with nekRS on JSC's flagship system JUWELS Booster is described. Finally, the interplay between JuMonC and LLview (a standard monitoring tool for HPC systems) is presented using a simple and secure JuMonC-LLview plugin, which collects performance metrics and enables their analysis in LLview. Overall, the portability and usefulness of JuMonC, together with its low performance impact, make it an important application for both current and future generations of exascale HPC systems.

This work is motivated by symmetric positive definite (SPD) matrices for which the best sparse approximate inverse (SPAI) with the prescribed nonzero pattern of A^k for some moderate value of k, e.g., 1, 2, 3, or 4, fails to capture essential features of the inverse, such as definiteness. In this context, we consider short-recurrence iterative methods for the computation of SPAIs, that is, methods with sparse matrix iterates whose nonzero structure is globally updated at each iteration based on short-recurrence relations. In particular, we consider the minimal residual (MR) method, its newly proposed variant enriched with one previous search direction, namely the locally optimal minimal residual (LOMR) method, and the conjugate gradient (CG) method with sparse matrix iterates and Frobenius inner products. We show that, for SPD matrices, the MR method converges linearly with rate (1-λ_min^2/tr(A^2))^1/2, irrespective of the initial guess. While LOMR inherits unconditional monotone convergence from MR, its observed convergence behavior is that of a monotonically decreasing lower envelope to the CG residual norm without the occasional spurious oscillations proper to CG. All three methods are implemented with practical dropping strategies to control the growth of nonzero patterns in the approximate inverse. Numerical experiments are performed where SPAIs are computed with a prescribed cap on density, and the performance of those SPAIs as preconditioners for CG solves is quantitatively assessed for each method.

Modern high-performance computing (HPC) increasingly relies on GPUs, but integrating GPU acceleration into complex scientific frameworks like OpenFOAM remains a challenge. Existing approaches either fully refactor the codebase or use plugin-based GPU solvers, each facing trade-offs between performance and development effort. In this work, we address the limitations of plugin-based GPU acceleration in OpenFOAM by proposing a repartitioning strategy that better balances CPU matrix assembly and GPU-based linear solves. We present a detailed computational model, describe a novel matrix repartitioning and update procedure, and evaluate its performance on large-scale CFD simulations. Our results show that the proposed method significantly mitigates oversubscription issues, improving solver performance and resource utilization in heterogeneous CPU-GPU environments.

Multigrid methods are popular iterative methods for solving large-scale sparse systems of linear equations. We present a mixed precision formulation of the multigrid V-cycle with general assumptions on the finite precision errors coming from the application of coarsest-level solver and smoothing. Inspired by existing analysis, we derive a bound on the relative finite precision error of the V-cycle which gives insight into how the finite precision errors from the individual components of the method may affect the overall finite precision error. We use the result to study V-cycle methods with smoothing based on incomplete Cholesky factorization. The results imply that in certain settings the precisions used for applying the IC smoothing can be significantly lower than the precision used for computing the residual, restriction, prolongation and correction on the concrete level. We perform numerical experiments using simulated floating point arithmetic with the MATLAB Advanpix toolbox as well as experiments computed on GPUs using the Ginkgo library. The experiments illustrate the theoretical findings and show that in the considered settings the IC smoothing can be applied in relatively low precisions, resulting in significant speedups (up to 1.43x) and energy savings (down to 71

Electrostatics simulations often employ Boundary Element Methods (BEM) that result in fully populated (dense) matrices. Direct solvers are typically used for solving of these dense linear systems. But iterative solvers can also be very effective when tunable accuracy is required. In this paper, we investigate to which extent offloading the iterative solver to a GPU accelerator can speed up the overall simulation. We use the Ginkgo library [5] as a solver backend for the Fortran-based CASOPT [17] simulations. We consider three production-relevant test cases and demonstrate the superiority of the GPU-accelerated CASOPT version.

The Joint Laboratory on Extreme-Scale Computing (JLESC) was initiated at the same time lossy compression for scientific data became an important topic for the scientific communities. The teams involved in the JLESC played and are still playing an important role in developing the research, techniques, methods, and technologies making lossy compression for scientific data a key tool for scientists and engineers. In this paper, we present the evolution of lossy compression for scientific data from 2015, describing the situation before the JLESC started, the evolution of this discipline in the past 8 years (until 2023) through the prism of the JLESC collaborations on this topic and some of the remaining open research questions.

The direct solution of batches of band linear systems in parallel is important for many applications. In this paper, we elaborate on three new GPU algorithms for the data-parallel direct solution of linear system batches, sharing a band structure. We develop algorithms for three matrix types: tridiagonal, small bandwidth, and wide bandwidth. We exploit the band structure of the matrix, and store it in an efficient fashion (LAPACK band matrix format) and ensure that the SIMD parallelism of the GPUs are maximized. We develop a panel-based factorization for wide-band matrices to ensure coalesced access (with column-major storage) while minimizing main memory traffic. For the tridiagonal solvers, to ensure a high level of concurrency, we adapt a divide-and-conquer approach and utilize co-operative group functionality to efficiently communicate between compute units (in registers) on GPUs. We implement these algorithms for NVIDIA GPUs and study the performance for varying matrix sizes (16 to 1024) and across a range of batch items (upto 1 × 10 6 ). We compare the performance of our implementations with the corresponding optimized vendor implementations (cuSPARSE and MKL), with the state-of-the-art GPU library MAGMA, and with the optimized LAPACK implementation provided by Intel MKL on Intel Skylake CPUs. We also showcase the effectiveness of our batched band solvers for matrices originating from XGC, a gyrokinetic Particle-In-Cell (PIC) application optimized for modeling the edge region plasma within a plasma physics application. We show that our implementations are on average ∼ 2× (for batched banded solvers, compared to MAGMA and MKL) to ∼ 3× (for batched tridiagonal solvers, compared to cuSPARSE) faster than the state-of-the-art and the vendor provided implementations.

Sparse linear algebra is a cornerstone of many scientific computing and machine learning applications. Python has become a popular choice for these applications due to its simplicity and ease of use. Yet high performance sparse kernels in Python remain limited in functionality, especially on modern CPU and GPU architectures. We present pyGinkgo, a lightweight and Pythonic interface to the Ginkgo library, offering high-performance sparse linear algebra support with platform portability across CUDA, HIP, and OpenMP backends. pyGinkgo bridges the gap between high-performance C++ backends and Python usability by exposing Ginkgo's capabilities via Pybind11 and a NumPy and PyTorch compatible interface. We benchmark pyGinkgo's performance against state-of-the-art Python libraries including SciPy, CuPy, PyTorch, and TensorFlow. Results across hardware from different vendors demonstrate that pyGinkgo consistently outperforms existing Python tools in both sparse matrix vector (SpMV) product and iterative solver performance, while maintaining performance parity with native Ginkgo C++ code. Our work positions pyGinkgo as a compelling backend for sparse machine learning models and scientific workflows.

The performance of the GMRES iterative solver on GPUs is limited by the GPU main memory bandwidth. Compressed Basis GMRES outperforms GMRES by storing the Krylov basis in low precision, thereby reducing the memory access. An open question is whether compression techniques that are more sophisticated than casting to low precision can enable large runtime savings while preserving the accuracy of the final results. This paper presents the lightweight in-register compressor FRSZ2 that can decompress at the bandwidth speed of a modern NVIDIA H100 GPU. In an experimental evaluation, we demonstrate using FRSZ2 instead of low precision for compression of the Krylov basis can bring larger runtime benefits without impacting final accuracy.

Large-scale simulations require efficient computation across the entire computing hierarchy. A challenge of the Exascale Computing Project (ECP) was to reconcile highly heterogeneous hardware with the myriad of applications that were required to run on these supercomputers. Mathematical software forms the backbone of almost all scientific applications, providing efficient abstractions and operations that are crucial to harness the performance of computing systems. Ginkgo is one such mathematical software library, nurtured by ECP, providing high-performance, user-friendly, and performance portable interfaces for applications in ECP and beyond. In this paper, we elaborate on Ginkgo's philosophy of high-performance software that is sustainable, reproducible, and easy to use. We showcase the wide feature set of solvers and preconditioners available in Ginkgo and the central concepts involved in their design. We elaborate on four different ECP software integrations: MFEM, PeleLM + SUNDIALS, XGC, and ExaSGD that use Ginkgo to accelerate their science runs. Performance studies of different problems from these applications highlight the effectiveness of Ginkgo and the benefits incurred by these ECP applications.

In order to understand cardiac arrhythmia, computer models for electrophysiology are essential. In the EuroHPC MicroCARD project, we adapt the current models and leverage modern computing resources to model diseased hearts and their microstructure accurately. Towards this objective, we develop a portable, highly efficient, and performing BDDC preconditioner and solver implementation, demonstrating scalability with over 90% efficiency on up to 100 GPUs.

style to avoid excess iterations for "easier" problems. Ginkgo integrates a simple but customizable stopping criterion for the residual norm and generally uses a pre-defined (relative or absolute) residual norm as the stopping criterion. In order to avoid the overhead of launching a kernel at every iteration, the iteration convergence and iteration control is part of the solver kernel. Each thread maintains its own copy of the iteration count.

The bulk-synchronous nature of the state-of-the-art iterative algorithms form a hurdle to their scalability. With increasing problem size and the need for hundreds of thousands of computing units, alleviating synchronization bottlenecks is important. Asynchronous methods, which remove these synchronization bottlenecks, and compute on latest available, possibly stale data are an alternative to the synchronous methods. The analysis of these asynchronous methods is particularly challenging due to the stochastic nature of their computation. We present a probabilistic model which can he used to understand the behaviour of these asynchronous methods. Using this model, we analyze two asynchronous iterative methods: asynchronous Richardson and the asynchronous multi-level Schwarz methods, and study the convergence behaviour of these methods from an empirical perspective. We also show that tunable parameters such as the relaxation parameter, w and additive weights in the asynchronous Richardson and asynchronous Schwarz methods, respectively, can be used to improve convergence of the methods, depending on the latency and network characteristics of the machine.