solver_static.hpp File Reference

Compile-time-shape solver of arrowhead (bordered block-diagonal) systems by static condensation, built on TiledLUpp. More...

#include <cmath>
#include <type_traits>
#include <utility>
#include <tdls/core/macros.hpp>
#include <tdls/core/math.hpp>
#include <tdls/core/target.hpp>
#include <tdls/solvers/arrowhead_lupp/config.hpp>
#include <tdls/solvers/arrowhead_lupp/storage_map.hpp>
#include <tdls/solvers/tiled_lupp/solver_static.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::detail::ArrowheadStreamedBlock< T, S, B, col, diag >
	Block k of a streamed solve of ArrowheadLUppSolverStatic (ArrowheadLUppSolverStatic::StreamedBlock): the working arrays of the solver themselves, which the caller fills through the accessors when the solver asks for the block; the factorization and the solve then run in place. Only a save is copied (D_k, and [A_kb | y_k] when the column-major check reads the forward stage), to store a refused block as given. Entries in logical order whatever the layout and the storage. The arrays come uninitialized, or holding another block: the caller writes every entry, zeros included, by assignment. More...
struct	tdls::detail::ArrowheadRecomputedBlock< T, S, B, col, diag >
	Block k of a recomputed solve (solve_inplace_recomputed, solve_inplace_recomputed_or_tangent): a StreamedBlock, the same accessors, and whether the solver asks for the block again, in its recovery. A provider taking a StreamedBlock& serves too; one taking a RecomputedBlock&, or auto&, reads the flag. The streamed entry points keep the StreamedBlock alone, so that their code is that of a solver without the recomputed ones. More...
struct	tdls::ArrowheadLUppSolverStatic< T, Shape, Config >
	Arrowhead solver with a compile-time shape. More...

Macros
#define	TDLS_ARROW_LUPP_A(idx)
	Stored entry at index idx of the matrix: contiguous under internal residency, strided otherwise.
#define	TDLS_ARROW_LUPP_Y(i)
	Entry i of the right-hand side and solution: contiguous under internal residency, strided otherwise.
#define	TDLS_ARROW_LUPP_YW(m, i)
	Entry i of column m of the right-hand sides and solutions: columns of dimension contiguous entries under internal residency, xcol_stride-strided columns otherwise. The function declares ycols, its number of columns: with one, m is ignored and the macro is TDLS_ARROW_LUPP_Y exactly.
#define	TDLS_ARROW_LUPP_W(e)
	Entry e of the workspace: contiguous under internal residency, strided otherwise.
#define	TDLS_ARROW_LUPP_F(r, c)
	Entry (r, c) of the front of the second pass, row-major at the start of the workspace, fcols columns per row.
#define	TDLS_ARROW_LUPP_SIDED(...)
	Under the side dispatch of nvc++ (TDLS_TARGET_GPU), when the host and the device choose differently (sided), returns the call, the variadic argument, run by the solver of the side of the code: SolverGpu on the device, SolverCpu on the host; nothing elsewhere and in constant expressions.
#define	TDLS_ARROW_LUPP_FOR(i, first, last, ...)
	for (int i = first; i < last; ++i) { body } under the unroll policy of the configuration: the two-branch if constexpr of the other families, the forced-unroll pragma lexically present in the first branch only. The body is the variadic argument.

Enumerations
enum class	tdls::detail::ArrowheadStreamedOutcome { solved , tangent , singular }
	Outcome of ArrowheadLUppSolverStatic::solve_inplace_streamed_or_tangent. More...
enum class	tdls::detail::ArrowheadRecomputedOutcome { solved , tangent , singular , deferred }
	Outcome of ArrowheadLUppSolverStatic::solve_inplace_recomputed and solve_inplace_recomputed_or_tangent. More...

Detailed Description

Compile-time-shape solver of arrowhead (bordered block-diagonal) systems by static condensation, built on TiledLUpp.

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 system has B border unknowns first, then NB blocks of S unknowns:

[ Abb  Ab0  Ab1 ... ] [xb]   [yb]
[ A0b  D0           ] [x0]   [y0]
[ A1b       D1      ] [x1] = [y1]
[ ...           ... ] [..]   [..]

Two blocks never couple. The solver reads only the entries of this pattern: the zeros are never read, so they need not be stored (see storage_map.hpp).

Algorithm. First pass, per block k:

1. D_k is copied into a local array and factored by TiledLUpp, a TiledLUppSolverStatic of dimension S.
2. The pivots are checked on the factors: no multiplier of the elimination of the block exceeds 1 / condensation_pivot_threshold, for its own rows (L) and for the border rows that share its columns (M = Cb U^-1). A block that fails, singular included, is refused and left untouched for the second pass.
3. Otherwise the block is condensed: its contribution goes to the Schur complement of the border, its recovery data in place of A_kb and y_k (see below; none for the canonical solves). Second pass, only when blocks were refused, in the workspace: each refused block is eliminated together with the border rows, by partial pivoting over its S rows and the B border rows. A border row may serve as pivot; a row of the block then takes its place in the border. A column whose best candidate falls below BlockSolver::singular_floor, zero or subnormal, makes the matrix singular. The border keeps B rows whatever the number of refused blocks. Its rows carry their couplings to the refused blocks still to come through a mixing matrix T (B x S under the diagonal coupling, B x B otherwise). Each refused block leaves x_k = c_k - G_k xb - H_k w, where w gathers the couplings A_bj x_j of the refused blocks after it. Then the border system (B x B) is solved by TiledLUpp, a TiledLUppSolverStatic of dimension B, and the blocks are recovered, the last one first.

Condensation of an accepted block, per layout:

• column-major, in one pass: the forward stage U^-T [A_kb | y_k] of the solve (TiledLUpp's forward_inplace_multirhs) holds the border multipliers of the check; the backward stage completes W_k = D_k^-1 [A_kb | y_k], stored in place of A_kb and y_k; the Schur complement receives -A_bk W_k; x_k = W_k(:, B) - W_k(:, 0:B) xb;
• row-major, diagonal coupling, or dense coupling under ArrowheadCondensation::MinTraffic (and column-major in several passes): the same with the whole solve, the check computing M = Cb U^-1 apart;
• row-major, dense coupling under ArrowheadCondensation::MinFlops, the default: V_k = L^-1 P [A_kb | y_k], the forward stage alone, is stored in place of A_kb and y_k, U in the upper triangle of D_k; the Schur complement receives -M V_k (M from the check: A_bk D_k^-1 = M L^-1 P); x_k = U^-1 (V_k(:, B) - V_k(:, 0:B) xb).

A refused block costs the second pass, never a failure: the solver succeeds whenever the full matrix is invertible, rounding aside. Like TiledLUpp, the solver does not detect non-finite entries: a block with one fails the pivot check and goes to the second pass, where they spread to the solution. A true return is no proof of a finite solution.

On exit y holds the solution. A holds the recovery data in place of A_kb, and in place of D_k for the refused blocks and, under the row-major layout with the dense coupling and MinFlops, for the accepted ones too; with schur_in_place, the coupled border rows of A hold the Schur complement of the accepted blocks. The other entries are only read. The canonical solves recover nothing and store no recovery data: they write in A only the Schur complement under schur_in_place and, streamed, the refused blocks. Without schur_in_place, a stored A thus comes out unchanged and serves the next columns.

The caller pre-offsets every remote pointer with the system index and passes a runtime stride, as for the other families, or a stride known at compile time as a template argument (see tdls::dynamic_stride); the internal_rhs / internal_work / internal_matrix template booleans declare a caller-local array instead (stride ignored). The workspace holds the front of the second pass: it is touched only when blocks are refused.

Streamed blocks. The *_streamed entry points take the blocks from a provider, called exactly once per block and in order when the solver condenses it, instead of reading them from A: the caller stores only the border, and the provider writes every entry of the block, zeros included, in the working arrays of the solver (StreamedBlock). A refused block is stored as given for the second pass. The arithmetic is that of the stored path, bit for bit when the compiler contracts multiply-adds alike in the two paths (-ffp-contract=off, or –fmad=false for nvcc, guarantees it; GCC with -march=native does not). The streamed path is selected at compile time: the stored entry points keep their signatures and their code. The entry points are noexcept: the provider must not throw.

Macro Definition Documentation

TDLS_ARROW_LUPP_A

#define TDLS_ARROW_LUPP_A

(

idx

)

Value:

A[internal_matrix ? unsigned(idx) : TDLS_STRIDED(idx, A_stride)]

Stored entry at index idx of the matrix: contiguous under internal residency, strided otherwise.

TDLS_ARROW_LUPP_F

#define TDLS_ARROW_LUPP_F	(		r,
			c )

Value:

TDLS_ARROW_LUPP_W((r) * fcols + (c))

Entry (r, c) of the front of the second pass, row-major at the start of the workspace, fcols columns per row.

TDLS_ARROW_LUPP_FOR

#define TDLS_ARROW_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__                                                                            \
        }                                                                                          \
    }

for (int i = first; i < last; ++i) { body } under the unroll policy of the configuration: the two-branch if constexpr of the other families, the forced-unroll pragma lexically present in the first branch only. The body is the variadic argument.

TDLS_ARROW_LUPP_W

#define TDLS_ARROW_LUPP_W

(

e

)

Value:

work[internal_work ? unsigned(e) : TDLS_STRIDED(e, work_stride)]

Entry e of the workspace: contiguous under internal residency, strided otherwise.

TDLS_ARROW_LUPP_Y

#define TDLS_ARROW_LUPP_Y

(

i

)

Value:

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

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

TDLS_ARROW_LUPP_YW

#define TDLS_ARROW_LUPP_YW	(		m,
			i )

Value:

    y[internal_rhs ? unsigned((ycols == 1 ? 0 : (m)) * dimension + (i))                            \
                   : TDLS_STRIDED(i, rhs_stride) +                                                 \
                         (ycols == 1 ? std::size_t(0) : TDLS_STRIDED_WIDE(m, xcol_stride))]

Entry i of column m of the right-hand sides and solutions: columns of dimension contiguous entries under internal residency, xcol_stride-strided columns otherwise. The function declares ycols, its number of columns: with one, m is ignored and the macro is TDLS_ARROW_LUPP_Y exactly.

Enumeration Type Documentation

ArrowheadRecomputedOutcome

enum class tdls::detail::ArrowheadRecomputedOutcome

strong

Outcome of ArrowheadLUppSolverStatic::solve_inplace_recomputed and solve_inplace_recomputed_or_tangent.

Enumerator
solved	y holds the solution
tangent	converged (solve_inplace_recomputed_or_tangent): x holds the tangent columns, y is untouched
singular	the matrix is singular; y is untouched
deferred	more blocks refused than the capacity of the store; y is untouched, for a solve that stores them (solve_inplace_streamed)

ArrowheadStreamedOutcome

enum class tdls::detail::ArrowheadStreamedOutcome

strong

Outcome of ArrowheadLUppSolverStatic::solve_inplace_streamed_or_tangent.

Enumerator
solved	y holds the solution (converged() returned false)
tangent	converged: x holds the canonical border columns, y intermediate values
singular	the matrix is singular