solver_static.hpp File Reference

CooperativeLUpp solver: LU with partial pivoting held in the registers of a group of threads, one system per group. More...

#include <type_traits>
#include <utility>
#include <tdls/core/group.hpp>
#include <tdls/core/macros.hpp>
#include <tdls/core/math.hpp>
#include <tdls/core/target.hpp>
#include <tdls/solvers/cooperative_lupp/config.hpp>

Include dependency graph for solver_static.hpp:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	tdls::CooperativeLUppSolverStatic< T, N, Config >
	Dense LU factorization with partial pivoting held in the registers of a group of threads, solving one NxN system per group. More...

Macros
#define	TDLS_COOP_LUPP_A(K, c)
	Element (row of slot K, c) of S, the stored matrix read column-major (A under the column-major layout, A^T under the row-major one): slice of the calling thread under internal residency, whole strided matrix otherwise. Consecutive threads hold consecutive rows of S, thus consecutive addresses.
#define	TDLS_COOP_LUPP_PIV(K)
	Pivot entry of the row of slot K: slice of the calling thread under internal residency, whole strided array otherwise.
#define	TDLS_COOP_LUPP_B(K)
	Right-hand-side entry of the row of slot K: slice of the calling thread under internal residency, whole strided vector otherwise.
#define	TDLS_COOP_LUPP_X(K)
	Solution entry of the row of slot K, addressed as TDLS_COOP_LUPP_B.
#define	TDLS_COOP_LUPP_Y(K)
	Entry of the row of slot K of the right-hand side solved in place, addressed as TDLS_COOP_LUPP_B.

Functions
template<typename Fn , int... Ks>
TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU constexpr void	tdls::detail::unroll_K_impl (Fn &fn, std::integer_sequence< int, Ks... >) noexcept
	Calls fn once per index of the sequence, in increasing order, with the index as a std::integral_constant.
template<int R, typename Fn >
TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU constexpr void	tdls::detail::unroll_K (Fn &&fn) noexcept
	Compile-time iteration over the row slots K = 0 .. R-1 of a thread.

Detailed Description

CooperativeLUpp solver: LU with partial pivoting held in the registers of a group of threads, one system per group.

Author

Tristan Chenaille

Copyright

Copyright (C) 2026 CEA. All rights reserved. This project is publicly released under the BSD 3-Clause License (see the LICENSE file). CEA may also distribute it under specific licensing conditions.

The elimination and substitution kernels of this file are derived from MAGMA 2.10.0 (Univ. of Tennessee, Knoxville; Univ. of California, Berkeley; Univ. of Colorado, Denver), file magmablas/zgesv_batched_small.cu, function zgesv_batched_small_device, authors Azzam Haidar and Ahmad Abdelfattah. Those portions remain subject to the MAGMA license (BSD 3-Clause), reproduced in full in the comment at the top of this file as required by it.

One system is solved by a group of threads_per_system threads. The solver factors S, the stored matrix read column-major: S = A under the column-major layout, S = A^T under the row-major one, the default. Thread tx holds rows tx, tx + threads_per_system, tx + 2 * threads_per_system and so on of S, rows_per_thread of them, in registers: rows of A under the column-major layout, columns of A under the row-major one. In both layouts the threads of a group thus read consecutive addresses, the access pattern the design of MAGMA rests on. The threads cooperate through a small workspace (shared memory on GPU) and a barrier, deduced or provided by the caller. The factorization is the right-looking LU with partial pivoting of MAGMA's register-blocked small-system kernel, applied to S: same operations, same order and, with Config.relative_pivot_threshold = 1, same pivot choice; a smaller threshold keeps the row in place more often, and the pivots then differ. Under the column-major layout, P A = L U and the solve is MAGMA's: L y = P b folded into the factorization, then U x = y. Under the row-major layout, P A^T = L U, so A = U^T L^T P: the pivoting interchanges columns of A, which is partial pivoting of A^T and as stable, and the solve is the mirror of MAGMA's: U^T z = b folded into the factorization, then L^T w = z, which divides by nothing since L has a unit diagonal, and x = P^T w. The rest of this documentation speaks of the rows of S. By default pivoting is logical: rows never move, each thread tracks the position of its rows in the pivoted order (MAGMA's rowid) and the pivot array records it. Config.row_interchange selects physical row interchanges instead, the scheme of LAPACK: the rows move between the threads, so that the row at position k always lives in slot k / threads_per_system of thread k % threads_per_system. Both schemes choose the same pivots and run the same operations in the same order. On a matrix they do not declare singular, their solutions are bitwise identical, up to the multiply-adds a compiler may fuse differently in each (see CooperativeLUppConfig::row_interchange). The mapping of rows to threads changes nothing in the arithmetic either, so every value of rows_per_thread produces bitwise-identical results on identical inputs. Config.relative_pivot_threshold below 1 relaxes the pivot choice, the same way in both schemes. The runtime dimension variant is CooperativeLUppSolverDynamic (solver_dynamic.hpp).

