Ilmu Komputer & AI editorial
Open AccessOA2026
Accelerating the Solving of Many Tiny General Linear Systems on GPUs: Application to Constitutive Laws
A specialized LUpp solver for tiny systems achieves up to 17.7x speedup over MAGMA in constitutive law evaluation on NVIDIA H100 GPUs
Tristan Chenaille; Francesca Cuteri; Rapha{ë}l Prat; Guillaume Latu; Thomas Helfer· 2026· DOI 10.48550/arXiv.2609.15217
The core problem
Many scientific and engineering applications require solving large numbers of independent linear systems on GPUs. While GPU-accelerated solvers are well-established for small to large systems, tiny systems—defined here as those with dimension below 32—remain challenging due to their low arithmetic intensity and the overhead of kernel launches. This challenge is particularly acute in the evaluation of constitutive laws, where millions of integration points are processed independently. Each constitutive update typically relies on a Newton iterative method, and each iteration requires the double-precision solution of a tiny general square linear system using LU factorization with partial pivoting (LUpp). The present study was conducted within a closed-source prototype, serving as a demonstrator for porting to NVIDIA GPUs the constitutive law evaluations currently provided on CPUs by TFEL/MFront, an open-source code generation tool for material knowledge. The authors compare several double-precision LUpp solvers, including implementations from GPU linear algebra libraries (cuSolverDx, MAGMA) and custom-designed CUDA kernels, first on large batches of standalone linear systems and then
Innovation
Experiments on an NVIDIA H100 GPU show that the specialized LUpp solvers proposed in this work can outperform existing state-of-the-art approaches for this class of workloads. The speedups achieved are up to 6.5x over cuSolverDx and up to 17.7x over MAGMA. The results indicate that the best LUpp solving strategy strongly depends on several factors, including system size and application context. For standalone linear systems, the custom kernels demonstrate significant advantages for dimensions below 32. Within the full constitutive-law evaluation workflow, the benefits are even more pronounced due to the ability to invoke the solver directly from device code, reducing kernel launch overhead and enabling better integration with the Newton loop. The study also finds that assigning several threads to each system is not always the most efficient strategy; for very tiny systems, a single thread per system may be preferable to avoid synchronization overhead. Register pressure and occupancy are critical for performance, and the custom kernels are tuned to balance these factors. The speedups are consistent across different batch sizes and system dimensions, with the largest gains observed f
Many scientific and engineering applications require solving large numbers of independent linear systems on GPUs. While GPU-accelerated solvers are well-established for small to large systems, tiny systems—defined here as those with dimension below 32—remain challenging due to their low arithmetic intensity and the overhead of kernel launches. This challenge is particularly acute in the evaluation of constitutive laws, where millions of integration points are processed independently. Each constitutive update typically relies on a Newton iterative method, and each iteration requires the double-precision solution of a tiny general square linear system using LU factorization with partial pivoting (LUpp). The present study was conducted within a closed-source prototype, serving as a demonstrator for porting to NVIDIA GPUs the constitutive law evaluations currently provided on CPUs by TFEL/MFront, an open-source code generation tool for material knowledge. The authors compare several double-precision LUpp solvers, including implementations from GPU linear algebra libraries (cuSolverDx, MAGMA) and custom-designed CUDA kernels, first on large batches of standalone linear systems and then within the full constitutive-law evaluation workflow. The study aims to identify the best LUpp solving strategy, which is shown to depend on system size and application context, and to discuss key aspects such as register pressure, occupancy, the ability to invoke the LUpp solver directly from device code, and whether assigning several threads to each system is most efficient.
The study evaluates several double-precision LUpp solvers for tiny linear systems on NVIDIA GPUs. The solvers include implementations from GPU linear algebra libraries (cuSolverDx and MAGMA) and custom-designed CUDA kernels. The comparison is performed in two phases: first on large batches of standalone linear systems, and then within the full constitutive-law evaluation workflow, where each integration point requires a sequence of distinct linear systems, one per iteration of its own Newton loop. The performance of each solver is assessed on an NVIDIA H100 GPU. Key factors analyzed include register pressure, occupancy, the ability to invoke the LUpp solver directly from device code, and the efficiency of assigning multiple threads per system. The custom kernels are designed to optimize these factors for tiny systems. The study also considers the overhead of kernel launches and the impact of batching. The LUpp algorithm solves a linear system by factoring , where is a permutation matrix, is lower triangular, and is upper triangular, then solving and via forward and backward substitution. For tiny systems, the factorization and substitution steps are performed with minimal synchronization. The performance metric is the total time to solution, and speedups are reported relative to cuSolverDx and MAGMA.
Why it matters
The study highlights that solving many tiny linear systems on GPUs requires specialized strategies that differ from those for larger systems. The key factors influencing performance are register pressure, occupancy, the ability to invoke the LUpp solver directly from device code, and the thread-to-system mapping. The custom CUDA kernels are designed to minimize register usage and maximize occupancy while avoiding unnecessary synchronization. The ability to call the solver from device code is crucial in the constitutive law workflow, as it allows the entire Newton iteration to be executed within a single kernel, eliminating intermediate global memory transfers and kernel launch overhead. The results show that for tiny systems, a single thread per system often outperforms multi-thread approaches due to lower synchronization costs. The study also notes that the optimal strategy depends on the specific application context and system size, and that a one-size-fits-all approach is not effective. The findings have implications for other applications requiring many tiny linear solves, such as implicit solvers in computational fluid dynamics and real-time simulation. The authors suggest that further optimizations could be achieved by exploiting the structure of the matrices or by using mixed-precision approaches. The prototype demonstrates the feasibility of porting TFEL/MFront constitutive law evaluations to GPUs, with significant potential for speedup in large-scale simulations.
Who should read this
CS practitioners and researchers
Opening member content…