Accelerating the Solving of Many Tiny General Linear Systems on GPUs: Application to Constitutive Laws
arXiv:2609.15217v1 [cs.DC] 14 Sep 2026
Tristan Chenaille1,2 , Francesca Cuteri1 , Raphaël Prat1 , Guillaume Latu1 , and Thomas Helfer1 CEA, DES, IRESNE, DEC, SESC, Cadarache, St-Paul-lez-Durance, France Aix-Marseille University, Mathematics and Computer Science Doctoral School, France [email protected]
1 2
Abstract. Many applications require solving large numbers of independent linear systems on GPUs. While this need is well addressed for small to large systems, tiny ones, understood here as systems of dimension below 32, remain challenging. This is especially relevant in constitutive law evaluation, where millions of integration points are handled independently, and where each constitutive update generally relies on a Newton iterative method. Each iteration then 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, which serves as a demonstrator for porting to NVIDIA GPUs constitutive law evaluations currently provided on CPUs by TFEL/MFront, an open-source code generation tool for material knowledge. We compare several double-precision LUpp solvers, including implementations from GPU linear algebra libraries as well as custom-designed CUDA kernels. The comparison is performed 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 study shows that the best LUpp solving strategy strongly depends on several factors including system size and application context. We discuss several key aspects, including register pressure, occupancy, the ability to invoke the LUpp solver directly from device code, and whether assigning several threads to each system is the most efficient strategy. 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, with speedups of up to 6.5× over cuSolverDx, and up to 17.7× over MAGMA. Keywords: GPU Computing · Computational Solid Mechanics · Constitutive Law Evaluation · Newton-Raphson Iterative Method · LU Factorization with Partial Pivoting · Batched Linear Solvers · Kernel Fusion
2
1
T. Chenaille et al.
Introduction
Graphics Processing Units (GPUs) have become a key component of modern high-performance computing systems due to their ability to accelerate many workloads. This trend has shaped the exascale era, with the main exascale systems relying heavily on GPU-based acceleration [1]. As a result, porting traditional CPU-oriented scientific applications to GPUs is now essential. Achieving high performance, however, is not straightforward and often requires substantial algorithmic and data structures rework. Attainable performance also depends on whether the underlying problem can expose enough fine-grained parallelism. Applications involving large numbers of independent computations, including the solution of many linear systems, are typically well suited to GPU execution. GPU-based linear algebra routines have been extensively developed, and efficient implementations are now available for a wide range of problem sizes. In particular, many solvers have demonstrated high performance for small to large linear systems. However, a different regime arises when considering tiny systems, i.e. those of dimension n ≤ 32. The amount of computation within these systems is limited, which reduces the available parallelism and makes it difficult to fully utilize GPU capabilities. Thus, overheads and implementation details that are negligible for larger systems can become dominant, and strategies that perform well at larger scales do not necessarily remain efficient [2]. This specialized regime calls for dedicated optimization strategies, especially in applications where such systems must be solved in large numbers. Such workloads arise in several fields, including reactive flow simulations [3], block-Jacobi preconditioning [4], and computational solid mechanics. In this work, we focus on the latter, and more specifically on evaluating stress-strain constitutive models at integration points. This task is currently performed on CPUs using the TFEL/MFront library [5], which compiles a high-level description of the constitutive law. To load the latter, global simulation codes use the MFrontGenericInterfaceSupport (MGIS) library [6], which drives the behavior evaluation. A representative use case is the high-performance simulation of nuclear fuel elements within the PLEIADES framework [7]: in the MFEM-MGIS application [8], to model the response of fuel under irradiation, phenomenas such as viscoplasticity, swelling, and damage are captured by constitutive models. In large-scale simulations, the material response must be computed at a vast number of integration points, each processed independently. At each point, it is typically obtained by solving a nonlinear problem using a Newton iterative method. Each iteration requires the solution of a general linear system using LUpp in double precision. Although each individual system is small to tiny (usually n ≤ 20), the very large number of independent integration points can make the overall cost significant, which makes this workload a natural candidate for GPU acceleration. It can also avoid costly host-device data transfers: when the surrounding global simulation runs on the GPU, keeping the constitutive evaluation on the CPU would require moving integration-point data between host and device at every iteration. Porting constitutive model evaluation to
Solving of Many Tiny General Linear Systems on GPUs
3
GPUs is therefore a step toward faster simulations, better hardware utilization, and a reduced energy footprint. In this work, we investigate the efficient solution of tiny general linear systems using LUpp on GPUs, with a particular focus on application-driven workloads arising from constitutive law evaluation. We study behavior integration as a batch of integration points on the GPU, decoupled from the surrounding global solution procedure. This focuses the comparison on the LUpp solving strategy and its interaction with the constitutive update, rather than on orchestration choices made by the global solver. When behavior integration is instead fused into a global assembly loop, the conclusions may differ. Our contributions are threefold. First, we design specialized LUpp solvers tailored to tiny systems, exploring different implementation strategies. These solvers are device-callable, meaning that they can be invoked directly from device code. Second, we provide a comprehensive comparison of these solvers against existing implementations from external libraries, both on large batches of standalone linear systems and within the full constitutive law evaluation workflow. Third, experimental results obtained on an NVIDIA H100 GPU demonstrate that the optimal approach strongly depends on both the system size and the application context, and that the proposed specialized LUpp solvers can outperform existing state-of-the-art ones by up to 17.7× in specific cases.
2
Background and Related Work
In this section, we briefly recall principles of Newton–Raphson method and LUbased solving, both of which are central to constitutive law evaluation under implicit schemes. We then review existing approaches for solving tiny linear systems on GPUs, and their adoption in constitutive modeling libraries. 2.1
Newton-Raphson Iterative Method
The Newton–Raphson method is a standard approach for solving nonlinear equations. Under standard assumptions, the method exhibits quadratic convergence when the iterate is sufficiently close to the solution [9]. In the vectorial case, given a nonlinear system F(x) = 0,
with x ∈ Rn and F : Rn → Rn ,
(1)
At iteration k, the method solves (typically using LUpp) the linearized system J(x(k) ) ∆x(k) = −F(x(k) ),
(2)
where J(x) = ∂F/∂x is the square, general Jacobian matrix of F, ∆x(k) is the correction, and F(x(k) ) is the residual. The solution is then updated as x(k+1) = x(k) + ∆x(k) .
(3)
Iterations stop when a prescribed convergence criterion is met, typically based on the residual norm, the correction norm, or a maximum number of iterations.
4
2.2
T. Chenaille et al.
LU Factorization, Triangular Solves, and Variants
LU factorization, also referred to as LU decomposition, consists in factorizing a square, general matrix A as the product of a lower triangular matrix L and an upper triangular matrix U . This is achieved through Gaussian elimination, which proceeds column by column: at each step, the current pivot row is used to eliminate entries below it, producing a reduced submatrix known as the trailing submatrix or Schur complement, on which the process is applied recursively. LU factorization is widely used to solve linear systems of the form Ax = b, since, once the LU factorization of A is available, the solution can be obtained by solving two simple triangular systems successively, namely Ly = b and U x = y. Numerical stability and robustness can be respectively improved and guaranteed by performing row interchanges during Gaussian elimination so that the pivot element is nonzero and as large as possible in magnitude. This strategy is known as partial pivoting. With partial pivoting, the factorization is written P A = LU , where P is a permutation matrix encoding row interchanges. This pivoting technique is widely used; while guaranteeing robustness, it offers a good compromise between numerical stability and computational cost [10]. The factorization step dominates the computational cost, requiring O(n3 ) floating-point operations (flops). Partial pivoting introduces an additional O(n2 ) comparison cost. The two triangular solves require O(n2 ) flops. From an algorithmic point of view, LU factorization is sequential across elimination steps, since each step depends on the previous ones. Parallelism can therefore only be exploited within the operations associated with each step. Different variants have been proposed to organize these operations [11], and can be grouped in two families. In right-looking (RL) variants, the transformations computed at a given elimination step are immediately applied to the trailing submatrix. This exposes substantial matrix-matrix parallelism and maps well to BLAS-3 operations. In left-looking (LL) variants, the current step is instead updated by reusing transformations computed during previous steps. Although they generally expose less immediate matrix-matrix parallelism than RL variants, they can reduce repeated accesses to matrix entries. Orthogonally to this scheduling choice, blocked variants partition the matrix to exploit data locality. The matrix partitioning is often one-dimensional, grouping consecutive columns into panels, but can also be two-dimensional, dividing the matrix into a grid of tiles [12]. All these design choices have a strong impact on data locality, available parallelism, and ultimately performance on different architectures.
2.3
State-of-the-Art GPU Optimization Strategies for Tiny LUpp
The scientific community has studied the efficient GPU implementation of batched partial-pivoting LU (LUpp), from which key optimization principles for the tiny regime (n ≤ 32) can be identified. These principles are largely shared with QR and Cholesky [2]. A central observation is that the factorization is largely memory-bound: data transfer costs dominate the O(n3 ) arithmetic [13,2].
Solving of Many Tiny General Linear Systems on GPUs
5
The classical strategy for large matrices relies on blocked algorithms that partition the matrix into panels and cast the trailing-matrix update as Level-3 BLAS to approach peak compute throughput. For tiny sizes, this approach becomes ineffective: the updates are too small to amortize their overhead [14]. Optimizing performance in this regime therefore requires size-aware kernel designs that minimize data access and movement. The thread-to-system mapping, i.e. how GPU threads are organized to process each independent system, is a central design choice. In the one-thread-persystem approach, a single thread processes an entire matrix typically held in its private register file, which eliminates all synchronization. This paradigm is generally used only for the smallest sizes of the tiny regime (n ≤ 8): each thread must store all n2 elements plus many temporaries, so local memory usage grows rapidly with n and degrades performance [15,16]. In the multi-thread approach, a group of threads cooperates on one matrix, each thread owning one row [2,13]. By spreading the matrix across threads, this design lowers per-thread register pressure, which usually avoids using local memory and increases multiprocessor occupancy [17]. Two storage strategies then arise. In the first one, the entire matrix is stored in shared memory, so every arithmetic operation depends on a shared-memory access. In the second one, register-blocking is used: each thread permanently holds its row in registers [13], with inter-thread communication relying on warp shuffles or on a small shared memory workspace. This significantly reduces data traffic. The LL and RL variants are then preferable in different LUpp cases: with the matrix in shared memory, the LL one reduces the number of accesses and the on-chip data movement [13], whereas with register-blocking, the RL one exposes more per-step parallelism and simplifies the LU’s partial pivoting. Several other optimizations are mentioned in the literature. Dedicated data layouts preserve coalesced memory accesses [13,18]. Tunable concurrency assigns multiple factorizations to a single thread group, increasing per-block workload [2,19]. Kernel fusion combines the factorization, triangular solves, and potentially application-level operations into a single kernel, eliminating redundant global memory round-trips [14,19]. Compile-time specialization via C++ templates enables loop unrolling, optimized register allocation, and aggressive instruction scheduling by the compiler [17,2]. Finally, specific to LUpp, logical pivoting records pivot indices in a vector and avoids the physical row interchange at each elimination step [2,19]. Our prototype builds on all these principles. 2.4
Existing GPU Libraries for Tiny LUpp Solvers
Leveraging some of these insights, a number of libraries implement LUpp on NVIDIA GPUs for solving large numbers of independent linear systems. These implementations can be divided into two categories. On the one hand, batched solvers, such as those provided in MAGMA [20] and cuBLAS [21], operate on large sets of systems but are typically invoked from the CPU. This forces the linear solve to be executed as a separate GPU phase, rather than being fused with application computations, which introduces additional overheads [22]. On the
6
T. Chenaille et al.
other hand, libraries such as cuSolverDx [23], KokkosKernels [24], and SNLS [25] provide device-callable solvers that can be invoked directly from within GPU kernels. This approach avoids host-device synchronization and enables tighter integration within application-driven GPU workflows. Within GPU-based constitutive law evaluation, whether provided by external libraries or developed in-house, linear solvers must be considered within a broader pipeline involving application-level constraints such as Newton iterations and constitutive updates. Figure 1 shows how the batched and device-callable approaches shape this workflow: the former leads to a CPU-driven evaluation, the latter to a fully GPU-driven one. The CPU-driven scheme breaks each Newton iteration into several kernel launches. This adds launch overheads, moves intermediate states through global GPU memory, and requires additional synchronization for convergence checks, which transfer only a single flag. The GPUdriven scheme avoids these costs by keeping the Newton loop inside a single GPU kernel. The next subsection reviews how these two paradigms are adopted in existing constitutive modeling libraries. Algorithm 1 Batched LUpp, CPU-driven constitutive evaluation send inputs: CPU → GPU while not all points converged do launch GPU kernel for all points i in parallel do assemble Fi , Ji ▷ Law launch batched GPU kernel(s) for all points i in parallel do solve Ji ∆xi = −Fi ▷ LUpp launch GPU kernel for all points i in parallel do xi ← xi + ∆xi check Newton convergence send convergence: GPU → CPU send outputs: GPU → CPU
Algorithm 2 Device-callable LUpp, GPU-driven constitutive evaluation send inputs: CPU → GPU launch GPU kernel for all points i in parallel do while point not converged do assemble Fi , Ji ▷ Law solve Ji ∆xi = −Fi ▷ LUpp xi ← xi + ∆xi check Newton convergence send outputs: GPU → CPU
Fig. 1: Simplified comparison of GPU-based constitutive law evaluation with a batched versus a device-callable LUpp solver (no tangent operator computation). 2.5
GPU-based Implicit Constitutive Law Evaluation
In computational mechanics, and in particular in constitutive law integration, implicit schemes based on Newton iterative methods are widely used for their numerical efficiency and their consistency with the global simulation procedure. Locally, they provide a robust framework for solving nonlinear constitutive equations, while their consistent linearization supplies the tangent operator needed by the global equilibrium solver [26]. In contrast, explicit schemes avoid such local iterative solves but may be subject to stronger stability constraints and are not always suitable in the same settings [27]. As a consequence, implicit
Solving of Many Tiny General Linear Systems on GPUs
7
constitutive updates typically require the repeated LUpp solution of small to tiny linear systems involving the Jacobian of the local nonlinear constitutive equations, until a prescribed convergence criterion is reached. Since convergence behavior depends on the local state, different integration points may require different numbers of Newton iterations, which can lead to irregular computational patterns in a parallel workflow. In some cases, this can result in significant load imbalance, making efficient GPU execution delicate [28]. On the batched side, several libraries and frameworks rely on a formulation in which each Newton step is applied uniformly to all integration points. This is notably the case for NEML2 [29], JAX-CPFEM [30], and the FEniCSx externaloperator framework together with its associated libraries [31]. MOOSE can also access this type of GPU path through its NEML2 integration [32]. These works cover a broad range of implicit constitutive models, including viscoplasticity, poroplasticity, and crystal plasticity. A structural cost of this batched paradigm is that integration points that have already converged still participate in subsequent Newton iterations until the last point in the batch reaches convergence. This cost is hard to quantify, as it highly depends on the input data and the constitutive law. Still, the number of Newton iterations across integration points generally follows a light-tailed distribution, which leads the libraries above to consider this cost to be acceptable to pay. On the device-callable side, other implementations keep the Newton loop local to each integration point and invoke the linear solver directly from within the GPU kernel. This strategy is adopted, for example, in ExaCMech [33], which relies on SNLS [25], and in AutoMat [34]. It has also been explored within Fierro [35]. MOOSE has also recently initiated a native Kokkos-based effort along similar lines [36]. Such applications span several classes of implicit constitutive laws, notably viscoplastic and crystal plasticity models. A drawback of this strategy is that it may introduce warp-level divergence. In the batched paradigm, the slowdown is driven by the slowest integration point in the entire batch. Here, the effect is confined to the warp level. It is therefore less penalizing. Taken together, and to the best of our knowledge, these works indicate that there is no consensus on a single dominant strategy for GPU-based implicit constitutive law evaluation. We therefore consider both batched and devicecallable LUpp solvers, described in Section 3 and compared in Section 4.
3
GPU LUpp Solvers and Integration Strategies
In this section, we describe the LUpp solver families implemented in our prototype. Some come from external libraries, while others are developed in-house. We classify them into three categories: batched solvers, multi-threaded devicecallable solvers, and single-threaded device-callable solvers. 3.1
Batched Solvers
At present, two libraries provide batched LUpp solvers for NVIDIA GPUs: cuBLAS, NVIDIA’s proprietary library, and MAGMA, which is open source. All these routines assume column-major storage and are multi-threaded: each
8
T. Chenaille et al.
system of dimension n ≤ 32 is processed by n active GPU threads, one per row, enabling coalesced accesses within each column. • cuBLAS: cuBLAS exposes two separate routines: cublasDgetrfBatched, which performs a RL LUpp factorization, and cublasDgetrsBatched, which performs the subsequent triangular solves. The three operands all reside in shared memory: matrix, right-hand side (RHS), and pivot vector. • MAGMA: In this work, MAGMA refers to magma_dgesv_batched, a unified routine combining LUpp factorization and triangular solves. Although MAGMA also provides magma_dgetrf_batched and magma_dgetrs_batched, the unified magma_dgesv_batched is more efficient throughout the tiny regime, dispatching for n ≤ 32 to a specialized RL, logical-pivoting, register-blocking kernel, named dgesv_batched_small_kernel. It uses register-blocking on the matrix and RHS, with the pivot vector in shared memory. 3.2
Device-callable Solvers
Multi-threaded Solvers In the multi-threaded approach, a group of threads cooperates on a single system. This follows the conventional GPU rationale: a single thread’s register limit is too low to hold a system without using local memory. The work is distributed across several cooperating threads instead. This lowers per-thread register pressure and can raise occupancy. The matrix is spread over these threads, held in their combined registers or in shared memory, and the threads synchronize to advance the factorization together. • MAGMA-custom (in-house): This solver reimplements MAGMA’s internal dgesv_batched_small_kernel. The numerical LUpp solve itself is unchanged: column-major storage, RL arithmetic, logical pivoting, and register-blocking are all kept as is. Only the thread-to-data mapping differs, so numerical results and stability are identical. The reimplementation extends the kernel in three ways. First, it is device-callable. MAGMA already implements this solve on the device, but the routine is internal to the compiled library and is not declared in any public header. Only its host-side launcher is exposed, so the solve cannot be called from a user kernel. Our version is callable from device code, so it can be fused inside a single Newton kernel. Second, it packs several systems per block. MAGMA assigns one system per block with n threads. When n < 32, this leaves 32 − n lanes of the warp idle. We launch full warps instead and fit as many systems as possible. For n = 12, a warp solves two systems and wastes only 8 lanes. Third, the number of rows per thread K is a tunable template parameter (1 to 6). MAGMA fixes K = 1, so a system of size n needs n threads. With K rows per thread a system needs T = ⌈n/K⌉ threads, and a warp holds ⌊32/T ⌋ systems. For n = 12 and K = 2, six threads solve one system, five systems fit in a warp, and only two lanes stay idle. When n is not a multiple of K, a phantom mechanism absorbs the remainder. The tail thread owns dead rows that are skipped at compile time. This adds no computational cost, and no register pressure, since every thread allocates the
Solving of Many Tiny General Linear Systems on GPUs
9
same K slots. Coalescing is preserved: rows are assigned round-robin, so within a system consecutive threads read consecutive memory along each column. • cuSolverDx-block: This solver is the RL gesv_partial_pivot routine of cuSolverDx, a proprietary device-callable NVIDIA library. It runs in the Block() execution mode, in which the threads of a block cooperate on the systems assigned to that block. All operands reside in shared memory and can benefit from coalesced accesses. The SystemsPerWarp parameter sets how many systems a warp processes. As with cuBLAS, the internal implementation remains opaque. Single-threaded Solvers In the single-threaded approach, each GPU thread solves an entire system on its own. The natural design keeps the whole matrix in the thread’s registers, but, as discussed in Section 2.3, this is only viable for small n (n ≲ 8). The following solvers therefore let their operands be placed in global or shared memory, used here not as a cooperation workspace, as in the multithreaded approach, but as a private extension of each thread’s register space. This reduces the unintended use of local memory. However, for large n (n ≳ 30 on GPUs with a 227KB per-block shared-memory limit such as the H100), the shared variant can force the block size below 32 to fit the shared-memory budget, leaving idle warp lanes and lowering effective thread-level occupancy. • cuSolverDx-thread: This solver is the RL gesv_partial_pivot routine of cuSolverDx running in the Thread() execution mode. All three operands can be placed in registers, shared memory, or global memory. In shared or global memory, only the RHS can benefit from coalesced accesses. • KokkosKernels: This solver processes each system using two open-source KokkosKernels functions: SerialGetrf and SerialGetrs. The factorization is a recursive RL, physical-pivoting LU, like LAPACK’s getrf2. Since device code cannot recurse cheaply, recursion is emulated by a fixed-size per-thread software stack. Passing the operands as strided Kokkos views lets each reside in registers, shared, or global memory, enabling coalesced access in the latter two cases. • SNLS: This solver uses SNLS’s open-source SNLS_LUP_Solve routine, a simple RL, logical-pivoting LUpp. While the pivot vector must stay in registers, matrix and RHS can be placed in registers, shared memory, or global memory. In shared or global memory, neither can benefit from coalesced accesses. • TFEL-baseline: This solver is used as the baseline in this work because it corresponds to the open-source tfel::math::TinyMatrixSolve routine, already used in production through the TFEL/MFront/MGIS framework. This routine was originally designed for CPU execution, where it is competitive with reference implementations. Its GPU port used here is direct and includes no GPU-specific optimization. It is a LL, logical-pivoting LUpp solver. All operands can be placed in registers, shared memory, or global memory, enabling coalesced access in the latter two cases through tfel::math::StridedCoalescedView objects. • Tiled (in-house): These solvers come in RL and LL variants. They split the matrix into square tiles of tunable size (≤ 6 to limit register spilling) and run Gaussian elimination on them rather than on scalar entries. The full matrix stays in shared memory (tiled-SHMEM) or global memory (tiled-DRAM), and only
10
T. Chenaille et al.
the tiles involved in the current step are loaded into registers. At each step, the diagonal tile is factored using scalar partial pivoting, then surrounding tiles are updated with a Schur complement according to the RL or LL schedule. To keep data local, the pivot row is sought and physically swapped inside that tile, and a lower tile is accessed for logical pivoting only when the best in-tile candidate falls below a tunable threshold. Tile-local pivoting admits slightly more pivot growth than textbook partial pivoting, a stability cost is quantified in Section 4.2. When the system dimension is not a multiple of the tile size, a phantom mechanism pads the trailing tile to full size and uses compile-time skips (cf Section 3.2). The RHS and pivot vector can independently stay in registers or follow the matrix in shared or global memory, where all operands can benefit from coalesced accesses. 3.3
Embedding LUpp Solvers in Constitutive Law Evaluation
Embedding a single-threaded device-callable LUpp solver is straightforward. The MFront-generated per-integration-point code is wrapped in a single CUDA kernel hosting the whole constitutive law evaluation, within which the only changes are placing the LUpp operands in the residency targeted by the new solver (registers, shared, or global memory, instead of TFEL’s default register-resident layout) and substituting the new solver’s entry point for the TinyMatrixSolve call. The Newton loop and law evaluation are reused as-is. A multi-threaded device-callable LUpp solver additionally requires partitioning material variables, Jacobian assembly, and convergence checks across the T threads of each cooperating group, as well as adding synchronization barriers. Embedding a batched LUpp solver requires the MFront-generated code refactor outlined in Figure 1. The originally fused per-point computation is split into up to ten device kernels per global CPU-driven Newton iteration, interleaved with the batched library calls. On top of it, this scheme pays a host-device synchronization per iteration for the global convergence flag check, and routes intermediate GPU-resident data through global memory across each kernel boundary.
4
Performance Evaluation and Discussion
4.1
Experimental Setup
Hardware: All experiments are run on the IDRIS’s Jean Zay supercomputer. Computation use an NVIDIA H100 80GB SXM GPU with driver 595. The hostdriven batched configurations use one core of an Intel Xeon Platinum 8468 CPU. Software stack: Our prototype is compiled with CUDA 13.0.3 and host GCC 14.2.0, in C++20. External solvers come from MAGMA 2.10.0, Kokkos-Kernels 5.1.0, cuSolverDx 0.4.0 with in-house patches, SNLS 0.4.4, and cuBLAS from CUDA 13.0.3 SDK. The TFEL/MFront version is the master branch (May 2026). All computations are in double precision, with no fast-math optimizations. Compilation: Device code is built with –expt-relaxed-constexpr and -O3. Solvers requiring cross-translation-unit device optimization additionally enable relocatable device code (-rdc=true) and device link-time optimization (-dlto). Host code is built with -O3 and -DTFEL_NO_RUNTIME_CHECK_BOUNDS.
Solving of Many Tiny General Linear Systems on GPUs
11
Kernel launch configuration: Every in-house kernel selects its launch config through a uniform heuristic. To maximize occupancy and minimize tail effect, we retain the block size maximizing O + W , with O the theoretical occupancy ratio and W the average wave fill ratio. Batched solvers (cuBLAS, MAGMA) use their library-defined launch configuration unchanged. Measurement protocol: Each measurement is repeated 11 times (one warmup plus ten timed runs), and we report the median. We verify numerical correctness by comparison of the solution against a LAPACK reference solution. Detailed GPU profiling was performed using Nsight Compute. Reporting convention: Some curves are labelled “best X”. At each problem instance, this label reports the variant of solver family X with the smallest kernel duration. This reflects some realistic production cases, where a size-aware runtime dispatcher selects among specializations of the same kernel family. 4.2
Pure LUpp Performance
Pure-LUpp performance is measured on batches of 105 independent linear systems with two input distributions, neither of which produces singular matrices. The default draws matrix and RHS entries uniformly from [−0.5, +0.5], mimicking the standard LINPACK benchmark configuration. A stress distribution instead draws from [−5 × 10−10 , +5 × 10−10 ] to force small magnitudes and out-of-tile pivoting in all our tiled solvers, which all use a 10−10 pivot acceptance threshold. The default distribution produces no out-of-tile pivots. Under the stress one, depending on the tile size, roughly half the systems trigger at least one, averaging 1.5 per triggering system. We quantify the stability cost (Section 3.2) using the per-system normwise backward error (BE). Because pivoting is initially confined to a tile, BE grows as tiles shrink. A single tile recovers textbook LUpp and its BE order, while TileSize=2 gives the worst BE. Under the default distribution, the worst-case BE metrics across all matrix sizes are 10−15 (median), 10−14 (mean), and 10−9 (maximum). Under the stress distribution, out-of-tile pivoting widens pivot search and improves stability. Worst-case metrics then drop to 10−16 (median), 10−16 (mean), and 10−13 (maximum). Figure 2 shows single-threaded throughput across n ∈ [2, 32]. The in-house Tiled family dominates the [4, 32] range, even under the stress distribution. All library-provided single-threaded solvers (SNLS, cuSolverDx-thread, KokkosKernels, TFEL-baseline) underperform: they conceptually stream the matrix scalar by scalar from their target residency into registers. Pivot indirections map directly onto memory accesses, and since pivot sequences differ across systems, these accesses are uncoalesced. Instead, tiled solvers improve data locality by caching full tiles in registers, and their static indexings guards prevent local memory usage [37]. Without out-of-tile pivoting, their pivot indirections stay in registers, sparing them the uncoalesced accesses the library solvers incur. Across the n ∈ [2, 32] range, the best Tiled-SHMEM variant issues a median 4× fewer shared-memory loads and 18× less local-memory traffic than TFEL-baseline in shared memory; the best Tiled-DRAM variant moves a median 3× fewer DRAM bytes than TFEL-baseline in DRAM.
12
T. Chenaille et al.
Fig. 2: Throughput of the single-threaded LUpp solvers on 105 linear systems (pure LUpp solving). Higher is better. Under the default distribution, the residency of the winning tiled variant is selected by a trade-off between access cost and parallelism. Tiled-SHMEM keeps matrix accesses on-chip, but its shared-memory footprint lowers occupancy in steps: 4 warps/SM up to n = 14, 3 from n = 15, 2 from n = 18, and 1 from n = 22. Tiled-DRAM pays off-chip accesses but keeps 8 warps/SM throughout. Up to n = 21, Tiled-SHMEM’s cheaper accesses dominate; from n = 22, its limited parallelism erases this edge and both residencies perform within ∼ 10%, with no consistent winner. The RL/LL choice follows the same trade-off and inverts with the residency. On Tiled-DRAM, LL wins because it reads each prior tile only when needed instead of repeatedly updating trailing tiles. On TiledSHMEM, shared-memory accesses are cheap, so LL gains little by reducing them. Meanwhile, it issues ∼ 30% more instructions than RL, which is why RL wins from n ≥ 14. For n ≤ 13, the matrix has too few tiles for either schedule to gain a measurable edge (< 5%). The optimal Tile Size (TS) follows the dimension while the matrix fits in one or two tiles (n ≤ 11), then locks to TS ∈ {4, 5}, regardless of residency and distribution. Larger tiles amortize more updates per tile fetched, but register pressure caps the benefit (Section 3.2). Under the stress distribution, out-of-tile pivoting fires frequently and the winners shift. Tiled-SHMEM now dominates the whole range: on-chip access remains cheaper for the now-frequent out-of-tile pivot loads, and it keeps the lead in the sub-warp NTPB regime it enters from n = 31 (Section 3.2), which stays shallow here (≤ 4 idle lanes). A deeper sub-warp, at larger n or with smaller shared-memory budgets, would favor Tiled-DRAM. The RL/LL choice shifts too: LL now wins only at some small sizes (n ≤ 10); RL dominates from n ≥ 11. This is because LL’s out-of-tile pivot search must replay all prior L·U tile
Solving of Many Tiny General Linear Systems on GPUs
13
updates on each candidate row in the lower tiles, since LL has not materialized them into the trailing matrix. Since each candidate row replays every lower tile, this overhead grows with the tile count. RL avoids it because its trailing matrix is already updated, so candidate evaluation reduces to a raw load.
Fig. 3: Throughput of the multi-threaded LUpp solvers on 105 linear systems (pure LUpp solving). Higher is better. Figure 3 reports the multi-threaded throughput, with the best singlethreaded Tiled curve from the previous figure overlaid as a cross-paradigm anchor. The most striking observation is that this single-threaded curve outperforms every multi-threaded solver for n ∈ [2, 21], and remains competitive for the larger dims. It is overtaken by best MAGMA-custom and MAGMA only from n = 22. It always closely matches or beats cuSolverDx-block, and is always strictly faster than cuBLAS. In pure LUpp, over roughly the two thirds of the tiny regime, Tiled gains more from its data locality than multi-threaded solvers gain from their lower register pressure and higher occupancy. Among multi-threaded solvers, best MAGMA-custom consistently outperforms best cuSolverDx-block, and outperforms MAGMA for n ∈ [2, 24]. The gain comes from MAGMA-custom’s multiple rows per thread, and multiple systems per warp. Both eliminate the idle-lane waste of MAGMA’s one-rowper-thread, one-system-per-block design. At n = 25, systems are too large for our packing optimizations to stay effective, and the two solvers track each other closely. From n ≥ 26, the best MAGMA-custom variant uses RowsPerThread=SystemsPerWarp=1. It is therefore conceptually equivalent to the original MAGMA algorithm from this point on. The residual perfor-
14
T. Chenaille et al.
mance gap reflects minor implementation differences (generalization overhead and launch configuration), for an average 12% throughput penalty. Turning to the batched solvers cuBLAS and MAGMA: cuBLAS is consistently slower than MAGMA for n ∈ [2, 32], confirming the benefit of registerblocking used by MAGMA. The latter is itself beaten by cuSolverDx-block for n ∈ [2, 13], a range where MAGMA’s idle-lane penalty is particularly high, with at least 19 of its 32 lanes idle. In this pure-LUpp setting, both batched solvers run as a single GPU kernel: they avoid the overheads they would incur in an application context such as the one studied in Section 4.3. In summary, thanks to its tile-granular register caching, the in-house Tiled family is the fastest option for n ∈ [4, 21] under both default and stress distributions. From there, while Tiled remains competitive, MAGMA-custom takes the lead up to n = 24. Beyond, MAGMA narrowly takes over. The next section examines how these solvers perform in our applicative case. 4.3
Full Constitutive Law Evaluation Performance
Full constitutive law evaluation performance is measured on a representative time step from a 3D simulation of the uniaxial compression of a nuclear fuel pellet, driven by an experimentally measured loading (compressive strain rate). The latter contains 105 finite-element quadrature points, integrated under a Norton viscoplasticity model. Each of these integration points yields a Newton problem with n = 12, with one linear system of this dimension to solve at every Newton iteration, plus one multi-RHS solve after convergence for the consistent tangent operator. Newton terminates on the first of the following: residual norm below 10−12 (always reached on this input), or 100 iterations (never reached). Convergence is highly homogeneous across the 105 points: 2.1 iterations on average, 3 at most. On this real data, Tiled solvers never block nor slow down Newton convergence. Their slight stability cost (Section 4.2) thus carries no penalty in this applicative setting. Furthermore, they trigger no out-of-tile pivoting.
Fig. 4: Throughput of all LUpp solvers on 105 integration points (Norton constitutive law evaluation). Systems of dimension n = 12. Higher is better.
Solving of Many Tiny General Linear Systems on GPUs
15
Figure 4 shows the throughput of every solver on this workload. The Tiled solver dominates by a wide margin. It is 4.7× faster than the second-best solver (best MAGMA-custom). It reaches 5.9× over the best library-provided devicecallable competitor (TFEL-baseline), and 17.7× over the best batched library (MAGMA). Strikingly, however, it exposes metrics typically associated with poor performance, with 8.9% occupancy and 255 registers per thread. Multi-threaded device-callable solvers perform proportionally worse in this applicative setting than in pure LUpp: at n = 12 in pure LUpp, the speedup of best Tiled over best MAGMA-custom was only 1.4×, and over best cuSolverDxblock was 2.6×. Here it grows to 4.7× and 6.5× respectively. Two factors compound. First, the multi-threaded paradigm pays additional synchronization barriers across the constitutive law evaluation. Second, the registers consumed by the Newton state and material-law variables raise local memory usage and push these solvers off the high-occupancy plateau they enjoyed in pure LUpp: best cuSolverDx-block drops from 48.2% to 18.1% achieved occupancy, and best MAGMA-custom drops from 23.8% to 11.9%. The best Tiled solver is bottlenecked only by its shared-memory footprint, so its occupancy stays at 8.9%. The batched solvers cuBLAS and MAGMA sit behind every device-callable alternative. In this applicative setting, they pay the paradigm overhead discussed in Section 2.4. One component of this overhead, the dragging of alreadyconverged points, stays marginal here thanks to Newton’s homogeneous convergence (2.1 iterations on average, 3 at most). Both solvers would degrade further still on workloads with broader convergence dispersion.
5
Conclusion
Efficient LUpp solving is a key component for porting implicit constitutive law evaluation to GPUs. We benchmarked batched, multi-threaded device-callable, and single-threaded device-callable strategies for tiny LUpp, including two inhouse solver families. Batched solvers are slowed down by applicative overheads and remain behind the TFEL library baseline. Among multi-threaded devicecallable solvers, our MAGMA-custom solver improves on MAGMA library in pure LUpp for linear systems of dimension n ∈ [2, 24], but brings limited acceleration in the full constitutive workflow because of register pressure and synchronization costs induced by the multi-threaded approach. The strongest gains come from our single-threaded Tiled solver, which exploits tile-level register-caching. On the studied constitutive workload, this solver reaches 5.9× over the TFEL baseline, 6.5× over cuSolverDx, and 17.7× over MAGMA. Further work includes evaluating performance on newer GPUs with other backends and constitutive laws. Studying the Tiled solver data locality on CPU is another direction. Integrating it into TFEL/MFront (open-source) may follow, requiring only minor changes to MFront’s code-generation logic. Acknowledgments. This work was granted access to HPC computing and storage resources from GENCI at CNRS-IDRIS under grant 2026-101137, on the H100 partition of the Jean Zay supercomputer. Disclosure of Interests. The author has no competing interests to declare that are relevant to the content of this article.
16
T. Chenaille et al.
References 1. TOP500: TOP500 List - November 2025. https://top500.org/lists/top500/ list/2025/11/ 2. Abdelfattah, A., et al.: Batched one-sided factorizations of tiny matrices using GPUs: challenges and countermeasures. J. Comput. Sci. 26, 226–236 (2018). https: //doi.org/10.1016/j.jocs.2018.01.005 3. Balos, C.J., et al.: SUNDIALS time integrators for exascale applications with many independent systems of ordinary differential equations. Int. J. High Perform. Comput. Appl. 39(1), 123–146 (2025). https://doi.org/10.1177/10943420241280060 4. Anzt, H., et al.: Variable-size batched LU for small matrices and its integration into block-Jacobi preconditioning. In: Proceedings of the 46th ICPP, pp. 91–100 (2017). https://doi.org/10.1109/ICPP.2017.18 5. Helfer, T., et al.: Introducing the open-source mfront code generator: Application to mechanical behaviours and material knowledge management within the PLEIADES fuel element modelling platform. Comput. Math. Appl. 70(5), 994–1023 (2015). https://doi.org/10.1016/j.camwa.2015.06.027 6. Helfer, T., et al.: The MFrontGenericInterfaceSupport project. J. Open Source Softw. 5(48), 2003 (2020). https://doi.org/10.21105/joss.02003 7. Bernaud, S., et al.: PLEIADES: a numerical framework dedicated to the multiphysics and multiscale nuclear fuel behaviour simulation. Ann. Nucl. Energy 205, 110577 (2024). https://doi.org/10.1016/j.anucene.2024.110577 8. Helfer, T., et al.: MFEM/MGIS, a HPC mini-application targeting nonlinear thermo-mechanical simulations of nuclear fuels at mesoscale. J. Open Source Softw. 10(108), 7719 (2025). https://doi.org/10.21105/joss.07719 9. Kelley, C.T.: Solving Nonlinear Equations with Newton’s Method. SIAM, Philadelphia (2003). https://doi.org/10.1137/1.9780898718898 10. Higham, N.J.: Accuracy and Stability of Numerical Algorithms. 2nd edn. SIAM, Philadelphia (2002) 11. Golub, G.H., Van Loan, C.F.: Matrix Computations. 4th edn. Johns Hopkins University Press, Baltimore (2013) 12. Buttari, A., et al.: A class of parallel tiled linear algebra algorithms for multicore architectures. Parallel Comput. 35(1), 38–53 (2009). https://doi.org/10.1016/j. parco.2008.10.002 13. Haidar, A., et al.: A guide for achieving high performance with very small matrices on GPU: a case study of batched LU and Cholesky factorizations. IEEE Trans. Parallel Distrib. Syst. 29(5), 973–984 (2018). https://doi.org/10.1109/TPDS.2017. 2783929 14. Abdelfattah, A., et al.: Progressive optimization of batched LU factorization on GPUs. In: Proc. 2019 IEEE High Perform. Extreme Comput. Conf. (HPEC), pp. 1–6 (2019). https://doi.org/10.1109/HPEC.2019.8916270 15. Villa, O., et al.: Power/performance trade-offs of small batched LU based solvers on GPUs. In: Proc. Euro-Par 2013. LNCS, vol. 8097, pp. 813–825. Springer, Berlin, Heidelberg (2013). https://doi.org/10.1007/978-3-642-40047-6_81 16. Lei, X., et al.: High-performance batched LU decomposition on GPU. In: Chen, C.-H., et al. (eds.) Applied Mathematics, Modeling and Computer Simulation, pp. 380–396. IOS Press (2022). https://doi.org/10.3233/ATDE221053 17. Anderson, M.J., et al.: A predictive model for solving small linear algebra problems in GPU registers. In: Proc. 2012 IEEE 26th IPDPS, pp. 2–13 (2012). https://doi. org/10.1109/IPDPS.2012.11
Solving of Many Tiny General Linear Systems on GPUs
17
18. Haidar, A., et al.: Batched matrix computations on hardware accelerators based on GPUs. Int. J. High Perform. Comput. Appl. 29(2), 193–208 (2015). https: //doi.org/10.1177/1094342014567546 19. Abdelfattah, A., et al.: Factorization and inversion of a million matrices using GPUs: challenges and countermeasures. Procedia Comput. Sci. 108, 606–615 (2017). https://doi.org/10.1016/j.procs.2017.05.250 20. Abdelfattah, A., et al.: MAGMA: Enabling exascale performance with accelerated BLAS and LAPACK for diverse GPU architectures. Int. J. High Perform. Comput. Appl. 38(5), 468–490 (2024). https://doi.org/10.1177/10943420241261960 21. NVIDIA: cuBLAS Library Homepage. https://developer.nvidia.com/cublas 22. Filipovic, J., et al.: Optimizing CUDA code by kernel fusion: application on BLAS. J. Supercomput. 71(10), 3934–3957 (2015). https://doi.org/10.1007/ s11227-015-1483-z 23. NVIDIA: cuSolverDx Library. https://docs.nvidia.com/cuda/cusolverdx/ 24. Rajamanickam, S., et al.: Kokkos Kernels: performance portable sparse/dense linear algebra and graph kernels. arXiv preprint arXiv:2103.11991 (2021). https: //doi.org/10.48550/arXiv.2103.11991 25. Wayne, B.M., et al.: SNLS: A C++ library for solving small non-linear systems. (2018). https://doi.org/10.11578/dc.20181217.9 26. Simo, J.C., Taylor, R.L.: Consistent tangent operators for rate-independent elastoplasticity. Comput. Methods Appl. Mech. Eng. 48(1), 101–118 (1985). https: //doi.org/10.1016/0045-7825(85)90070-2 27. de Souza Neto, E.A., et al.: Computational Methods for Plasticity: Theory and Applications. Wiley, Chichester (2008) 28. Savage, D.J., Knezevic, M.: Computer implementations of iterative and noniterative crystal plasticity solvers on high performance graphics hardware. Comput. Mech. 56(4), 677–690 (2015). https://doi.org/10.1007/s00466-015-1194-6 29. Hu, T., Messner, M.C.: NEML2: An efficient and modular multiphysics constitutive modeling library for hybrid computing environments. SoftwareX 31, 102302 (2025). https://doi.org/10.1016/j.softx.2025.102302 30. Hu, F., et al.: Efficient GPU-computing simulation platform JAX-CPFEM for differentiable crystal plasticity finite element method. npj Comput. Mater. 11, 49 (2025). https://doi.org/10.1038/s41524-025-01528-2 31. Latyshev, A., et al.: Expressing general constitutive models in FEniCSx using external operators and algorithmic automatic differentiation. J. Theor. Comput. Appl. Mech. (2025). https://doi.org/10.46298/jtcam.14449 32. Harbour, L., et al.: 4.0 MOOSE: Enabling massively parallel Multiphysics simulation. SoftwareX (2024). https://doi.org/10.1016/j.softx.2024.101754 33. Barton, N.R., et al.: LLNL/ExaCMech. (2018). https://doi.org/10.11578/dc. 20190809.2 34. Blühdorn, J., et al.: AutoMat: automatic differentiation for generalized standard materials on GPUs. Comput. Mech. 69(2), 589–613 (2022). https://doi.org/10. 1007/s00466-021-02105-2 35. Morgan, N., et al.: Enabling parallel performance and portability of solid mechanics simulations across CPU and GPU architectures. Information 15(11), 716 (2024). https://doi.org/10.3390/info15110716 36. Idaho National Lab.: GPU capabilities in the MOOSE framework through Kokkos. GitHub issue #30655 (2025). https://github.com/idaholab/moose/issues/30655 37. Bourgeois, R.: Five basic performance advice for porting kernels to the GPU. Github personal blog. https://rbourgeois33.github.io/posts/post1/ #precautions-when-using-static-arrays-for-temporary-storage