Calling convention. Every thread of the group calls the same entry point, with the same template arguments, the same operands (or, for an internal operand, its own slice), the same workspace and its own rank tx in [0, threads_per_system). The solver never sees a global thread index: callers pre-offset every remote pointer with the system index and pass a runtime stride, as for the TiledLUpp solvers, or a stride known at compile time as a template argument (see tdls::dynamic_stride).

Barrier. The last argument, sync, is the barrier of the group, and most calls omit it. By default, tdls::AutoSync, the solver deduces it from the compilation target: nothing with one thread per system; on GPU, the lanes of the group within their warp or wavefront, checked on entry; on CPU, a barrier held in the workspace. A caller passes its own barrier only where the compiler asks for one, or for a group that the deduced barrier does not know, such as a group larger than a warp. Any callable taking no argument serves, invoked by every thread of the group, that makes the memory writes of each thread visible to all of them, in the workspace and in the operands: a block, SYCL group, Kokkos team or CPU thread barrier. It is then used as is, without any check. make_sync returns the deduced barrier, for a caller that also needs it for its own exchanges. tdls::NoSync, no barrier at all, is refused at compile time for a group of several threads. core/group.hpp details these cases, and what the solver checks.

Fixed barrier sequence. Every entry point executes a sequence of barriers set by N, the row interchanges, the layout and the entry point alone: it never depends on the data, not even on a singular verdict, which is returned at the end instead of by an early exit. A barrier wider than the group (a whole warp, a work-group, a team) is therefore valid as well, provided every thread of its scope makes the same sequence of calls. Between two barriers, no workspace element is written by one thread while another thread accesses it, so the solver is free of data races even when the threads of a group do not run in lockstep.

Inputs on entry. Each thread reads only the entries of its own rows of S and of the right-hand side: the matrix entries S(tx + K * threads_per_system, c) and the right-hand-side entries tx + K * threads_per_system. A caller that builds them in that thread needs no barrier before the call: under the row-major layout, the columns of A of the thread.

Results on return. Every entry point ends with a barrier. When a thread returns, every thread of the group has written its results (factored rows, pivot entries, solution entries), and they are visible to the whole group, whatever the residency of the operands: a thread may read the solution entries of the others at once. The workspace is free again on return as well, so the caller may reuse it immediately, for instance to exchange data between the threads.

Residency. Each operand has an internal mode, selected by the internal_rhs / internal_piv / internal_matrix template booleans: the operand is then the slice of the calling thread, a plain caller-local array, and the stride is ignored. In external mode the operand is the whole object, in memory reachable by every thread of the group, walked with the stride.

• matrix: external element (r, c) of A lives at A[((r)*N+(c))*A_stride] under the default row-major layout, A[((c)*N+(r))*A_stride] under the column-major one. Element (v, c) of S thus lives at A[((c)*N+(v))*A_stride] in both. The internal slice holds the rows_per_thread rows of S of the thread: element c of row tx + K * threads_per_system lives at A[c*rows_per_thread+K]. It is the block of A made of these vectors, stored in the layout of Config.layout: the N x rows_per_thread block of the columns of the thread, row-major, or the rows_per_thread x N block of its rows, column-major. With one thread per system, the slice is the whole matrix in that layout.
• pivot: under logical row interchanges, external entry r holds the position of physical row r of S in the pivoted order. This maps physical to logical rows, the inverse of the TiledLUpp convention. Under physical ones, entry k holds the original index of the row of S at position k, as in the TiledLUpp solvers. Under the row-major layout, the rows of S are the columns of A. The internal slice holds the entries of the rows of the thread, entry K for row tx + K * threads_per_system.
• right-hand side and solution: external entry r is row r; the internal slice holds entry tx + K * threads_per_system at index K.

Element offsets are computed in 64-bit arithmetic, as in the TiledLUpp solvers. The single form base[e*stride] covers every batched layout: AoS (stride 1), SoA (stride = batch size), AoSoA (stride = W).

Factored format. S is factored in place, in the stored array. Under logical row interchanges, the factorization stays in the physical rows of S: physical row r holds, before the column of its pivot step, the multipliers of L, and from that column on, its row of U. Under physical ones, the rows end in the pivoted order, as in LAPACK: row k holds the multipliers of L before column k, and its row of U from column k on. Both formats hold the same rows, bit for bit, at different places. Under the row-major layout, these rows of S are the columns of the stored array. The diagonal holds the pivots themselves, not their reciprocals. A factorization produced here must be consumed by the substitution routines of this solver, under the same configuration.

Workspace. Under logical row interchanges, 3 * N elements per group, MAGMA's sB, sx and dsx, in this order. Under physical ones, 2 * N + 2 * threads_per_system + 5 elements: the candidates of the threads and their positions, the magnitude of the row in place, then two exchange records of N + 2 elements, the rows that leave and enter the position of the step. The substitutions use the first 2 * N + 2 elements at most. Two more elements come first, those of the deduced CPU barrier: on CPU, they must be zero before the first call of a group of several threads, and the barrier then leaves the counter at zero and the flag at 0 or 1 between calls. They keep their place whatever the dimension and the configuration, so that a workspace may serve other solves. workspace_size counts them all. The workspace pointer carries no restrict qualifier: the other threads of the group write the workspace between two barriers, and a restrict pointer would let the compiler reuse a value read before the barrier. Every other operand keeps the qualifier, since each thread only accesses its own rows of it.

