solver_static.hpp File Reference

TiledLUpp solver: LU with partial pivoting on tile grids, one thread / work-item per system. More...

#include <tdls/core/math.hpp>
#include <tdls/core/target.hpp>
#include <type_traits>
#include <tdls/solvers/tiled_lupp/config.hpp>
#include <tdls/solvers/tiled_lupp/tile_operations.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::TiledLUppSolverStatic< T, N, Config >
	Tiled dense LU factorization with logical partial pivoting and out-of-tile pivot recovery, solving one NxN system per call. More...

Macros
#define	TDLS_LUPP_A(r, c)
	Element (r, c) of S, the stored matrix read row-major (A or A^T, see the file documentation): contiguous under internal residency, strided otherwise.
#define	TDLS_LUPP_PIV(i)
	Pivot entry i: contiguous under internal residency, strided otherwise.
#define	TDLS_LUPP_X(i)
	Entry i of the solution vector: contiguous under internal residency, strided otherwise.
#define	TDLS_LUPP_B(i)
	Entry i of the right-hand side: contiguous under internal residency, strided otherwise.
#define	TDLS_LUPP_XW(w, i)
	Entry i of column w of a multi right-hand-side block: pass_width contiguous columns under internal residency, xcol_stride-strided columns in remote memory. pass_width = 1 collapses to TDLS_LUPP_X exactly.
#define	TDLS_LUPP_BW(w, i)
	Entry i of column w of a multi right-hand-side block of b, addressed exactly as TDLS_LUPP_XW (b and x share both strides).
#define	TDLS_LUPP_Y(i)
	Entry i of the fused right-hand side of solve_inplace (under the row-major layout y follows the matrix rows through pivoting).
#define	TDLS_LUPP_FOR(i, first, last, ...)
	Loop of the physical row interchanges, unrolled under Config.unroll_loops: the bounds are compile-time constants, the triangle selected inside the body, so that every index is known at compile time.

Detailed Description

TiledLUpp solver: LU with partial pivoting on tile grids, one thread / work-item per system.

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 matrix is split into a grid of tile_size x tile_size tiles; Gaussian elimination runs on tiles instead of scalars. The tiles of the current step are copied into local arrays, the working set kept in registers. The full matrix can stay in remote memory (shared, global or plain host memory), walked with a stride, or be a caller-local array (the residency booleans). Pivoting is logical by default (Config.row_interchange): piv[] maps logical row -> physical row, rows are physically swapped only inside the diagonal tile, and each tile goes back to the physical rows that piv gives, an index known at run time only. A caller-local matrix and its pivots are thus indexed at run time: on GPU they live in the local memory of the thread, which the caches serve, and only the working tiles stay in registers. Inside the tile, the row in place keeps the pivot when it reaches both Config.relative_pivot_threshold times the best in-tile candidate and Config.oot_pivot_threshold (threshold pivoting, in-tile only). When the best in-tile pivot falls below Config.oot_pivot_threshold, the search extends below the tile (out-of-tile pivoting) and candidate values are corrected on the fly for the eliminations they have not received yet.

Under Config.row_interchange = RowInterchange::Physical, the whole matrix is one tile and the rows are interchanged as in LAPACK, through a test of every unrolled row index against the pivot row: every index is known at compile time, so that a caller-local matrix, its pivots and its right-hand sides stay in registers on GPU, for matrices small enough to fit there. The pivots and the factorization are those of the logical scheme over a single tile; the pivot array holds the interchanges, piv[j] the row exchanged with row j at step j.

Layout. The solver factors S, the stored matrix read row-major: element (r,c) of S lives at A[((r)*N+(c))*A_stride]. S = A under the row-major layout, the default, and S = A^T under the column-major one. The tiles are thus read in address order whatever the layout. Under the row-major layout, P A = L U and the solve is L y = P b, folded into the factorization by solve_inplace, then U x = y. Under the column-major layout, P A^T = L U, so A = U^T L^T P. The pivoting interchanges columns of A: this is partial pivoting of A^T, as stable as partial pivoting of A. The solve is the mirror of the row-major one. First U^T z = b, folded into the factorization by solve_inplace. Then L^T w = z, which divides by nothing since L has a unit diagonal, and x = P^T w. Both layouts run the same factorization on the same bytes. A stored column-major gives, bitwise, the factored array and the pivots of the row-major solver given A^T. The layout is a compile-time choice. The rest of this documentation speaks of the rows of S.