Departures from the MAGMA kernel. The arithmetic is MAGMA's, but for the relative threshold of the pivot choice; the other changes are structural.

• rows_per_thread rows per thread instead of one, so that several systems share a warp, phantom rows skipped (at compile time in a slot where no thread holds a real row, at run time otherwise).
• Device-callable entry points with a deduced or caller-provided barrier, instead of a kernel with one system per block.
• The back substitution keeps the solution entry x_i in a register of the thread owning entry i. MAGMA writes it to sB[i] in the very barrier interval where every thread reads sB[i]: correct in lockstep execution only, a data race otherwise.
• Under logical row interchanges, the pivot output is the logical position of each row (rowid), not the LAPACK swap sequence, and the matrix keeps its physical row order: a separate substitution can then rebuild its state from the factored rows alone. The swap sequence (sipiv) is not computed.
• Under the row-major layout, MAGMA's factorization of A^T and the mirror solve U^T z = b, L^T w = z, x = P^T w: the threads hold columns of A, so that a group still reads consecutive addresses. Each substitution step takes one barrier instead of two, the exchange being double-buffered.
• Physical row interchanges as an option, with the same pivots and the same operations: one candidate per thread instead of one magnitude per row, two barriers per column instead of three, and no row moved while the pivot is in place. Under them a zero pivot row is not scaled to zero: the factorization of a singular matrix is unspecified, and a pivot reciprocal of 1 keeps it finite.
• A relative threshold of the pivot choice, threshold pivoting: 1 keeps the choice of MAGMA and LAPACK, a smaller value keeps the row in place more often, and the pivots and the results then differ.
• A pivot below Config.singular_floor is singular, instead of an exactly zero one. The default floor only adds the subnormal pivots, the criterion of every TDLS solver.
• Real scalar types: the magnitude of the pivot search is the absolute value, as in the real-precision variants of MAGMA.

Macro Definition Documentation

TDLS_COOP_LUPP_A

#define TDLS_COOP_LUPP_A	(		K,
			c )

Value:

    A[internal_matrix ? unsigned((c) * rows_per_thread + (K))                                      \
                      : TDLS_STRIDED((c) * N + tx + (K) * threads_per_system, A_stride)]

Element (row of slot K, c) of S, the stored matrix read column-major (A under the column-major layout, A^T under the row-major one): slice of the calling thread under internal residency, whole strided matrix otherwise. Consecutive threads hold consecutive rows of S, thus consecutive addresses.

TDLS_COOP_LUPP_B

#define TDLS_COOP_LUPP_B

(

K

)

Value:

    b[internal_rhs ? unsigned(K) : TDLS_STRIDED(tx + (K) * threads_per_system, rhs_stride)]

Right-hand-side entry of the row of slot K: slice of the calling thread under internal residency, whole strided vector otherwise.

TDLS_COOP_LUPP_PIV

#define TDLS_COOP_LUPP_PIV

(

K

)

Value:

    piv[internal_piv ? unsigned(K) : TDLS_STRIDED(tx + (K) * threads_per_system, piv_stride)]

Pivot entry of the row of slot K: slice of the calling thread under internal residency, whole strided array otherwise.

TDLS_COOP_LUPP_X

#define TDLS_COOP_LUPP_X

(

K

)

Value:

    x[internal_rhs ? unsigned(K) : TDLS_STRIDED(tx + (K) * threads_per_system, rhs_stride)]

Solution entry of the row of slot K, addressed as TDLS_COOP_LUPP_B.

TDLS_COOP_LUPP_Y

#define TDLS_COOP_LUPP_Y

(

K

)

Value:

    y[internal_rhs ? unsigned(K) : TDLS_STRIDED(tx + (K) * threads_per_system, rhs_stride)]

Entry of the row of slot K of the right-hand side solved in place, addressed as TDLS_COOP_LUPP_B.

Function Documentation

unroll_K()

template<int R, typename Fn >

TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU constexpr void tdls::detail::unroll_K ( Fn && fn )

constexprnoexcept

Compile-time iteration over the row slots K = 0 .. R-1 of a thread.

Inside fn, K is a constant expression, unlike the index of a runtime loop: the slot classification goes through if constexpr, and the register rows rA[K] are indexed at compile time even where no unroll pragma applies. The solver passes it only lambdas that never call the barrier, the one operation that may throw, hence noexcept.

Template Parameters

R	number of slots
Fn	callable type

Parameters

[in]

fn

callable, invoked as fn(std::integral_constant<int, K>{})

unroll_K_impl()

template<typename Fn , int... Ks>

TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU constexpr void tdls::detail::unroll_K_impl ( Fn & fn, std::integer_sequence< int, Ks... >  )

constexprnoexcept

Calls fn once per index of the sequence, in increasing order, with the index as a std::integral_constant.

Template Parameters

Fn	callable type
Ks	indices

Parameters

[in]

fn

callable, invoked as fn(std::integral_constant<int, K>{})