Addressing convention (backend-agnostic, no accessor objects). The solver never sees a thread index: callers pre-offset every remote pointer with the lane/system index and pass a runtime stride, or a stride known at compile time as a template argument (see tdls::dynamic_stride). The pivot array and the right-hand sides follow the same convention. The flat element index (for the matrix, N*N) must stay below 2^31; element offsets, index times stride, are computed in 64-bit arithmetic. Each of the three arrays has an internal residency mode (plain caller-local array, stride ignored) selected by the internal_rhs / internal_piv / internal_matrix template booleans. The ternary on a constexpr bool costs nothing. The single form base[e*stride] covers every batched layout without touching the solver: AoS (stride 1), SoA (stride = batch size), AoSoA (stride = W).

Tile grid. full_tiles = N / tile_size full tiles per dimension, plus a last tile of extent last_tile_tail = N - full_tiles * tile_size when N is not a multiple of tile_size; tile_size may exceed N, the grid then being that single partial tile (full_tiles = 0, last_tile_tail = N). Every sweep is a runtime loop over the full tiles followed by an if constexpr (last_tile_tail > 0) epilogue instantiated with the tail extent; phantom slots of partial tiles are skipped at compile time by the extent template parameters.

Factored format. The diagonal of the factored matrix holds the RECIPROCALS of the U pivots (consumers multiply; each division is paid once per pivot). This deliberately differs from the LAPACK getrf convention. A factorization produced here must be consumed by the substitution routines of this solver.

Originally developed and validated in the tfelGPU project.

Macro Definition Documentation

TDLS_LUPP_A

#define TDLS_LUPP_A	(		r,
			c )

Value:

    A[internal_matrix ? unsigned((r) * N + (c)) : TDLS_STRIDED((r) * N + (c), A_stride)]

Element (r, c) of S, the stored matrix read row-major (A or A^T, see the file documentation): contiguous under internal residency, strided otherwise.

TDLS_LUPP_B

#define TDLS_LUPP_B

(

i

)

Value:

b[internal_rhs ? unsigned(i) : TDLS_STRIDED(i, rhs_stride)]

Entry i of the right-hand side: contiguous under internal residency, strided otherwise.

TDLS_LUPP_BW

#define TDLS_LUPP_BW	(		w,
			i )

Value:

    b[internal_rhs ? unsigned((w) * N + (i))                                                       \
                   : TDLS_STRIDED(i, rhs_stride) + TDLS_STRIDED_WIDE(w, xcol_stride)]

Entry i of column w of a multi right-hand-side block of b, addressed exactly as TDLS_LUPP_XW (b and x share both strides).

TDLS_LUPP_FOR

#define TDLS_LUPP_FOR	(		i,
			first,
			last,
			... )

Value:

    if constexpr (Config.unroll_loops) {                                                           \
        TDLS_UNROLL_FORCE                                                                          \
        for (int i = (first); i < (last); ++i) {                                                   \
            __VA_ARGS__                                                                            \
        }                                                                                          \
    } else {                                                                                       \
        for (int i = (first); i < (last); ++i) {                                                   \
            __VA_ARGS__                                                                            \
        }                                                                                          \
    }

Loop of the physical row interchanges, unrolled under Config.unroll_loops: the bounds are compile-time constants, the triangle selected inside the body, so that every index is known at compile time.

TDLS_LUPP_PIV

#define TDLS_LUPP_PIV

(

i

)

Value:

piv[internal_piv ? unsigned(i) : TDLS_STRIDED(i, piv_stride)]

Pivot entry i: contiguous under internal residency, strided otherwise.

TDLS_LUPP_X

#define TDLS_LUPP_X

(

i

)

Value:

x[internal_rhs ? unsigned(i) : TDLS_STRIDED(i, rhs_stride)]

Entry i of the solution vector: contiguous under internal residency, strided otherwise.

TDLS_LUPP_XW

#define TDLS_LUPP_XW	(		w,
			i )

Value:

    x[internal_rhs ? unsigned((w) * N + (i))                                                       \
                   : TDLS_STRIDED(i, rhs_stride) + TDLS_STRIDED_WIDE(w, xcol_stride)]

Entry i of column w of a multi right-hand-side block: pass_width contiguous columns under internal residency, xcol_stride-strided columns in remote memory. pass_width = 1 collapses to TDLS_LUPP_X exactly.

TDLS_LUPP_Y

#define TDLS_LUPP_Y

(

i

)

Value:

y[internal_rhs ? unsigned(i) : TDLS_STRIDED(i, rhs_stride)]

Entry i of the fused right-hand side of solve_inplace (under the row-major layout y follows the matrix rows through pivoting).