solver_static.hpp

Go to the documentation of this file.

#ifndef TDLS_SOLVERS_TILED_LUPP_SOLVER_STATIC_HPP
#define TDLS_SOLVERS_TILED_LUPP_SOLVER_STATIC_HPP
 
 
 
 
 
 
#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>
 
 
 
// The flow analysis of older gcc releases (observed on the 11.4 of
// Ubuntu 22.04) cannot prove that the replay row of the out-of-tile
// search is written before being read: L_row[p] is only read for
// p < t, and slot t is written at the end of iteration t, a
// loop-carried property validated by the oracle suites. The resulting
// -Wmaybe-uninitialized reports are spurious.
//
// Under UBSan instrumentation (-fsanitize=undefined) at -O2 and
// above, gcc additionally loses the value ranges of the loop indices
// and reports impossible subscripts through -Warray-bounds: the
// unsigned cast of a provably non-negative index read as a wrapped
// negative, or a loop counter assumed below its own loop guard.
// Observed from gcc 12 to gcc 16; plain builds and clang are clean,
// and the exercised paths are certified free of out-of-bounds
// accesses by the constexpr suite.
//
// Both suppressions are scoped to this header and to those warnings
// only, exactly as in solver_dynamic.hpp.
#if defined(TDLS_DETAIL_COMPILER_GCC)
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
#pragma GCC diagnostic ignored "-Warray-bounds"
#endif
 
 
 
// clang reports a forced unrolling that the optimizer could not perform
// through -Wpass-failed. The unrolling requested by TDLS_UNROLL_FORCE is
// a performance hint: a failed hint does not affect correctness. The
// suppression is scoped to this header and to that warning only.
#if defined(TDLS_DETAIL_COMPILER_CLANG)
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wpass-failed"
#endif
 
namespace tdls {
 
 
 
/* Addressing macros. They expand inside member functions where the
   parameter names (A, A_stride, piv, piv_stride, x, b, y, rhs_stride,
   xcol_stride), the residency template booleans and the Config value
   are in scope. #undef'd at the end of this header. */
 
#define TDLS_LUPP_A(r, c)                                                                          \
    A[internal_matrix ? unsigned((r) * N + (c)) : TDLS_STRIDED((r) * N + (c), A_stride)]
#define TDLS_LUPP_PIV(i) piv[internal_piv ? unsigned(i) : TDLS_STRIDED(i, piv_stride)]
#define TDLS_LUPP_X(i) x[internal_rhs ? unsigned(i) : TDLS_STRIDED(i, rhs_stride)]
#define TDLS_LUPP_B(i) b[internal_rhs ? unsigned(i) : TDLS_STRIDED(i, rhs_stride)]
#define TDLS_LUPP_XW(w, i)                                                                         \
    x[internal_rhs ? unsigned((w) * N + (i))                                                       \
                   : TDLS_STRIDED(i, rhs_stride) + TDLS_STRIDED_WIDE(w, xcol_stride)]
#define TDLS_LUPP_BW(w, i)                                                                         \
    b[internal_rhs ? unsigned((w) * N + (i))                                                       \
                   : TDLS_STRIDED(i, rhs_stride) + TDLS_STRIDED_WIDE(w, xcol_stride)]
#define TDLS_LUPP_Y(i) y[internal_rhs ? unsigned(i) : TDLS_STRIDED(i, rhs_stride)]
#define TDLS_LUPP_FOR(i, first, last, ...)                                                         \
    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__                                                                            \
        }                                                                                          \
    }
 
 
 
template<typename T, int N, TiledLUppConfig<T> Config = TiledLUppConfig<T>{}>
struct TiledLUppSolverStatic {
 
    static constexpr int tile_size = Config.tile_size; 
    static constexpr Schedule schedule =
        Config.schedule; 
 
    static constexpr T oot_pivot_threshold = Config.oot_pivot_threshold;
    static constexpr T relative_pivot_threshold = Config.relative_pivot_threshold;
    static constexpr T singular_floor = Config.singular_floor;
    static constexpr bool transposed = Config.layout == MatrixLayout::ColMajor;
    static constexpr bool physical = Config.row_interchange == RowInterchange::Physical;
 
    static_assert(N >= 1, "TiledLUppSolverStatic: N must be >= 1");
    static_assert(Config.row_interchange != RowInterchange::Auto,
                  "TiledLUppSolverStatic: row_interchange must be Logical or Physical (Auto is "
                  "the default of the blocks of ArrowheadLUpp only)");
    static_assert(
        Config.oot_pivot_threshold.is_finite() && Config.singular_floor.is_finite(),
        "TiledLUppSolverStatic: oot_pivot_threshold and singular_floor must be finite (and fit a "
        "63-bit mantissa)");
    static_assert(Config.relative_pivot_threshold.is_finite(),
                  "TiledLUppSolverStatic: relative_pivot_threshold must be finite (and fit a "
                  "63-bit mantissa)");
    static_assert(relative_pivot_threshold > T(0) && relative_pivot_threshold <= T(1),
                  "TiledLUppSolverStatic: relative_pivot_threshold must be in (0, 1]");
    static_assert(singular_floor > T(0), "TiledLUppSolverStatic: singular_floor must be positive");
    static_assert(singular_floor <= oot_pivot_threshold,
                  "TiledLUppSolverStatic: singular_floor must not exceed oot_pivot_threshold (the "
                  "floor applies to the out-of-tile recovery path)");
    static_assert(tile_size >= 1, "TiledLUppSolverStatic: tile size must be >= 1");
 
    static constexpr int full_tiles = N / tile_size; 
    static constexpr int last_tile_tail = N - full_tiles * tile_size;
    static constexpr int num_tiles = full_tiles + (last_tile_tail > 0 ? 1 : 0);
 
    using Operations = TiledLUppTileOperations<T, tile_size, Config.unroll_loops>;
 
    /* =====================================================================
       Remote <-> register tile movement.
       `load/store_tile` go through a cached physical-row segment (the RL
       schedule keeps one per tile); `load/store_tile_piv` read the
       permutation inline, once per row (LL schedule and the substitution
       do this; it adds no integer registers).
       ===================================================================== */
 
    template<int row_extent, int col_extent, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_tile(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT prow,
              const int col0, T* TDLS_RESTRICT t) noexcept {
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = 0; i < row_extent; ++i) {
                TDLS_UNROLL_FORCE
                for (int j = 0; j < col_extent; ++j)
                    t[i * tile_size + j] = TDLS_LUPP_A(prow[i], col0 + j);
            }
        } else {
            for (int i = 0; i < row_extent; ++i) {
                for (int j = 0; j < col_extent; ++j)
                    t[i * tile_size + j] = TDLS_LUPP_A(prow[i], col0 + j);
            }
        }
    }
 
    template<int row_extent, int col_extent, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    store_tile(T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT prow,
               const int col0, const T* TDLS_RESTRICT t) noexcept {
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = 0; i < row_extent; ++i) {
                TDLS_UNROLL_FORCE
                for (int j = 0; j < col_extent; ++j)
                    TDLS_LUPP_A(prow[i], col0 + j) = t[i * tile_size + j];
            }
        } else {
            for (int i = 0; i < row_extent; ++i) {
                for (int j = 0; j < col_extent; ++j)
                    TDLS_LUPP_A(prow[i], col0 + j) = t[i * tile_size + j];
            }
        }
    }
 
    template<int row_extent, int col_extent, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_tile_piv(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                  const int piv_stride, const int row0, const int col0,
                  T* TDLS_RESTRICT t) noexcept {
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = 0; i < row_extent; ++i) {
                const int phys = TDLS_LUPP_PIV(row0 + i);
                TDLS_UNROLL_FORCE
                for (int j = 0; j < col_extent; ++j)
                    t[i * tile_size + j] = TDLS_LUPP_A(phys, col0 + j);
            }
        } else {
            for (int i = 0; i < row_extent; ++i) {
                const int phys = TDLS_LUPP_PIV(row0 + i);
                for (int j = 0; j < col_extent; ++j)
                    t[i * tile_size + j] = TDLS_LUPP_A(phys, col0 + j);
            }
        }
    }
 
    template<int row_extent, int col_extent, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    store_tile_piv(T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                   const int piv_stride, const int row0, const int col0,
                   const T* TDLS_RESTRICT t) noexcept {
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = 0; i < row_extent; ++i) {
                const int phys = TDLS_LUPP_PIV(row0 + i);
                TDLS_UNROLL_FORCE
                for (int j = 0; j < col_extent; ++j)
                    TDLS_LUPP_A(phys, col0 + j) = t[i * tile_size + j];
            }
        } else {
            for (int i = 0; i < row_extent; ++i) {
                const int phys = TDLS_LUPP_PIV(row0 + i);
                for (int j = 0; j < col_extent; ++j)
                    TDLS_LUPP_A(phys, col0 + j) = t[i * tile_size + j];
            }
        }
    }
 
    template<int row_extent, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_tile_piv_lower(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                        const int piv_stride, const int row0, const int col0,
                        T* TDLS_RESTRICT t) noexcept {
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = 1; i < row_extent; ++i) {
                const int phys = TDLS_LUPP_PIV(row0 + i);
                TDLS_UNROLL_FORCE
                for (int j = 0; j < i; ++j)
                    t[i * tile_size + j] = TDLS_LUPP_A(phys, col0 + j);
            }
        } else {
            for (int i = 1; i < row_extent; ++i) {
                const int phys = TDLS_LUPP_PIV(row0 + i);
                for (int j = 0; j < i; ++j)
                    t[i * tile_size + j] = TDLS_LUPP_A(phys, col0 + j);
            }
        }
    }
 
    template<int row_extent, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_tile_piv_upper(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                        const int piv_stride, const int row0, const int col0,
                        T* TDLS_RESTRICT t) noexcept {
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = 0; i < row_extent; ++i) {
                const int phys = TDLS_LUPP_PIV(row0 + i);
                TDLS_UNROLL_FORCE
                for (int j = i; j < row_extent; ++j)
                    t[i * tile_size + j] = TDLS_LUPP_A(phys, col0 + j);
            }
        } else {
            for (int i = 0; i < row_extent; ++i) {
                const int phys = TDLS_LUPP_PIV(row0 + i);
                for (int j = i; j < row_extent; ++j)
                    t[i * tile_size + j] = TDLS_LUPP_A(phys, col0 + j);
            }
        }
    }
 
    /* =====================================================================
       Diagonal-tile factorization with out-of-tile pivoting.
       Shared by both schedules; the only schedule-dependent part is the
       replay a candidate row needs before it can be compared:
         RL: the trailing matrix is already Schur-updated, so a candidate
             only misses the current tile's columns t < c.
         LL: rows below additionally miss every prior tile's L*U
             contribution (columns [0, k0)), replayed term by term.
       ===================================================================== */
 
    template<int k_extent, bool internal_piv, bool internal_matrix, bool oot_diagnostics,
             bool fuse_rhs, bool internal_rhs>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    factor_diag_column(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
                       const int piv_stride, const int k0, T* TDLS_RESTRICT tile, int& oot_count,
                       const int c, T* TDLS_RESTRICT y, const int rhs_stride) noexcept {
 
        const int gc = k0 + c; // global column
 
        // In-tile pivot search (rows c..k_extent of the register tile)
        int best_r = c;
        T best     = detail::abs(tile[c * tile_size + c]);
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int r = c + 1; r < k_extent; ++r) {
                const T v = detail::abs(tile[r * tile_size + c]);
                if (v > best) {
                    best   = v;
                    best_r = r;
                }
            }
        } else {
            for (int r = c + 1; r < k_extent; ++r) {
                const T v = detail::abs(tile[r * tile_size + c]);
                if (v > best) {
                    best   = v;
                    best_r = r;
                }
            }
        }
 
        // Threshold pivoting inside the tile: the row in place keeps the
        // pivot when it reaches both relative_pivot_threshold times the
        // best in-tile candidate and oot_pivot_threshold. The out-of-tile
        // recovery below keeps its own absolute trigger. A threshold of 1
        // compiles the test out: the best in-tile candidate wins.
        if constexpr (relative_pivot_threshold < T(1)) {
            const T in_place = detail::abs(tile[c * tile_size + c]);
            if (in_place >= relative_pivot_threshold * best && in_place >= oot_pivot_threshold)
                best_r = c;
        }
 
        int piv_row; // winning global (logical) row
 
        if (best >= oot_pivot_threshold) {
            piv_row = k0 + best_r;
        } else if constexpr (k_extent < tile_size) {
            // Trailing tile: no rows below to recover from. Diagnostic
            // order: singularity verdict first, then count the weak pivot
            // (full tiles count before the verdict).
            if (best < singular_floor) return false;
            if constexpr (oot_diagnostics) ++oot_count;
            piv_row = k0 + best_r;
        } else {
            // Out-of-tile recovery: scan the rows below the tile and
            // evaluate each candidate as if it had received the
            // eliminations it is missing, keeping the best (or, with
            // Config.oot_first_acceptable, the first to reach the threshold).
            // Cold path, never unroll-annotated.
            if constexpr (oot_diagnostics) ++oot_count;
            T gbest       = best;
            int gbest_row = k0 + best_r;
 
            for (int row = k0 + k_extent; row < N; ++row) {
                const int phys = TDLS_LUPP_PIV(row);
 
                T corrected = TDLS_LUPP_A(phys, gc);
                if constexpr (schedule == Schedule::LeftLooking) {
                    for (int bj0 = 0; bj0 < k0; bj0 += tile_size)
                        for (int p = 0; p < tile_size; ++p)
                            corrected -= TDLS_LUPP_A(phys, bj0 + p) *
                                         TDLS_LUPP_A(TDLS_LUPP_PIV(bj0 + p), gc);
                }
 
                if (c > 0) {
                    T L_row[tile_size];
                    for (int t = 0; t < c; ++t) {
                        T a_t = TDLS_LUPP_A(phys, k0 + t);
                        if constexpr (schedule == Schedule::LeftLooking) {
                            for (int bj0 = 0; bj0 < k0; bj0 += tile_size)
                                for (int p = 0; p < tile_size; ++p)
                                    a_t -= TDLS_LUPP_A(phys, bj0 + p) *
                                           TDLS_LUPP_A(TDLS_LUPP_PIV(bj0 + p), k0 + t);
                        }
                        for (int p = 0; p < t; ++p)
                            a_t -= L_row[p] * tile[p * tile_size + t];
                        L_row[t] = a_t * tile[t * tile_size + t]; // diag holds 1/pivot
                        corrected -= L_row[t] * tile[t * tile_size + c];
                    }
                }
 
                const T v = detail::abs(corrected);
                if (v > gbest) {
                    gbest     = v;
                    gbest_row = row;
                }
 
                // First-acceptable out-of-tile pivot: a candidate that
                // reaches the threshold already beats the sub-threshold
                // in-tile pivot, so stop scanning. The running max above is
                // kept as the fallback when no candidate is acceptable.
                if constexpr (Config.oot_first_acceptable)
                    if (v >= oot_pivot_threshold) break;
            }
 
            if (gbest < singular_floor) return false;
            piv_row = gbest_row;
        }
 
        {
            // Swap the permutation entries unconditionally: when
            // piv_row == gc the exchanges write back the same values
            // (bitwise no-ops), and removing the comparison removes a
            // divergent branch from the hot path (the lanes of a warp can
            // pick different pivots).
            const int tmp          = TDLS_LUPP_PIV(gc);
            TDLS_LUPP_PIV(gc)      = TDLS_LUPP_PIV(piv_row);
            TDLS_LUPP_PIV(piv_row) = tmp;
 
            if constexpr (fuse_rhs && !transposed) {
                // The fused RHS follows the rows through pivoting. For an
                // internal y the swap is predicated over the unrolled index
                // range: y must never be dynamically indexed or it is
                // demoted to local memory. Under the column-major layout y
                // is indexed by the columns of S, which pivoting never moves.
                if constexpr (internal_rhs) {
                    if constexpr (Config.unroll_loops) {
                        TDLS_UNROLL_FORCE
                        for (int r = 0; r < N; ++r) {
                            if (r == piv_row) {
                                const T ty      = TDLS_LUPP_Y(gc);
                                TDLS_LUPP_Y(gc) = TDLS_LUPP_Y(r);
                                TDLS_LUPP_Y(r)  = ty;
                            }
                        }
                    } else {
                        for (int r = 0; r < N; ++r) {
                            if (r == piv_row) {
                                const T ty      = TDLS_LUPP_Y(gc);
                                TDLS_LUPP_Y(gc) = TDLS_LUPP_Y(r);
                                TDLS_LUPP_Y(r)  = ty;
                            }
                        }
                    }
                } else {
                    const T ty           = TDLS_LUPP_Y(gc);
                    TDLS_LUPP_Y(gc)      = TDLS_LUPP_Y(piv_row);
                    TDLS_LUPP_Y(piv_row) = ty;
                }
            }
 
            if (piv_row < k0 + k_extent) {
                Operations::template swap_rows<k_extent>(tile, c, piv_row - k0);
            } else {
                // Cross-tile swap: pull the new row into the tile and
                // replay everything it missed: prior tiles (LL), then
                // the current tile's factored columns, on the FULL row.
                const int phys = TDLS_LUPP_PIV(gc);
                if constexpr (Config.unroll_loops) {
                    TDLS_UNROLL_FORCE
                    for (int j = 0; j < k_extent; ++j) {
                        tile[c * tile_size + j] = TDLS_LUPP_A(phys, k0 + j);
                        if constexpr (schedule == Schedule::LeftLooking) {
                            for (int bj0 = 0; bj0 < k0; bj0 += tile_size)
                                for (int p = 0; p < tile_size; ++p)
                                    tile[c * tile_size + j] -=
                                        TDLS_LUPP_A(phys, bj0 + p) *
                                        TDLS_LUPP_A(TDLS_LUPP_PIV(bj0 + p), k0 + j);
                        }
                    }
                } else {
                    for (int j = 0; j < k_extent; ++j) {
                        tile[c * tile_size + j] = TDLS_LUPP_A(phys, k0 + j);
                        if constexpr (schedule == Schedule::LeftLooking) {
                            for (int bj0 = 0; bj0 < k0; bj0 += tile_size)
                                for (int p = 0; p < tile_size; ++p)
                                    tile[c * tile_size + j] -=
                                        TDLS_LUPP_A(phys, bj0 + p) *
                                        TDLS_LUPP_A(TDLS_LUPP_PIV(bj0 + p), k0 + j);
                        }
                    }
                }
 
                if (c > 0) {
                    T L_row[tile_size];
                    for (int t = 0; t < c; ++t) {
                        T a_t = tile[c * tile_size + t];
                        for (int p = 0; p < t; ++p)
                            a_t -= L_row[p] * tile[p * tile_size + t];
                        L_row[t] = a_t * tile[t * tile_size + t]; // diag holds 1/pivot
                        tile[c * tile_size + t] = L_row[t];
                    }
                    for (int j = c; j < k_extent; ++j) {
                        for (int t = 0; t < c; ++t)
                            tile[c * tile_size + j] -= L_row[t] * tile[t * tile_size + j];
                    }
                }
            }
        }
 
        Operations::template eliminate_column<k_extent, k_extent>(tile, c);
        return true;
    }
 
    template<int k_extent, bool internal_piv, bool internal_matrix, bool oot_diagnostics,
             bool fuse_rhs = false, bool internal_rhs = true>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    factor_diag_tile(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
                     const int piv_stride, const int k0, T* TDLS_RESTRICT tile, int& oot_count,
                     T* TDLS_RESTRICT y = nullptr, const int rhs_stride = 1) noexcept {
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int c = 0; c < k_extent; ++c) {
                if (!factor_diag_column<k_extent, internal_piv, internal_matrix, oot_diagnostics,
                                        fuse_rhs, internal_rhs>(A, A_stride, piv, piv_stride, k0,
                                                                tile, oot_count, c, y, rhs_stride))
                    return false;
            }
        } else {
            for (int c = 0; c < k_extent; ++c) {
                if (!factor_diag_column<k_extent, internal_piv, internal_matrix, oot_diagnostics,
                                        fuse_rhs, internal_rhs>(A, A_stride, piv, piv_stride, k0,
                                                                tile, oot_count, c, y, rhs_stride))
                    return false;
            }
        }
        return true;
    }
 
    /* =====================================================================
       RIGHT-LOOKING schedule. Step k: factor the diagonal tile, then push
       its factors into the trailing matrix (TRSM right over the row panel,
       TRSM down + Schur complement over the rows below). The physical-row
       segments pk/pi are cached per tile. Their entries are known at run
       time only, an internal permutation included: a caller-local matrix
       is indexed through them, so on GPU it lives in local memory.
       ===================================================================== */
 
    template<int k_extent, int j_extent, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    rl_trsm_right_one(T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT pk,
                      const T* TDLS_RESTRICT tile, const int j0) noexcept {
        T Akj[tile_size * tile_size];
        load_tile<k_extent, j_extent, internal_matrix>(A, A_stride, pk, j0, Akj);
        Operations::template trsm_left_unit<k_extent, j_extent>(tile, Akj);
        store_tile<k_extent, j_extent, internal_matrix>(A, A_stride, pk, j0, Akj);
    }
 
    template<int k_extent, int j_extent, bool internal_matrix, bool internal_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    rl_trsm_right_push_one(T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT pk,
                           const T* TDLS_RESTRICT tile, const int j0, const int k0,
                           T* TDLS_RESTRICT y, const int rhs_stride) noexcept {
        T Akj[tile_size * tile_size];
        load_tile<k_extent, j_extent, internal_matrix>(A, A_stride, pk, j0, Akj);
        Operations::template trsm_left_unit<k_extent, j_extent>(tile, Akj);
        store_tile<k_extent, j_extent, internal_matrix>(A, A_stride, pk, j0, Akj);
        transposed_push_one<k_extent, j_extent, 1, internal_rhs>(Akj, y, rhs_stride, 0, k0, j0);
    }
 
    template<int k_extent, int i_extent, int j_extent, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    rl_schur_one(T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT pk,
                 const int* TDLS_RESTRICT pi, const T* TDLS_RESTRICT Aik, const int j0) noexcept {
        T Aij[tile_size * tile_size];
        load_tile<i_extent, j_extent, internal_matrix>(A, A_stride, pi, j0, Aij);
 
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int p = 0; p < k_extent; ++p) {
                T Akj_row[tile_size];
                TDLS_UNROLL_FORCE
                for (int j = 0; j < j_extent; ++j)
                    Akj_row[j] = TDLS_LUPP_A(pk[p], j0 + j);
                TDLS_UNROLL_FORCE
                for (int i = 0; i < i_extent; ++i) {
                    const T L_ip = Aik[i * tile_size + p];
                    TDLS_UNROLL_FORCE
                    for (int j = 0; j < j_extent; ++j)
                        Aij[i * tile_size + j] -= L_ip * Akj_row[j];
                }
            }
        } else {
            for (int p = 0; p < k_extent; ++p) {
                T Akj_row[tile_size];
                for (int j = 0; j < j_extent; ++j)
                    Akj_row[j] = TDLS_LUPP_A(pk[p], j0 + j);
                for (int i = 0; i < i_extent; ++i) {
                    const T L_ip = Aik[i * tile_size + p];
                    for (int j = 0; j < j_extent; ++j)
                        Aij[i * tile_size + j] -= L_ip * Akj_row[j];
                }
            }
        }
 
        store_tile<i_extent, j_extent, internal_matrix>(A, A_stride, pi, j0, Aij);
    }
 
    template<int k_extent, int i_extent, bool internal_piv, bool internal_matrix,
             bool fuse_rhs = false, bool internal_rhs = true>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    rl_update_row_one(T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                      const int piv_stride, const int* TDLS_RESTRICT pk,
                      const T* TDLS_RESTRICT tile, const int k, const int i0,
                      T* TDLS_RESTRICT y = nullptr, const int rhs_stride = 1) noexcept {
        const int k0 = k * tile_size;
 
        int pi[tile_size];
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = 0; i < i_extent; ++i)
                pi[i] = TDLS_LUPP_PIV(i0 + i);
        } else {
            for (int i = 0; i < i_extent; ++i)
                pi[i] = TDLS_LUPP_PIV(i0 + i);
        }
 
        // TRSM down: Aik := Aik * U^-1
        T Aik[tile_size * tile_size];
        load_tile<i_extent, k_extent, internal_matrix>(A, A_stride, pi, k0, Aik);
        Operations::template trsm_right<k_extent, i_extent>(tile, Aik);
        store_tile<i_extent, k_extent, internal_matrix>(A, A_stride, pi, k0, Aik);
 
        // Fused forward substitution: push the solved y_k segment into this
        // row block while its L panel sits in registers (this is the whole
        // point of solve_inplace: the separate forward pass reloads vanish).
        // Under the column-major layout the forward pass runs on U^T and is
        // pushed from the U row panel instead (rl_trsm_right_push_one).
        if constexpr (fuse_rhs && !transposed) {
            if constexpr (Config.unroll_loops) {
                TDLS_UNROLL_FORCE
                for (int r = 0; r < i_extent; ++r) {
                    T sum = T(0);
                    TDLS_UNROLL_FORCE
                    for (int j = 0; j < k_extent; ++j)
                        sum += Aik[r * tile_size + j] * TDLS_LUPP_Y(k0 + j);
                    TDLS_LUPP_Y(i0 + r) -= sum;
                }
            } else {
                for (int r = 0; r < i_extent; ++r) {
                    T sum = T(0);
                    for (int j = 0; j < k_extent; ++j)
                        sum += Aik[r * tile_size + j] * TDLS_LUPP_Y(k0 + j);
                    TDLS_LUPP_Y(i0 + r) -= sum;
                }
            }
        }
 
        // Schur sweep over the trailing columns
        for (int j = k + 1; j < full_tiles; ++j)
            rl_schur_one<k_extent, i_extent, tile_size, internal_matrix>(A, A_stride, pk, pi, Aik,
                                                                         j * tile_size);
        if constexpr (last_tile_tail > 0 && k_extent == tile_size)
            rl_schur_one<k_extent, i_extent, last_tile_tail, internal_matrix>(
                A, A_stride, pk, pi, Aik, full_tiles * tile_size);
    }
 
    template<int k_extent, bool internal_piv, bool internal_matrix, bool oot_diagnostics,
             bool fuse_rhs = false, bool internal_rhs = true>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    rl_step(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv, const int piv_stride,
            const int k, int& oot_count, T* TDLS_RESTRICT y = nullptr,
            const int rhs_stride = 1) noexcept {
        const int k0 = k * tile_size;
 
        T tile[tile_size * tile_size];
        load_tile_piv<k_extent, k_extent, internal_piv, internal_matrix>(A, A_stride, piv,
                                                                         piv_stride, k0, k0, tile);
 
        if (!factor_diag_tile<k_extent, internal_piv, internal_matrix, oot_diagnostics, fuse_rhs,
                              internal_rhs>(A, A_stride, piv, piv_stride, k0, tile, oot_count, y,
                                            rhs_stride))
            return false;
 
        // Physical rows of the tile after the swaps of this step
        int pk[tile_size];
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = 0; i < k_extent; ++i)
                pk[i] = TDLS_LUPP_PIV(k0 + i);
        } else {
            for (int i = 0; i < k_extent; ++i)
                pk[i] = TDLS_LUPP_PIV(k0 + i);
        }
 
        store_tile<k_extent, k_extent, internal_matrix>(A, A_stride, pk, k0, tile);
 
        // Fused forward substitution: this tile's y segment is final from
        // here on: solve it while the tile is in registers, unit-lower
        // (row-major layout) or with U^T (column-major layout).
        if constexpr (fuse_rhs) {
            if constexpr (transposed) {
                transposed_solve_tile_upper<k_extent, 1, internal_rhs>(tile, y, rhs_stride, 0, k0);
            } else if constexpr (Config.unroll_loops) {
                TDLS_UNROLL_FORCE
                for (int kk = 0; kk < k_extent; ++kk) {
                    TDLS_UNROLL_FORCE
                    for (int i = kk + 1; i < k_extent; ++i)
                        TDLS_LUPP_Y(k0 + i) -= tile[i * tile_size + kk] * TDLS_LUPP_Y(k0 + kk);
                }
            } else {
                for (int kk = 0; kk < k_extent; ++kk) {
                    for (int i = kk + 1; i < k_extent; ++i)
                        TDLS_LUPP_Y(k0 + i) -= tile[i * tile_size + kk] * TDLS_LUPP_Y(k0 + kk);
                }
            }
        }
 
        // TRSM right over the row panel; under the column-major layout,
        // solve_inplace pushes the solved z_k segment along
        for (int j = k + 1; j < full_tiles; ++j) {
            if constexpr (fuse_rhs && transposed)
                rl_trsm_right_push_one<k_extent, tile_size, internal_matrix, internal_rhs>(
                    A, A_stride, pk, tile, j * tile_size, k0, y, rhs_stride);
            else
                rl_trsm_right_one<k_extent, tile_size, internal_matrix>(A, A_stride, pk, tile,
                                                                        j * tile_size);
        }
        if constexpr (last_tile_tail > 0 && k_extent == tile_size) {
            if constexpr (fuse_rhs && transposed)
                rl_trsm_right_push_one<k_extent, last_tile_tail, internal_matrix, internal_rhs>(
                    A, A_stride, pk, tile, full_tiles * tile_size, k0, y, rhs_stride);
            else
                rl_trsm_right_one<k_extent, last_tile_tail, internal_matrix>(
                    A, A_stride, pk, tile, full_tiles * tile_size);
        }
 
        // TRSM down + Schur over the rows below
        for (int i = k + 1; i < full_tiles; ++i)
            rl_update_row_one<k_extent, tile_size, internal_piv, internal_matrix, fuse_rhs,
                              internal_rhs>(A, A_stride, piv, piv_stride, pk, tile, k,
                                            i * tile_size, y, rhs_stride);
        if constexpr (last_tile_tail > 0 && k_extent == tile_size)
            rl_update_row_one<k_extent, last_tile_tail, internal_piv, internal_matrix, fuse_rhs,
                              internal_rhs>(A, A_stride, piv, piv_stride, pk, tile, k,
                                            full_tiles * tile_size, y, rhs_stride);
        return true;
    }
 
    /* =====================================================================
       LEFT-LOOKING schedule. Step k: pull every prior tile's L*U
       contribution into the current block column/row, then factor the
       diagonal tile and triangular-solve its panels. Each tile of the
       trailing matrix is written once per factorization (vs once per step
       for RL) at the cost of replaying prior tiles on each visit. The
       permutation is read inline (one read per row), never cached in
       segments.
       ===================================================================== */
 
    template<int row_extent, int col_extent, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    ll_correct_tile(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                    const int piv_stride, const int row0, const int col0, const int k0,
                    T* TDLS_RESTRICT t) noexcept {
        for (int bj0 = 0; bj0 < k0; bj0 += tile_size) {
            T Lt[tile_size * tile_size];
            load_tile_piv<row_extent, tile_size, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, row0, bj0, Lt);
            T Ut[tile_size * tile_size];
            load_tile_piv<tile_size, col_extent, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, bj0, col0, Ut);
            Operations::template gemm_sub<row_extent, col_extent, tile_size>(t, Lt, Ut);
        }
    }
 
    template<int k_extent, int i_extent, bool internal_piv, bool internal_matrix,
             bool fuse_rhs = false, bool internal_rhs = true>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    ll_update_below_one(T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                        const int piv_stride, const T* TDLS_RESTRICT tile, const int k0,
                        const int i0, T* TDLS_RESTRICT y = nullptr,
                        const int rhs_stride = 1) noexcept {
        T B[tile_size * tile_size];
        load_tile_piv<i_extent, k_extent, internal_piv, internal_matrix>(A, A_stride, piv,
                                                                         piv_stride, i0, k0, B);
        ll_correct_tile<i_extent, k_extent, internal_piv, internal_matrix>(
            A, A_stride, piv, piv_stride, i0, k0, k0, B);
        Operations::template trsm_right<k_extent, i_extent>(tile, B);
        store_tile_piv<i_extent, k_extent, internal_piv, internal_matrix>(A, A_stride, piv,
                                                                          piv_stride, i0, k0, B);
 
        // Fused forward substitution: B is the final L(i,k) panel, so push the
        // solved y_k segment into this row block while it is in registers.
        // Under the column-major layout the forward pass runs on U^T and is
        // pushed from the U row panel instead (ll_update_right_push_one).
        if constexpr (fuse_rhs && !transposed) {
            if constexpr (Config.unroll_loops) {
                TDLS_UNROLL_FORCE
                for (int r = 0; r < i_extent; ++r) {
                    T sum = T(0);
                    TDLS_UNROLL_FORCE
                    for (int j = 0; j < k_extent; ++j)
                        sum += B[r * tile_size + j] * TDLS_LUPP_Y(k0 + j);
                    TDLS_LUPP_Y(i0 + r) -= sum;
                }
            } else {
                for (int r = 0; r < i_extent; ++r) {
                    T sum = T(0);
                    for (int j = 0; j < k_extent; ++j)
                        sum += B[r * tile_size + j] * TDLS_LUPP_Y(k0 + j);
                    TDLS_LUPP_Y(i0 + r) -= sum;
                }
            }
        }
    }
 
    template<int k_extent, int j_extent, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    ll_update_right_one(T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                        const int piv_stride, const T* TDLS_RESTRICT tile, const int k0,
                        const int j0) noexcept {
        T B[tile_size * tile_size];
        load_tile_piv<k_extent, j_extent, internal_piv, internal_matrix>(A, A_stride, piv,
                                                                         piv_stride, k0, j0, B);
        ll_correct_tile<k_extent, j_extent, internal_piv, internal_matrix>(
            A, A_stride, piv, piv_stride, k0, j0, k0, B);
        Operations::template trsm_left_unit<k_extent, j_extent>(tile, B);
        store_tile_piv<k_extent, j_extent, internal_piv, internal_matrix>(A, A_stride, piv,
                                                                          piv_stride, k0, j0, B);
    }
 
    template<int k_extent, int j_extent, bool internal_piv, bool internal_matrix, bool internal_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    ll_update_right_push_one(T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                             const int piv_stride, const T* TDLS_RESTRICT tile, const int k0,
                             const int j0, T* TDLS_RESTRICT y, const int rhs_stride) noexcept {
        T B[tile_size * tile_size];
        load_tile_piv<k_extent, j_extent, internal_piv, internal_matrix>(A, A_stride, piv,
                                                                         piv_stride, k0, j0, B);
        ll_correct_tile<k_extent, j_extent, internal_piv, internal_matrix>(
            A, A_stride, piv, piv_stride, k0, j0, k0, B);
        Operations::template trsm_left_unit<k_extent, j_extent>(tile, B);
        store_tile_piv<k_extent, j_extent, internal_piv, internal_matrix>(A, A_stride, piv,
                                                                          piv_stride, k0, j0, B);
        transposed_push_one<k_extent, j_extent, 1, internal_rhs>(B, y, rhs_stride, 0, k0, j0);
    }
 
    template<int k_extent, bool internal_piv, bool internal_matrix, bool oot_diagnostics,
             bool fuse_rhs = false, bool internal_rhs = true>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    ll_step(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv, const int piv_stride,
            const int k, int& oot_count, T* TDLS_RESTRICT y = nullptr,
            const int rhs_stride = 1) noexcept {
        const int k0 = k * tile_size;
 
        T tile[tile_size * tile_size];
        load_tile_piv<k_extent, k_extent, internal_piv, internal_matrix>(A, A_stride, piv,
                                                                         piv_stride, k0, k0, tile);
        ll_correct_tile<k_extent, k_extent, internal_piv, internal_matrix>(
            A, A_stride, piv, piv_stride, k0, k0, k0, tile);
 
        if (!factor_diag_tile<k_extent, internal_piv, internal_matrix, oot_diagnostics, fuse_rhs,
                              internal_rhs>(A, A_stride, piv, piv_stride, k0, tile, oot_count, y,
                                            rhs_stride))
            return false;
 
        store_tile_piv<k_extent, k_extent, internal_piv, internal_matrix>(A, A_stride, piv,
                                                                          piv_stride, k0, k0, tile);
 
        // Fused forward substitution: this tile's y segment is final from
        // here on: solve it while the tile is in registers, unit-lower
        // (row-major layout) or with U^T (column-major layout).
        if constexpr (fuse_rhs) {
            if constexpr (transposed) {
                transposed_solve_tile_upper<k_extent, 1, internal_rhs>(tile, y, rhs_stride, 0, k0);
            } else if constexpr (Config.unroll_loops) {
                TDLS_UNROLL_FORCE
                for (int kk = 0; kk < k_extent; ++kk) {
                    TDLS_UNROLL_FORCE
                    for (int i = kk + 1; i < k_extent; ++i)
                        TDLS_LUPP_Y(k0 + i) -= tile[i * tile_size + kk] * TDLS_LUPP_Y(k0 + kk);
                }
            } else {
                for (int kk = 0; kk < k_extent; ++kk) {
                    for (int i = kk + 1; i < k_extent; ++i)
                        TDLS_LUPP_Y(k0 + i) -= tile[i * tile_size + kk] * TDLS_LUPP_Y(k0 + kk);
                }
            }
        }
 
        // L panel below the diagonal
        for (int i = k + 1; i < full_tiles; ++i)
            ll_update_below_one<k_extent, tile_size, internal_piv, internal_matrix, fuse_rhs,
                                internal_rhs>(A, A_stride, piv, piv_stride, tile, k0, i * tile_size,
                                              y, rhs_stride);
        if constexpr (last_tile_tail > 0 && k_extent == tile_size)
            ll_update_below_one<k_extent, last_tile_tail, internal_piv, internal_matrix, fuse_rhs,
                                internal_rhs>(A, A_stride, piv, piv_stride, tile, k0,
                                              full_tiles * tile_size, y, rhs_stride);
 
        // U row panel right of the diagonal; under the column-major layout,
        // solve_inplace pushes the solved z_k segment along
        for (int j = k + 1; j < full_tiles; ++j) {
            if constexpr (fuse_rhs && transposed)
                ll_update_right_push_one<k_extent, tile_size, internal_piv, internal_matrix,
                                         internal_rhs>(A, A_stride, piv, piv_stride, tile, k0,
                                                       j * tile_size, y, rhs_stride);
            else
                ll_update_right_one<k_extent, tile_size, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, tile, k0, j * tile_size);
        }
        if constexpr (last_tile_tail > 0 && k_extent == tile_size) {
            if constexpr (fuse_rhs && transposed)
                ll_update_right_push_one<k_extent, last_tile_tail, internal_piv, internal_matrix,
                                         internal_rhs>(A, A_stride, piv, piv_stride, tile, k0,
                                                       full_tiles * tile_size, y, rhs_stride);
            else
                ll_update_right_one<k_extent, last_tile_tail, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, tile, k0, full_tiles * tile_size);
        }
        return true;
    }
 
    /* =====================================================================
       Physical row interchanges (Config.row_interchange == Physical): the
       whole matrix is one tile, the rows exchanged as in LAPACK. Every
       loop runs over a full range with the triangle selected inside: nvcc
       does not unroll a loop whose trip count depends on an outer index,
       and a run-time index would send a caller-local array to the stack.
       ===================================================================== */
 
    template<int w, bool internal_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    physical_interchange_rhs(T* TDLS_RESTRICT x, const int rhs_stride,
                             [[maybe_unused]] const std::size_t xcol_stride, const int j,
                             const int r) noexcept {
        TDLS_LUPP_FOR(
            i, 0, N, if (i > j && i == r) {
                TDLS_LUPP_FOR(c, 0, w, const T t = TDLS_LUPP_XW(c, j);
                              TDLS_LUPP_XW(c, j) = TDLS_LUPP_XW(c, i); TDLS_LUPP_XW(c, i) = t;)
            })
    }
 
    template<int n>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    physical_interchange_local(T* TDLS_RESTRICT X, const int j, const int r, const int cs,
                               const int rs) noexcept {
        TDLS_LUPP_FOR(
            i, 0, N, if (i > j && i == r) {
                TDLS_LUPP_FOR(c, 0, n, const T t = X[j * rs + c * cs];
                              X[j * rs + c * cs] = X[i * rs + c * cs]; X[i * rs + c * cs] = t;)
            })
    }
 
    template<int w, bool internal_rhs, bool internal_piv, bool reverse>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    physical_permute(const int* TDLS_RESTRICT piv, [[maybe_unused]] const int piv_stride,
                     T* TDLS_RESTRICT x, const int rhs_stride,
                     const std::size_t xcol_stride) noexcept {
        TDLS_LUPP_FOR(
            t, 0, N, const int j = reverse ? N - 1 - t : t; if (TDLS_LUPP_PIV(j) != j) {
                if constexpr (internal_rhs)
                    physical_interchange_local<w>(x, j, TDLS_LUPP_PIV(j), N, 1);
                else
                    physical_interchange_rhs<w, internal_rhs>(x, rhs_stride, xcol_stride, j,
                                                              TDLS_LUPP_PIV(j));
            })
    }
 
    template<bool internal_piv, bool internal_matrix, bool check>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    physical_step(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
                  const int piv_stride, const int j, [[maybe_unused]] const T u, bool& singular,
                  int& weak) noexcept {
        T best = detail::abs(TDLS_LUPP_A(j, j));
        int r  = j;
        TDLS_LUPP_FOR(
            i, 0, N, const T a = detail::abs(TDLS_LUPP_A(i, j)); if (i > j && a > best) {
                best = a;
                r    = i;
            })
        if constexpr (relative_pivot_threshold < T(1)) {
            const T in_place = detail::abs(TDLS_LUPP_A(j, j));
            if (in_place >= relative_pivot_threshold * best && in_place >= oot_pivot_threshold)
                r = j;
        }
        const bool pivot_ok = best >= singular_floor;
        singular            = singular || !pivot_ok;
        weak += best < oot_pivot_threshold ? 1 : 0;
        TDLS_LUPP_PIV(j) = r;
        // Rarely taken on matrices that keep their rows: uniform across a warp.
        if (r != j) {
            TDLS_LUPP_FOR(
                i, 0, N, if (i > j && i == r) {
                    TDLS_LUPP_FOR(c, 0, N, const T t = TDLS_LUPP_A(j, c);
                                  TDLS_LUPP_A(j, c) = TDLS_LUPP_A(i, c); TDLS_LUPP_A(i, c) = t;)
                })
        }
        const T inv       = pivot_ok ? T(1) / TDLS_LUPP_A(j, j) : T(0);
        TDLS_LUPP_A(j, j) = inv;
        bool ok           = true;
        TDLS_LUPP_FOR(
            i, 0, N, if (i > j) {
                const T l         = TDLS_LUPP_A(i, j) * inv;
                TDLS_LUPP_A(i, j) = l;
                if constexpr (check) ok = ok && u * detail::abs(l) <= T(1);
                TDLS_LUPP_FOR(c, 0, N, if (c > j) TDLS_LUPP_A(i, c) -= l * TDLS_LUPP_A(j, c);)
            })
        return ok;
    }
 
    template<bool internal_piv, bool internal_matrix>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    physical_factorize(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
                       const int piv_stride, int& weak) noexcept {
        bool singular = false;
        weak          = 0;
        TDLS_LUPP_FOR(j, 0, N,
                      (void)physical_step<internal_piv, internal_matrix, false>(
                          A, A_stride, piv, piv_stride, j, T(0), singular, weak);)
        return !singular;
    }
 
    template<bool internal_piv, bool internal_matrix>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr int
    physical_factorize_checked(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
                               const int piv_stride, const T u) noexcept {
        bool ok = true, singular = false;
        int weak = 0;
        TDLS_LUPP_FOR(j, 0, N,
                      ok = physical_step<internal_piv, internal_matrix, true>(
                               A, A_stride, piv, piv_stride, j, u, singular, weak) &&
                           ok;)
        return singular ? 2 : ok ? 0 : 1;
    }
 
    template<int w, bool internal_rhs, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    physical_forward(const T* TDLS_RESTRICT A, const int A_stride,
                     [[maybe_unused]] const int* TDLS_RESTRICT piv,
                     [[maybe_unused]] const int piv_stride, T* TDLS_RESTRICT x,
                     const int rhs_stride, const std::size_t xcol_stride) noexcept {
        if constexpr (!transposed) {
            physical_permute<w, internal_rhs, internal_piv, false>(piv, piv_stride, x, rhs_stride,
                                                                   xcol_stride);
            TDLS_LUPP_FOR(j, 0, N,
                          TDLS_LUPP_FOR(
                              i, 0, N, if (i > j) {
                                  const T l = TDLS_LUPP_A(i, j);
                                  TDLS_LUPP_FOR(c, 0, w,
                                                TDLS_LUPP_XW(c, i) -= l * TDLS_LUPP_XW(c, j);)
                              }))
        } else {
            TDLS_LUPP_FOR(
                j, 0, N,
                TDLS_LUPP_FOR(c, 0, w, TDLS_LUPP_XW(c, j) *= TDLS_LUPP_A(j, j);) TDLS_LUPP_FOR(
                    q, 0, N, if (q > j) {
                        const T uq = TDLS_LUPP_A(j, q);
                        TDLS_LUPP_FOR(c, 0, w, TDLS_LUPP_XW(c, q) -= uq * TDLS_LUPP_XW(c, j);)
                    }))
        }
    }
 
    template<int w, bool internal_rhs, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    physical_backward(const T* TDLS_RESTRICT A, const int A_stride,
                      [[maybe_unused]] const int* TDLS_RESTRICT piv,
                      [[maybe_unused]] const int piv_stride, T* TDLS_RESTRICT x,
                      const int rhs_stride, const std::size_t xcol_stride) noexcept {
        if constexpr (!transposed) {
            TDLS_LUPP_FOR(
                t, 0, N, const int j = N - 1 - t; TDLS_LUPP_FOR(
                    c, 0, w, T v = TDLS_LUPP_XW(c, j);
                    TDLS_LUPP_FOR(q, 0, N, if (q > j) v -= TDLS_LUPP_A(j, q) * TDLS_LUPP_XW(c, q);)
                        TDLS_LUPP_XW(c, j) = v * TDLS_LUPP_A(j, j);))
        } else {
            TDLS_LUPP_FOR(t, 0, N, const int j = N - 1 - t; TDLS_LUPP_FOR(
                              q, 0, N, if (q > j) {
                                  const T l = TDLS_LUPP_A(q, j);
                                  TDLS_LUPP_FOR(c, 0, w,
                                                TDLS_LUPP_XW(c, j) -= l * TDLS_LUPP_XW(c, q);)
                              }))
            physical_permute<w, internal_rhs, internal_piv, true>(piv, piv_stride, x, rhs_stride,
                                                                  xcol_stride);
        }
    }
 
    template<int nrhs, int pass_width, int stages, bool internal_rhs, bool internal_piv,
             bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    physical_multirhs(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                      const int piv_stride, T* TDLS_RESTRICT x, const int rhs_stride,
                      const std::size_t xcol_stride) noexcept {
        constexpr int cpp    = (pass_width <= 0 || pass_width >= nrhs) ? nrhs : pass_width;
        constexpr int passes = nrhs / cpp;
        constexpr int tail   = nrhs % cpp;
        for (int p = 0; p < passes; ++p) {
            const int c0 = p * cpp;
            const std::size_t off =
                internal_rhs ? unsigned(c0) * unsigned(N) : TDLS_STRIDED_WIDE(c0, xcol_stride);
            if constexpr ((stages & 1) != 0)
                physical_forward<cpp, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x + off, rhs_stride, xcol_stride);
            if constexpr ((stages & 2) != 0)
                physical_backward<cpp, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x + off, rhs_stride, xcol_stride);
        }
        if constexpr (tail > 0) {
            constexpr int c0 = nrhs - tail;
            const std::size_t off =
                internal_rhs ? unsigned(c0) * unsigned(N) : TDLS_STRIDED_WIDE(c0, xcol_stride);
            if constexpr ((stages & 1) != 0)
                physical_forward<tail, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x + off, rhs_stride, xcol_stride);
            if constexpr ((stages & 2) != 0)
                physical_backward<tail, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x + off, rhs_stride, xcol_stride);
        }
    }
 
    /* =====================================================================
       FACTORIZE - public entry point.
       ===================================================================== */
 
    template<bool internal_piv, bool internal_matrix, bool oot_diagnostics = true,
             bool fuse_rhs = false, bool internal_rhs = true, int A_static_stride = dynamic_stride,
             int piv_static_stride = dynamic_stride, int rhs_static_stride = dynamic_stride>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    factorize(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv, const int piv_stride,
              int& oot_count, T* TDLS_RESTRICT y = nullptr, const int rhs_stride = 1) noexcept {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return factorize<internal_piv, internal_matrix, oot_diagnostics, fuse_rhs,
                             internal_rhs>(A, detail::stride_or<A_static_stride>(A_stride), piv,
                                           detail::stride_or<piv_static_stride>(piv_stride),
                                           oot_count, y,
                                           detail::stride_or<rhs_static_stride>(rhs_stride));
        } else {
            if constexpr (physical) {
                int weak      = 0;
                const bool ok = physical_factorize<internal_piv, internal_matrix>(A, A_stride, piv,
                                                                                  piv_stride, weak);
                if constexpr (oot_diagnostics) oot_count = weak;
                if (!ok) return false;
                if constexpr (fuse_rhs)
                    physical_forward<1, internal_rhs, internal_piv, internal_matrix>(
                        A, A_stride, piv, piv_stride, y, rhs_stride, 0);
                return true;
            }
 
            if constexpr (oot_diagnostics) oot_count = 0;
 
            if constexpr (Config.unroll_loops) {
                TDLS_UNROLL_FORCE
                for (int i = 0; i < N; ++i)
                    TDLS_LUPP_PIV(i) = i;
            } else {
                for (int i = 0; i < N; ++i)
                    TDLS_LUPP_PIV(i) = i;
            }
 
            if constexpr (schedule == Schedule::RightLooking) {
                for (int k = 0; k < full_tiles; ++k)
                    if (!rl_step<tile_size, internal_piv, internal_matrix, oot_diagnostics,
                                 fuse_rhs, internal_rhs>(A, A_stride, piv, piv_stride, k, oot_count,
                                                         y, rhs_stride))
                        return false;
                if constexpr (last_tile_tail > 0)
                    if (!rl_step<last_tile_tail, internal_piv, internal_matrix, oot_diagnostics,
                                 fuse_rhs, internal_rhs>(A, A_stride, piv, piv_stride, full_tiles,
                                                         oot_count, y, rhs_stride))
                        return false;
            } else {
                for (int k = 0; k < full_tiles; ++k)
                    if (!ll_step<tile_size, internal_piv, internal_matrix, oot_diagnostics,
                                 fuse_rhs, internal_rhs>(A, A_stride, piv, piv_stride, k, oot_count,
                                                         y, rhs_stride))
                        return false;
                if constexpr (last_tile_tail > 0)
                    if (!ll_step<last_tile_tail, internal_piv, internal_matrix, oot_diagnostics,
                                 fuse_rhs, internal_rhs>(A, A_stride, piv, piv_stride, full_tiles,
                                                         oot_count, y, rhs_stride))
                        return false;
            }
 
            return true;
        }
    }
 
    template<bool internal_piv, bool internal_matrix, int A_static_stride = dynamic_stride,
             int piv_static_stride = dynamic_stride>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    factorize(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
              const int piv_stride) noexcept {
        int unused = 0;
        return factorize<internal_piv, internal_matrix, false, false, true, A_static_stride,
                         piv_static_stride>(A, A_stride, piv, piv_stride, unused);
    }
 
    /* =====================================================================
       SUBSTITUTION - schedule-independent.
       Forward: unit-lower solve tile by tile, each solved segment pushed
       into the tiles below it. Backward: each segment first pulls the
       trailing contributions, then upper-solves its diagonal tile.
       ===================================================================== */
 
    template<int k_extent, int m_extent, int pass_width, bool internal_rhs, bool internal_piv,
             bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    fwd_push_one(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                 const int piv_stride, T* TDLS_RESTRICT x, const int rhs_stride,
                 const std::size_t xcol_stride, const int k0, const int m0) noexcept {
        T Lmk[tile_size * tile_size];
        load_tile_piv<m_extent, k_extent, internal_piv, internal_matrix>(A, A_stride, piv,
                                                                         piv_stride, m0, k0, Lmk);
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = 0; i < m_extent; ++i) {
                TDLS_UNROLL_FORCE
                for (int w = 0; w < pass_width; ++w) {
                    T sum = T(0);
                    TDLS_UNROLL_FORCE
                    for (int j = 0; j < k_extent; ++j)
                        sum += Lmk[i * tile_size + j] * TDLS_LUPP_XW(w, k0 + j);
                    TDLS_LUPP_XW(w, m0 + i) -= sum;
                }
            }
        } else {
            for (int i = 0; i < m_extent; ++i) {
                for (int w = 0; w < pass_width; ++w) {
                    T sum = T(0);
                    for (int j = 0; j < k_extent; ++j)
                        sum += Lmk[i * tile_size + j] * TDLS_LUPP_XW(w, k0 + j);
                    TDLS_LUPP_XW(w, m0 + i) -= sum;
                }
            }
        }
    }
 
    template<int k_extent, int pass_width, bool internal_rhs, bool internal_piv,
             bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    fwd_step(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
             const int piv_stride, T* TDLS_RESTRICT x, const int rhs_stride,
             const std::size_t xcol_stride, const int k) noexcept {
        const int k0 = k * tile_size;
 
        T Lkk[tile_size * tile_size];
        load_tile_piv_lower<k_extent, internal_piv, internal_matrix>(A, A_stride, piv, piv_stride,
                                                                     k0, k0, Lkk);
 
        // In-tile unit-lower solve
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int kk = 0; kk < k_extent; ++kk) {
                TDLS_UNROLL_FORCE
                for (int i = kk + 1; i < k_extent; ++i) {
                    TDLS_UNROLL_FORCE
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, k0 + i) -=
                            Lkk[i * tile_size + kk] * TDLS_LUPP_XW(w, k0 + kk);
                }
            }
        } else {
            for (int kk = 0; kk < k_extent; ++kk) {
                for (int i = kk + 1; i < k_extent; ++i) {
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, k0 + i) -=
                            Lkk[i * tile_size + kk] * TDLS_LUPP_XW(w, k0 + kk);
                }
            }
        }
 
        // Push into the tiles below
        for (int m = k + 1; m < full_tiles; ++m)
            fwd_push_one<k_extent, tile_size, pass_width, internal_rhs, internal_piv,
                         internal_matrix>(A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride,
                                          k0, m * tile_size);
        if constexpr (last_tile_tail > 0 && k_extent == tile_size)
            fwd_push_one<k_extent, last_tile_tail, pass_width, internal_rhs, internal_piv,
                         internal_matrix>(A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride,
                                          k0, full_tiles * tile_size);
    }
 
    template<int k_extent, int m_extent, int pass_width, bool internal_rhs, bool internal_piv,
             bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    bwd_pull_one(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                 const int piv_stride, T* TDLS_RESTRICT x, const int rhs_stride,
                 const std::size_t xcol_stride, const int k0, const int m0) noexcept {
        T Ukm[tile_size * tile_size];
        load_tile_piv<k_extent, m_extent, internal_piv, internal_matrix>(A, A_stride, piv,
                                                                         piv_stride, k0, m0, Ukm);
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = 0; i < k_extent; ++i) {
                TDLS_UNROLL_FORCE
                for (int w = 0; w < pass_width; ++w) {
                    T sum = T(0);
                    TDLS_UNROLL_FORCE
                    for (int j = 0; j < m_extent; ++j)
                        sum += Ukm[i * tile_size + j] * TDLS_LUPP_XW(w, m0 + j);
                    TDLS_LUPP_XW(w, k0 + i) -= sum;
                }
            }
        } else {
            for (int i = 0; i < k_extent; ++i) {
                for (int w = 0; w < pass_width; ++w) {
                    T sum = T(0);
                    for (int j = 0; j < m_extent; ++j)
                        sum += Ukm[i * tile_size + j] * TDLS_LUPP_XW(w, m0 + j);
                    TDLS_LUPP_XW(w, k0 + i) -= sum;
                }
            }
        }
    }
 
    template<int k_extent, int pass_width, bool internal_rhs, bool internal_piv,
             bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    bwd_step(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
             const int piv_stride, T* TDLS_RESTRICT x, const int rhs_stride,
             const std::size_t xcol_stride, const int k) noexcept {
        const int k0 = k * tile_size;
 
        // Pull the trailing contributions
        for (int m = k + 1; m < full_tiles; ++m)
            bwd_pull_one<k_extent, tile_size, pass_width, internal_rhs, internal_piv,
                         internal_matrix>(A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride,
                                          k0, m * tile_size);
        if constexpr (last_tile_tail > 0 && k_extent == tile_size)
            bwd_pull_one<k_extent, last_tile_tail, pass_width, internal_rhs, internal_piv,
                         internal_matrix>(A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride,
                                          k0, full_tiles * tile_size);
 
        // In-tile upper solve
        T Ukk[tile_size * tile_size];
        load_tile_piv_upper<k_extent, internal_piv, internal_matrix>(A, A_stride, piv, piv_stride,
                                                                     k0, k0, Ukk);
 
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int kk = k_extent - 1; kk >= 0; --kk) {
                TDLS_UNROLL_FORCE
                for (int w = 0; w < pass_width; ++w)
                    TDLS_LUPP_XW(w, k0 + kk) *= Ukk[kk * tile_size + kk]; // diag holds 1/pivot
                TDLS_UNROLL_FORCE
                for (int i = 0; i < kk; ++i) {
                    TDLS_UNROLL_FORCE
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, k0 + i) -=
                            Ukk[i * tile_size + kk] * TDLS_LUPP_XW(w, k0 + kk);
                }
            }
        } else {
            for (int kk = k_extent - 1; kk >= 0; --kk) {
                for (int w = 0; w < pass_width; ++w)
                    TDLS_LUPP_XW(w, k0 + kk) *= Ukk[kk * tile_size + kk]; // diag holds 1/pivot
                for (int i = 0; i < kk; ++i) {
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, k0 + i) -=
                            Ukk[i * tile_size + kk] * TDLS_LUPP_XW(w, k0 + kk);
                }
            }
        }
    }
 
    template<int pass_width, bool internal_rhs, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    bwd_only(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
             const int piv_stride, T* TDLS_RESTRICT x, const int rhs_stride,
             const std::size_t xcol_stride) noexcept {
        if constexpr (last_tile_tail > 0)
            bwd_step<last_tile_tail, pass_width, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride, full_tiles);
        for (int k = full_tiles - 1; k >= 0; --k)
            bwd_step<tile_size, pass_width, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride, k);
    }
 
    template<int pass_width, bool internal_rhs, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    fwd_bwd(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
            const int piv_stride, T* TDLS_RESTRICT x, const int rhs_stride,
            const std::size_t xcol_stride) noexcept {
        for (int k = 0; k < full_tiles; ++k)
            fwd_step<tile_size, pass_width, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride, k);
        if constexpr (last_tile_tail > 0)
            fwd_step<last_tile_tail, pass_width, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride, full_tiles);
 
        bwd_only<pass_width, internal_rhs, internal_piv, internal_matrix>(
            A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
    }
 
    /* =====================================================================
       TRANSPOSED SUBSTITUTION - column-major layout. With P S = L U and
       S = A^T, A = U^T L^T P, so A x = b is U^T z = b (forward, the U
       tiles read transposed, the diagonal holding 1/pivot), then
       L^T w = z (backward, the L tiles read transposed, unit diagonal),
       then x = P^T w. The right-hand side is indexed by the columns of S
       and is never permuted before the sweeps. The tile kernels below are
       shared by the fused forward pass of solve_inplace and by the
       separate passes, which therefore agree bitwise when the compiler
       contracts multiply-adds alike in both (-ffp-contract=off, or
       --fmad=false for nvcc).
       ===================================================================== */
 
    template<int k_extent, int pass_width, bool internal_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_solve_tile_upper(const T* TDLS_RESTRICT t, T* TDLS_RESTRICT x, const int rhs_stride,
                                const std::size_t xcol_stride, const int k0) noexcept {
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int kk = 0; kk < k_extent; ++kk) {
                TDLS_UNROLL_FORCE
                for (int w = 0; w < pass_width; ++w)
                    TDLS_LUPP_XW(w, k0 + kk) *= t[kk * tile_size + kk]; // diag holds 1/pivot
                TDLS_UNROLL_FORCE
                for (int i = kk + 1; i < k_extent; ++i) {
                    TDLS_UNROLL_FORCE
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, k0 + i) -= t[kk * tile_size + i] * TDLS_LUPP_XW(w, k0 + kk);
                }
            }
        } else {
            for (int kk = 0; kk < k_extent; ++kk) {
                for (int w = 0; w < pass_width; ++w)
                    TDLS_LUPP_XW(w, k0 + kk) *= t[kk * tile_size + kk]; // diag holds 1/pivot
                for (int i = kk + 1; i < k_extent; ++i) {
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, k0 + i) -= t[kk * tile_size + i] * TDLS_LUPP_XW(w, k0 + kk);
                }
            }
        }
    }
 
    template<int k_extent, int m_extent, int pass_width, bool internal_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_push_one(const T* TDLS_RESTRICT u, T* TDLS_RESTRICT x, const int rhs_stride,
                        const std::size_t xcol_stride, const int k0, const int m0) noexcept {
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = 0; i < m_extent; ++i) {
                TDLS_UNROLL_FORCE
                for (int w = 0; w < pass_width; ++w) {
                    T sum = T(0);
                    TDLS_UNROLL_FORCE
                    for (int j = 0; j < k_extent; ++j)
                        sum += u[j * tile_size + i] * TDLS_LUPP_XW(w, k0 + j);
                    TDLS_LUPP_XW(w, m0 + i) -= sum;
                }
            }
        } else {
            for (int i = 0; i < m_extent; ++i) {
                for (int w = 0; w < pass_width; ++w) {
                    T sum = T(0);
                    for (int j = 0; j < k_extent; ++j)
                        sum += u[j * tile_size + i] * TDLS_LUPP_XW(w, k0 + j);
                    TDLS_LUPP_XW(w, m0 + i) -= sum;
                }
            }
        }
    }
 
    template<int k_extent, int m_extent, int pass_width, bool internal_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_pull_one(const T* TDLS_RESTRICT l, T* TDLS_RESTRICT x, const int rhs_stride,
                        const std::size_t xcol_stride, const int k0, const int m0) noexcept {
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = 0; i < k_extent; ++i) {
                TDLS_UNROLL_FORCE
                for (int w = 0; w < pass_width; ++w) {
                    T sum = T(0);
                    TDLS_UNROLL_FORCE
                    for (int j = 0; j < m_extent; ++j)
                        sum += l[j * tile_size + i] * TDLS_LUPP_XW(w, m0 + j);
                    TDLS_LUPP_XW(w, k0 + i) -= sum;
                }
            }
        } else {
            for (int i = 0; i < k_extent; ++i) {
                for (int w = 0; w < pass_width; ++w) {
                    T sum = T(0);
                    for (int j = 0; j < m_extent; ++j)
                        sum += l[j * tile_size + i] * TDLS_LUPP_XW(w, m0 + j);
                    TDLS_LUPP_XW(w, k0 + i) -= sum;
                }
            }
        }
    }
 
    template<int k_extent, int pass_width, bool internal_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_solve_tile_lower_unit(const T* TDLS_RESTRICT l, T* TDLS_RESTRICT x,
                                     const int rhs_stride, const std::size_t xcol_stride,
                                     const int k0) noexcept {
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int kk = k_extent - 1; kk > 0; --kk) {
                TDLS_UNROLL_FORCE
                for (int i = 0; i < kk; ++i) {
                    TDLS_UNROLL_FORCE
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, k0 + i) -= l[kk * tile_size + i] * TDLS_LUPP_XW(w, k0 + kk);
                }
            }
        } else {
            for (int kk = k_extent - 1; kk > 0; --kk) {
                for (int i = 0; i < kk; ++i) {
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, k0 + i) -= l[kk * tile_size + i] * TDLS_LUPP_XW(w, k0 + kk);
                }
            }
        }
    }
 
    template<int k_extent, int pass_width, bool internal_rhs, bool internal_piv,
             bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_fwd_step(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                        const int piv_stride, T* TDLS_RESTRICT x, const int rhs_stride,
                        const std::size_t xcol_stride, const int k) noexcept {
        const int k0 = k * tile_size;
        T Ukk[tile_size * tile_size];
        load_tile_piv_upper<k_extent, internal_piv, internal_matrix>(A, A_stride, piv, piv_stride,
                                                                     k0, k0, Ukk);
        transposed_solve_tile_upper<k_extent, pass_width, internal_rhs>(Ukk, x, rhs_stride,
                                                                        xcol_stride, k0);
        for (int m = k + 1; m < full_tiles; ++m) {
            T Ukm[tile_size * tile_size];
            load_tile_piv<k_extent, tile_size, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, k0, m * tile_size, Ukm);
            transposed_push_one<k_extent, tile_size, pass_width, internal_rhs>(
                Ukm, x, rhs_stride, xcol_stride, k0, m * tile_size);
        }
        if constexpr (last_tile_tail > 0 && k_extent == tile_size) {
            T Ukm[tile_size * tile_size];
            load_tile_piv<k_extent, last_tile_tail, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, k0, full_tiles * tile_size, Ukm);
            transposed_push_one<k_extent, last_tile_tail, pass_width, internal_rhs>(
                Ukm, x, rhs_stride, xcol_stride, k0, full_tiles * tile_size);
        }
    }
 
    template<int k_extent, int pass_width, bool internal_rhs, bool internal_piv,
             bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_bwd_step(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                        const int piv_stride, T* TDLS_RESTRICT x, const int rhs_stride,
                        const std::size_t xcol_stride, const int k) noexcept {
        const int k0 = k * tile_size;
        for (int m = k + 1; m < full_tiles; ++m) {
            T Lmk[tile_size * tile_size];
            load_tile_piv<tile_size, k_extent, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, m * tile_size, k0, Lmk);
            transposed_pull_one<k_extent, tile_size, pass_width, internal_rhs>(
                Lmk, x, rhs_stride, xcol_stride, k0, m * tile_size);
        }
        if constexpr (last_tile_tail > 0 && k_extent == tile_size) {
            T Lmk[tile_size * tile_size];
            load_tile_piv<last_tile_tail, k_extent, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, full_tiles * tile_size, k0, Lmk);
            transposed_pull_one<k_extent, last_tile_tail, pass_width, internal_rhs>(
                Lmk, x, rhs_stride, xcol_stride, k0, full_tiles * tile_size);
        }
        T Lkk[tile_size * tile_size];
        load_tile_piv_lower<k_extent, internal_piv, internal_matrix>(A, A_stride, piv, piv_stride,
                                                                     k0, k0, Lkk);
        transposed_solve_tile_lower_unit<k_extent, pass_width, internal_rhs>(Lkk, x, rhs_stride,
                                                                             xcol_stride, k0);
    }
 
    template<int pass_width, bool internal_rhs, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_bwd_only(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                        const int piv_stride, T* TDLS_RESTRICT x, const int rhs_stride,
                        const std::size_t xcol_stride) noexcept {
        if constexpr (last_tile_tail > 0)
            transposed_bwd_step<last_tile_tail, pass_width, internal_rhs, internal_piv,
                                internal_matrix>(A, A_stride, piv, piv_stride, x, rhs_stride,
                                                 xcol_stride, full_tiles);
        for (int k = full_tiles - 1; k >= 0; --k)
            transposed_bwd_step<tile_size, pass_width, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride, k);
    }
 
    template<int pass_width, bool internal_rhs, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_fwd_bwd(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                       const int piv_stride, T* TDLS_RESTRICT x, const int rhs_stride,
                       const std::size_t xcol_stride) noexcept {
        for (int k = 0; k < full_tiles; ++k)
            transposed_fwd_step<tile_size, pass_width, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride, k);
        if constexpr (last_tile_tail > 0)
            transposed_fwd_step<last_tile_tail, pass_width, internal_rhs, internal_piv,
                                internal_matrix>(A, A_stride, piv, piv_stride, x, rhs_stride,
                                                 xcol_stride, full_tiles);
 
        transposed_bwd_only<pass_width, internal_rhs, internal_piv, internal_matrix>(
            A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
    }
 
    template<int pass_width, bool internal_rhs, bool internal_piv>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    unpermute_cycles(const int* TDLS_RESTRICT piv, const int piv_stride, T* TDLS_RESTRICT x,
                     const int rhs_stride, const std::size_t xcol_stride) noexcept {
        if constexpr (N <= 64) {
            using mask_t   = std::conditional_t<(N <= 32), unsigned, unsigned long long>;
            mask_t visited = mask_t(0);
            if constexpr (Config.unroll_loops) {
                TDLS_UNROLL_FORCE
                for (int s = 0; s < N; ++s) {
                    if ((visited >> s) & mask_t(1)) continue;
                    T carry[pass_width];
                    TDLS_UNROLL_FORCE
                    for (int w = 0; w < pass_width; ++w)
                        carry[w] = TDLS_LUPP_XW(w, s);
                    int dst = TDLS_LUPP_PIV(s);
                    while (dst != s) {
                        TDLS_UNROLL_FORCE
                        for (int w = 0; w < pass_width; ++w) {
                            const T displaced    = TDLS_LUPP_XW(w, dst);
                            TDLS_LUPP_XW(w, dst) = carry[w];
                            carry[w]             = displaced;
                        }
                        visited |= mask_t(1) << dst;
                        dst = TDLS_LUPP_PIV(dst);
                    }
                    TDLS_UNROLL_FORCE
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, s) = carry[w];
                    visited |= mask_t(1) << s;
                }
            } else {
                for (int s = 0; s < N; ++s) {
                    if ((visited >> s) & mask_t(1)) continue;
                    T carry[pass_width];
                    for (int w = 0; w < pass_width; ++w)
                        carry[w] = TDLS_LUPP_XW(w, s);
                    int dst = TDLS_LUPP_PIV(s);
                    while (dst != s) {
                        for (int w = 0; w < pass_width; ++w) {
                            const T displaced    = TDLS_LUPP_XW(w, dst);
                            TDLS_LUPP_XW(w, dst) = carry[w];
                            carry[w]             = displaced;
                        }
                        visited |= mask_t(1) << dst;
                        dst = TDLS_LUPP_PIV(dst);
                    }
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, s) = carry[w];
                    visited |= mask_t(1) << s;
                }
            }
        } else {
            for (int s = 0; s < N; ++s) {
                int probe = TDLS_LUPP_PIV(s);
                while (probe > s)
                    probe = TDLS_LUPP_PIV(probe);
                if (probe != s) continue;
 
                T carry[pass_width];
                for (int w = 0; w < pass_width; ++w)
                    carry[w] = TDLS_LUPP_XW(w, s);
                int dst = TDLS_LUPP_PIV(s);
                while (dst != s) {
                    for (int w = 0; w < pass_width; ++w) {
                        const T displaced    = TDLS_LUPP_XW(w, dst);
                        TDLS_LUPP_XW(w, dst) = carry[w];
                        carry[w]             = displaced;
                    }
                    dst = TDLS_LUPP_PIV(dst);
                }
                for (int w = 0; w < pass_width; ++w)
                    TDLS_LUPP_XW(w, s) = carry[w];
            }
        }
    }
 
    template<int pass_width, bool internal_piv>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    unpermute_predicated(const int* TDLS_RESTRICT piv, const int piv_stride,
                         T* TDLS_RESTRICT x) noexcept {
        constexpr bool internal_rhs       = true;
        constexpr int rhs_stride          = 1;
        constexpr std::size_t xcol_stride = N;
        T w_copy[pass_width * N];
        int p[N];
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = 0; i < pass_width * N; ++i)
                w_copy[i] = x[i];
            TDLS_UNROLL_FORCE
            for (int i = 0; i < N; ++i)
                p[i] = TDLS_LUPP_PIV(i);
            TDLS_UNROLL_FORCE
            for (int j = 0; j < N; ++j) {
                TDLS_UNROLL_FORCE
                for (int w = 0; w < pass_width; ++w) {
                    T v = T(0);
                    TDLS_UNROLL_FORCE
                    for (int i = 0; i < N; ++i)
                        if (p[i] == j) v = w_copy[w * N + i];
                    TDLS_LUPP_XW(w, j) = v;
                }
            }
        } else {
            for (int i = 0; i < pass_width * N; ++i)
                w_copy[i] = x[i];
            for (int i = 0; i < N; ++i)
                p[i] = TDLS_LUPP_PIV(i);
            for (int j = 0; j < N; ++j) {
                for (int w = 0; w < pass_width; ++w) {
                    T v = T(0);
                    for (int i = 0; i < N; ++i)
                        if (p[i] == j) v = w_copy[w * N + i];
                    TDLS_LUPP_XW(w, j) = v;
                }
            }
        }
    }
 
    template<int pass_width, bool internal_rhs, bool internal_piv>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    unpermute(const int* TDLS_RESTRICT piv, const int piv_stride, T* TDLS_RESTRICT x,
              const int rhs_stride, const std::size_t xcol_stride) noexcept {
        if constexpr (internal_rhs)
            unpermute_predicated<pass_width, internal_piv>(piv, piv_stride, x);
        else
            unpermute_cycles<pass_width, internal_rhs, internal_piv>(piv, piv_stride, x, rhs_stride,
                                                                     xcol_stride);
    }
 
    /* =====================================================================
       SUBSTITUTION - public entry points.
       ===================================================================== */
 
    template<bool internal_rhs, bool internal_piv, bool internal_matrix,
             int A_static_stride = dynamic_stride, int piv_static_stride = dynamic_stride,
             int rhs_static_stride = dynamic_stride>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
               const int piv_stride, const T* TDLS_RESTRICT b, T* TDLS_RESTRICT x,
               const int rhs_stride) noexcept {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return substitute<internal_rhs, internal_piv, internal_matrix>(
                A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), b, x,
                detail::stride_or<rhs_static_stride>(rhs_stride));
        } else {
            if constexpr (physical) {
                TDLS_LUPP_FOR(i, 0, N, TDLS_LUPP_X(i) = TDLS_LUPP_B(i);)
                physical_multirhs<1, 0, 3, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x, rhs_stride, 0);
                return;
            }
            if constexpr (transposed) {
                // Column-major layout: b is consumed unpermuted, the solution is
                // permuted at the end (x = P^T w).
                if constexpr (Config.unroll_loops) {
                    TDLS_UNROLL_FORCE
                    for (int i = 0; i < N; ++i)
                        TDLS_LUPP_X(i) = TDLS_LUPP_B(i);
                } else {
                    for (int i = 0; i < N; ++i)
                        TDLS_LUPP_X(i) = TDLS_LUPP_B(i);
                }
                transposed_fwd_bwd<1, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x, rhs_stride, 0);
                unpermute<1, internal_rhs, internal_piv>(piv, piv_stride, x, rhs_stride, 0);
            } else {
                if constexpr (internal_rhs) {
                    // Predicated gather: b[piv[i]] would index the internal b with
                    // a runtime value and demote it to local memory. The equality
                    // sweep keeps every index compile-time (measured: local memory
                    // eliminated, identical values).
                    if constexpr (Config.unroll_loops) {
                        TDLS_UNROLL_FORCE
                        for (int i = 0; i < N; ++i) {
                            const int p = TDLS_LUPP_PIV(i);
                            T v         = T(0);
                            TDLS_UNROLL_FORCE
                            for (int j = 0; j < N; ++j)
                                if (p == j) v = b[j];
                            TDLS_LUPP_X(i) = v;
                        }
                    } else {
                        for (int i = 0; i < N; ++i) {
                            const int p = TDLS_LUPP_PIV(i);
                            T v         = T(0);
                            for (int j = 0; j < N; ++j)
                                if (p == j) v = b[j];
                            TDLS_LUPP_X(i) = v;
                        }
                    }
                } else {
                    if constexpr (Config.unroll_loops) {
                        TDLS_UNROLL_FORCE
                        for (int i = 0; i < N; ++i)
                            TDLS_LUPP_X(i) = TDLS_LUPP_B(TDLS_LUPP_PIV(i));
                    } else {
                        for (int i = 0; i < N; ++i)
                            TDLS_LUPP_X(i) = TDLS_LUPP_B(TDLS_LUPP_PIV(i));
                    }
                }
                fwd_bwd<1, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x, rhs_stride, 0);
            }
        }
    }
 
    template<bool internal_rhs, bool internal_piv, bool internal_matrix,
             int A_static_stride = dynamic_stride, int piv_static_stride = dynamic_stride,
             int rhs_static_stride = dynamic_stride>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_canonical(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                         const int piv_stride, const int col, T* TDLS_RESTRICT x,
                         const int rhs_stride) noexcept {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return substitute_canonical<internal_rhs, internal_piv, internal_matrix>(
                A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), col, x,
                detail::stride_or<rhs_static_stride>(rhs_stride));
        } else {
            if constexpr (physical) {
                TDLS_LUPP_FOR(i, 0, N, TDLS_LUPP_X(i) = i == col ? T(1) : T(0);)
                physical_multirhs<1, 0, 3, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x, rhs_stride, 0);
                return;
            }
            substitute_canonical_multirhs_pass<1, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, col, x, rhs_stride, 0);
        }
    }
 
    template<int pass_width, bool internal_rhs, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_canonical_multirhs_pass(const T* TDLS_RESTRICT A, const int A_stride,
                                       const int* TDLS_RESTRICT piv, const int piv_stride,
                                       const int col0, T* TDLS_RESTRICT x, const int rhs_stride,
                                       const std::size_t xcol_stride) noexcept {
        if constexpr (transposed) {
            // Column-major layout: e_col unpermuted, the solution permuted at
            // the end (x = P^T w).
            if constexpr (Config.unroll_loops) {
                TDLS_UNROLL_FORCE
                for (int i = 0; i < N; ++i) {
                    TDLS_UNROLL_FORCE
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, i) = (i == col0 + w) ? T(1) : T(0);
                }
            } else {
                for (int i = 0; i < N; ++i) {
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, i) = (i == col0 + w) ? T(1) : T(0);
                }
            }
            transposed_fwd_bwd<pass_width, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
            unpermute<pass_width, internal_rhs, internal_piv>(piv, piv_stride, x, rhs_stride,
                                                              xcol_stride);
        } else {
            if constexpr (Config.unroll_loops) {
                TDLS_UNROLL_FORCE
                for (int i = 0; i < N; ++i) {
                    const int p = TDLS_LUPP_PIV(i);
                    TDLS_UNROLL_FORCE
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, i) = (p == col0 + w) ? T(1) : T(0);
                }
            } else {
                for (int i = 0; i < N; ++i) {
                    const int p = TDLS_LUPP_PIV(i);
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, i) = (p == col0 + w) ? T(1) : T(0);
                }
            }
            fwd_bwd<pass_width, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
        }
    }
 
    template<int nrhs, bool internal_rhs, bool internal_piv, bool internal_matrix,
             int pass_width = 0, int A_static_stride = dynamic_stride,
             int piv_static_stride = dynamic_stride, int rhs_static_stride = dynamic_stride>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_canonical_multirhs(const T* TDLS_RESTRICT A, const int A_stride,
                                  const int* TDLS_RESTRICT piv, const int piv_stride,
                                  const int col0, T* TDLS_RESTRICT x, const int rhs_stride,
                                  const std::size_t xcol_stride) noexcept {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return substitute_canonical_multirhs<nrhs, internal_rhs, internal_piv, internal_matrix,
                                                 pass_width>(
                A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), col0, x,
                detail::stride_or<rhs_static_stride>(rhs_stride), xcol_stride);
        } else {
            if constexpr (physical) {
                TDLS_LUPP_FOR(
                    c, 0, nrhs,
                    TDLS_LUPP_FOR(i, 0, N, TDLS_LUPP_XW(c, i) = i == col0 + c ? T(1) : T(0);))
                physical_multirhs<nrhs, pass_width, 3, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
                return;
            }
            static_assert(nrhs >= 1, "tdls: nrhs must be at least 1");
            static_assert(pass_width >= 0, "tdls: pass_width must not be negative");
            constexpr int columns_per_pass =
                (pass_width <= 0 || pass_width >= nrhs) ? nrhs : pass_width;
            for (int c0 = 0; c0 + columns_per_pass <= nrhs; c0 += columns_per_pass) {
                const std::size_t off =
                    internal_rhs ? unsigned(c0) * unsigned(N) : TDLS_STRIDED_WIDE(c0, xcol_stride);
                substitute_canonical_multirhs_pass<columns_per_pass, internal_rhs, internal_piv,
                                                   internal_matrix>(
                    A, A_stride, piv, piv_stride, col0 + c0, x + off, rhs_stride, xcol_stride);
            }
            if constexpr (nrhs % columns_per_pass > 0) {
                constexpr int c0 = nrhs - nrhs % columns_per_pass;
                const std::size_t off =
                    internal_rhs ? unsigned(c0) * unsigned(N) : TDLS_STRIDED_WIDE(c0, xcol_stride);
                substitute_canonical_multirhs_pass<nrhs % columns_per_pass, internal_rhs,
                                                   internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, col0 + c0, x + off, rhs_stride, xcol_stride);
            }
        }
    }
 
    template<int pass_width, bool internal_rhs, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_multirhs_pass(const T* TDLS_RESTRICT A, const int A_stride,
                             const int* TDLS_RESTRICT piv, const int piv_stride,
                             const T* TDLS_RESTRICT b, T* TDLS_RESTRICT x, const int rhs_stride,
                             const std::size_t xcol_stride) noexcept {
        if constexpr (transposed) {
            // Column-major layout, as in substitute.
            if constexpr (Config.unroll_loops) {
                TDLS_UNROLL_FORCE
                for (int i = 0; i < N; ++i) {
                    TDLS_UNROLL_FORCE
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, i) = TDLS_LUPP_BW(w, i);
                }
            } else {
                for (int i = 0; i < N; ++i) {
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, i) = TDLS_LUPP_BW(w, i);
                }
            }
            transposed_fwd_bwd<pass_width, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
            unpermute<pass_width, internal_rhs, internal_piv>(piv, piv_stride, x, rhs_stride,
                                                              xcol_stride);
        } else {
            if constexpr (internal_rhs) {
                // Predicated gather, exactly as in substitute: b[piv[i]] would
                // index the internal block with a runtime value and demote it
                // to local memory.
                if constexpr (Config.unroll_loops) {
                    TDLS_UNROLL_FORCE
                    for (int i = 0; i < N; ++i) {
                        const int p = TDLS_LUPP_PIV(i);
                        TDLS_UNROLL_FORCE
                        for (int w = 0; w < pass_width; ++w) {
                            T v = T(0);
                            TDLS_UNROLL_FORCE
                            for (int j = 0; j < N; ++j)
                                if (p == j) v = TDLS_LUPP_BW(w, j);
                            TDLS_LUPP_XW(w, i) = v;
                        }
                    }
                } else {
                    for (int i = 0; i < N; ++i) {
                        const int p = TDLS_LUPP_PIV(i);
                        for (int w = 0; w < pass_width; ++w) {
                            T v = T(0);
                            for (int j = 0; j < N; ++j)
                                if (p == j) v = TDLS_LUPP_BW(w, j);
                            TDLS_LUPP_XW(w, i) = v;
                        }
                    }
                }
            } else {
                if constexpr (Config.unroll_loops) {
                    TDLS_UNROLL_FORCE
                    for (int i = 0; i < N; ++i) {
                        const int p = TDLS_LUPP_PIV(i);
                        TDLS_UNROLL_FORCE
                        for (int w = 0; w < pass_width; ++w)
                            TDLS_LUPP_XW(w, i) = TDLS_LUPP_BW(w, p);
                    }
                } else {
                    for (int i = 0; i < N; ++i) {
                        const int p = TDLS_LUPP_PIV(i);
                        for (int w = 0; w < pass_width; ++w)
                            TDLS_LUPP_XW(w, i) = TDLS_LUPP_BW(w, p);
                    }
                }
            }
            fwd_bwd<pass_width, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
        }
    }
 
    template<int nrhs, bool internal_rhs, bool internal_piv, bool internal_matrix,
             int pass_width = 0, int A_static_stride = dynamic_stride,
             int piv_static_stride = dynamic_stride, int rhs_static_stride = dynamic_stride>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_multirhs(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                        const int piv_stride, const T* TDLS_RESTRICT b, T* TDLS_RESTRICT x,
                        const int rhs_stride, const std::size_t xcol_stride) noexcept {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return substitute_multirhs<nrhs, internal_rhs, internal_piv, internal_matrix,
                                       pass_width>(
                A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), b, x,
                detail::stride_or<rhs_static_stride>(rhs_stride), xcol_stride);
        } else {
            if constexpr (physical) {
                TDLS_LUPP_FOR(c, 0, nrhs,
                              TDLS_LUPP_FOR(i, 0, N, TDLS_LUPP_XW(c, i) = TDLS_LUPP_BW(c, i);))
                physical_multirhs<nrhs, pass_width, 3, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
                return;
            }
            static_assert(nrhs >= 1, "tdls: nrhs must be at least 1");
            static_assert(pass_width >= 0, "tdls: pass_width must not be negative");
            constexpr int columns_per_pass =
                (pass_width <= 0 || pass_width >= nrhs) ? nrhs : pass_width;
            for (int c0 = 0; c0 + columns_per_pass <= nrhs; c0 += columns_per_pass) {
                const std::size_t off =
                    internal_rhs ? unsigned(c0) * unsigned(N) : TDLS_STRIDED_WIDE(c0, xcol_stride);
                substitute_multirhs_pass<columns_per_pass, internal_rhs, internal_piv,
                                         internal_matrix>(A, A_stride, piv, piv_stride, b + off,
                                                          x + off, rhs_stride, xcol_stride);
            }
            if constexpr (nrhs % columns_per_pass > 0) {
                constexpr int c0 = nrhs - nrhs % columns_per_pass;
                const std::size_t off =
                    internal_rhs ? unsigned(c0) * unsigned(N) : TDLS_STRIDED_WIDE(c0, xcol_stride);
                substitute_multirhs_pass<nrhs % columns_per_pass, internal_rhs, internal_piv,
                                         internal_matrix>(A, A_stride, piv, piv_stride, b + off,
                                                          x + off, rhs_stride, xcol_stride);
            }
        }
    }
 
    template<bool internal_rhs, bool internal_piv, bool internal_matrix,
             int A_static_stride = dynamic_stride, int piv_static_stride = dynamic_stride,
             int rhs_static_stride = dynamic_stride>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_inplace(const T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT piv,
                       const int piv_stride, T* TDLS_RESTRICT x, const int rhs_stride) noexcept {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return substitute_inplace<internal_rhs, internal_piv, internal_matrix>(
                A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), x,
                detail::stride_or<rhs_static_stride>(rhs_stride));
        } else {
            if constexpr (physical) {
                physical_multirhs<1, 0, 3, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x, rhs_stride, 0);
                return;
            }
            if constexpr (transposed) {
                // Column-major layout: the sweeps on the unpermuted x, then the
                // inverse cycle rotation (x = P^T w).
                transposed_fwd_bwd<1, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x, rhs_stride, 0);
                unpermute_cycles<1, internal_rhs, internal_piv>(piv, piv_stride, x, rhs_stride, 0);
            } else {
                if constexpr (N <= 64) {
                    // Bitmask cycle decomposition: one visited bit per row, the
                    // whole state in a single 32- or 64-bit register.
                    using mask_t   = std::conditional_t<(N <= 32), unsigned, unsigned long long>;
                    mask_t visited = mask_t(0);
                    if constexpr (Config.unroll_loops) {
                        TDLS_UNROLL_FORCE
                        for (int s = 0; s < N; ++s) {
                            if ((visited >> s) & mask_t(1)) continue;
                            const T tmp = TDLS_LUPP_X(s);
                            int cur     = s;
                            int nxt     = TDLS_LUPP_PIV(cur);
                            while (nxt != s) {
                                TDLS_LUPP_X(cur) = TDLS_LUPP_X(nxt);
                                visited |= mask_t(1) << cur;
                                cur = nxt;
                                nxt = TDLS_LUPP_PIV(cur);
                            }
                            TDLS_LUPP_X(cur) = tmp;
                            visited |= mask_t(1) << cur;
                        }
                    } else {
                        for (int s = 0; s < N; ++s) {
                            if ((visited >> s) & mask_t(1)) continue;
                            const T tmp = TDLS_LUPP_X(s);
                            int cur     = s;
                            int nxt     = TDLS_LUPP_PIV(cur);
                            while (nxt != s) {
                                TDLS_LUPP_X(cur) = TDLS_LUPP_X(nxt);
                                visited |= mask_t(1) << cur;
                                cur = nxt;
                                nxt = TDLS_LUPP_PIV(cur);
                            }
                            TDLS_LUPP_X(cur) = tmp;
                            visited |= mask_t(1) << cur;
                        }
                    }
                } else {
                    // Cycle-leader scan (no visited storage): a cycle is rotated
                    // only when reached from its smallest index, detected by
                    // walking the orbit. Values and order match the mask path.
                    for (int s = 0; s < N; ++s) {
                        int probe = TDLS_LUPP_PIV(s);
                        while (probe > s)
                            probe = TDLS_LUPP_PIV(probe);
                        if (probe != s) continue;
 
                        const T tmp = TDLS_LUPP_X(s);
                        int cur     = s;
                        int nxt     = TDLS_LUPP_PIV(cur);
                        while (nxt != s) {
                            TDLS_LUPP_X(cur) = TDLS_LUPP_X(nxt);
                            cur              = nxt;
                            nxt              = TDLS_LUPP_PIV(cur);
                        }
                        TDLS_LUPP_X(cur) = tmp;
                    }
                }
                fwd_bwd<1, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x, rhs_stride, 0);
            }
        }
    }
 
    template<int pass_width, bool internal_rhs, bool internal_piv>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL static constexpr void
    permute_cycles(const int* TDLS_RESTRICT piv, const int piv_stride, T* TDLS_RESTRICT x,
                   const int rhs_stride, const std::size_t xcol_stride) noexcept {
        if constexpr (N <= 64) {
            // Bitmask cycle decomposition, as in substitute_inplace; tmp
            // widens to one entry per column.
            using mask_t   = std::conditional_t<(N <= 32), unsigned, unsigned long long>;
            mask_t visited = mask_t(0);
            if constexpr (Config.unroll_loops) {
                TDLS_UNROLL_FORCE
                for (int s = 0; s < N; ++s) {
                    if ((visited >> s) & mask_t(1)) continue;
                    T tmp[pass_width];
                    TDLS_UNROLL_FORCE
                    for (int w = 0; w < pass_width; ++w)
                        tmp[w] = TDLS_LUPP_XW(w, s);
                    int cur = s;
                    int nxt = TDLS_LUPP_PIV(cur);
                    while (nxt != s) {
                        TDLS_UNROLL_FORCE
                        for (int w = 0; w < pass_width; ++w)
                            TDLS_LUPP_XW(w, cur) = TDLS_LUPP_XW(w, nxt);
                        visited |= mask_t(1) << cur;
                        cur = nxt;
                        nxt = TDLS_LUPP_PIV(cur);
                    }
                    TDLS_UNROLL_FORCE
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, cur) = tmp[w];
                    visited |= mask_t(1) << cur;
                }
            } else {
                for (int s = 0; s < N; ++s) {
                    if ((visited >> s) & mask_t(1)) continue;
                    T tmp[pass_width];
                    for (int w = 0; w < pass_width; ++w)
                        tmp[w] = TDLS_LUPP_XW(w, s);
                    int cur = s;
                    int nxt = TDLS_LUPP_PIV(cur);
                    while (nxt != s) {
                        for (int w = 0; w < pass_width; ++w)
                            TDLS_LUPP_XW(w, cur) = TDLS_LUPP_XW(w, nxt);
                        visited |= mask_t(1) << cur;
                        cur = nxt;
                        nxt = TDLS_LUPP_PIV(cur);
                    }
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, cur) = tmp[w];
                    visited |= mask_t(1) << cur;
                }
            }
        } else {
            // Cycle-leader scan, as in substitute_inplace.
            for (int s = 0; s < N; ++s) {
                int probe = TDLS_LUPP_PIV(s);
                while (probe > s)
                    probe = TDLS_LUPP_PIV(probe);
                if (probe != s) continue;
 
                T tmp[pass_width];
                for (int w = 0; w < pass_width; ++w)
                    tmp[w] = TDLS_LUPP_XW(w, s);
                int cur = s;
                int nxt = TDLS_LUPP_PIV(cur);
                while (nxt != s) {
                    for (int w = 0; w < pass_width; ++w)
                        TDLS_LUPP_XW(w, cur) = TDLS_LUPP_XW(w, nxt);
                    cur = nxt;
                    nxt = TDLS_LUPP_PIV(cur);
                }
                for (int w = 0; w < pass_width; ++w)
                    TDLS_LUPP_XW(w, cur) = tmp[w];
            }
        }
    }
 
    template<int pass_width, bool internal_rhs, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_inplace_multirhs_pass(const T* TDLS_RESTRICT A, const int A_stride,
                                     const int* TDLS_RESTRICT piv, const int piv_stride,
                                     T* TDLS_RESTRICT x, const int rhs_stride,
                                     const std::size_t xcol_stride) noexcept {
        if constexpr (transposed) {
            // Column-major layout, as in substitute_inplace.
            transposed_fwd_bwd<pass_width, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
            unpermute_cycles<pass_width, internal_rhs, internal_piv>(piv, piv_stride, x, rhs_stride,
                                                                     xcol_stride);
        } else {
            permute_cycles<pass_width, internal_rhs, internal_piv>(piv, piv_stride, x, rhs_stride,
                                                                   xcol_stride);
            fwd_bwd<pass_width, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
        }
    }
 
    template<int pass_width, bool internal_rhs, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    forward_inplace_multirhs_pass(const T* TDLS_RESTRICT A, const int A_stride,
                                  const int* TDLS_RESTRICT piv, const int piv_stride,
                                  T* TDLS_RESTRICT x, const int rhs_stride,
                                  const std::size_t xcol_stride) noexcept {
        if constexpr (transposed) {
            for (int k = 0; k < full_tiles; ++k)
                transposed_fwd_step<tile_size, pass_width, internal_rhs, internal_piv,
                                    internal_matrix>(A, A_stride, piv, piv_stride, x, rhs_stride,
                                                     xcol_stride, k);
            if constexpr (last_tile_tail > 0)
                transposed_fwd_step<last_tile_tail, pass_width, internal_rhs, internal_piv,
                                    internal_matrix>(A, A_stride, piv, piv_stride, x, rhs_stride,
                                                     xcol_stride, full_tiles);
        } else {
            permute_cycles<pass_width, internal_rhs, internal_piv>(piv, piv_stride, x, rhs_stride,
                                                                   xcol_stride);
            for (int k = 0; k < full_tiles; ++k)
                fwd_step<tile_size, pass_width, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride, k);
            if constexpr (last_tile_tail > 0)
                fwd_step<last_tile_tail, pass_width, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride, full_tiles);
        }
    }
 
    template<int pass_width, bool internal_rhs, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    backward_inplace_multirhs_pass(const T* TDLS_RESTRICT A, const int A_stride,
                                   const int* TDLS_RESTRICT piv, const int piv_stride,
                                   T* TDLS_RESTRICT x, const int rhs_stride,
                                   const std::size_t xcol_stride) noexcept {
        if constexpr (transposed) {
            transposed_bwd_only<pass_width, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
            unpermute_cycles<pass_width, internal_rhs, internal_piv>(piv, piv_stride, x, rhs_stride,
                                                                     xcol_stride);
        } else {
            bwd_only<pass_width, internal_rhs, internal_piv, internal_matrix>(
                A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
        }
    }
 
    template<int nrhs, bool internal_rhs, bool internal_piv, bool internal_matrix,
             int pass_width = 0, int A_static_stride = dynamic_stride,
             int piv_static_stride = dynamic_stride, int rhs_static_stride = dynamic_stride>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_inplace_multirhs(const T* TDLS_RESTRICT A, const int A_stride,
                                const int* TDLS_RESTRICT piv, const int piv_stride,
                                T* TDLS_RESTRICT x, const int rhs_stride,
                                const std::size_t xcol_stride) noexcept {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return substitute_inplace_multirhs<nrhs, internal_rhs, internal_piv, internal_matrix,
                                               pass_width>(
                A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), x,
                detail::stride_or<rhs_static_stride>(rhs_stride), xcol_stride);
        } else {
            if constexpr (physical) {
                physical_multirhs<nrhs, pass_width, 3, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
                return;
            }
            static_assert(nrhs >= 1, "tdls: nrhs must be at least 1");
            static_assert(pass_width >= 0, "tdls: pass_width must not be negative");
            constexpr int columns_per_pass =
                (pass_width <= 0 || pass_width >= nrhs) ? nrhs : pass_width;
            for (int c0 = 0; c0 + columns_per_pass <= nrhs; c0 += columns_per_pass) {
                const std::size_t off =
                    internal_rhs ? unsigned(c0) * unsigned(N) : TDLS_STRIDED_WIDE(c0, xcol_stride);
                substitute_inplace_multirhs_pass<columns_per_pass, internal_rhs, internal_piv,
                                                 internal_matrix>(A, A_stride, piv, piv_stride,
                                                                  x + off, rhs_stride, xcol_stride);
            }
            if constexpr (nrhs % columns_per_pass > 0) {
                constexpr int c0 = nrhs - nrhs % columns_per_pass;
                const std::size_t off =
                    internal_rhs ? unsigned(c0) * unsigned(N) : TDLS_STRIDED_WIDE(c0, xcol_stride);
                substitute_inplace_multirhs_pass<nrhs % columns_per_pass, internal_rhs,
                                                 internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x + off, rhs_stride, xcol_stride);
            }
        }
    }
 
    template<int nrhs, bool internal_rhs, bool internal_piv, bool internal_matrix,
             int pass_width = 0, int A_static_stride = dynamic_stride,
             int piv_static_stride = dynamic_stride, int rhs_static_stride = dynamic_stride>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    forward_inplace_multirhs(const T* TDLS_RESTRICT A, const int A_stride,
                             const int* TDLS_RESTRICT piv, const int piv_stride, T* TDLS_RESTRICT x,
                             const int rhs_stride, const std::size_t xcol_stride) noexcept {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return forward_inplace_multirhs<nrhs, internal_rhs, internal_piv, internal_matrix,
                                            pass_width>(
                A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), x,
                detail::stride_or<rhs_static_stride>(rhs_stride), xcol_stride);
        } else {
            if constexpr (physical) {
                physical_multirhs<nrhs, pass_width, 1, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
                return;
            }
            static_assert(nrhs >= 1, "tdls: nrhs must be at least 1");
            static_assert(pass_width >= 0, "tdls: pass_width must not be negative");
            constexpr int columns_per_pass =
                (pass_width <= 0 || pass_width >= nrhs) ? nrhs : pass_width;
            for (int c0 = 0; c0 + columns_per_pass <= nrhs; c0 += columns_per_pass) {
                const std::size_t off =
                    internal_rhs ? unsigned(c0) * unsigned(N) : TDLS_STRIDED_WIDE(c0, xcol_stride);
                forward_inplace_multirhs_pass<columns_per_pass, internal_rhs, internal_piv,
                                              internal_matrix>(A, A_stride, piv, piv_stride,
                                                               x + off, rhs_stride, xcol_stride);
            }
            if constexpr (nrhs % columns_per_pass > 0) {
                constexpr int c0 = nrhs - nrhs % columns_per_pass;
                const std::size_t off =
                    internal_rhs ? unsigned(c0) * unsigned(N) : TDLS_STRIDED_WIDE(c0, xcol_stride);
                forward_inplace_multirhs_pass<nrhs % columns_per_pass, internal_rhs, internal_piv,
                                              internal_matrix>(A, A_stride, piv, piv_stride,
                                                               x + off, rhs_stride, xcol_stride);
            }
        }
    }
 
    template<int nrhs, bool internal_rhs, bool internal_piv, bool internal_matrix,
             int pass_width = 0, int A_static_stride = dynamic_stride,
             int piv_static_stride = dynamic_stride, int rhs_static_stride = dynamic_stride>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    backward_inplace_multirhs(const T* TDLS_RESTRICT A, const int A_stride,
                              const int* TDLS_RESTRICT piv, const int piv_stride,
                              T* TDLS_RESTRICT x, const int rhs_stride,
                              const std::size_t xcol_stride) noexcept {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return backward_inplace_multirhs<nrhs, internal_rhs, internal_piv, internal_matrix,
                                             pass_width>(
                A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), x,
                detail::stride_or<rhs_static_stride>(rhs_stride), xcol_stride);
        } else {
            if constexpr (physical) {
                physical_multirhs<nrhs, pass_width, 2, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
                return;
            }
            static_assert(nrhs >= 1, "tdls: nrhs must be at least 1");
            static_assert(pass_width >= 0, "tdls: pass_width must not be negative");
            constexpr int columns_per_pass =
                (pass_width <= 0 || pass_width >= nrhs) ? nrhs : pass_width;
            for (int c0 = 0; c0 + columns_per_pass <= nrhs; c0 += columns_per_pass) {
                const std::size_t off =
                    internal_rhs ? unsigned(c0) * unsigned(N) : TDLS_STRIDED_WIDE(c0, xcol_stride);
                backward_inplace_multirhs_pass<columns_per_pass, internal_rhs, internal_piv,
                                               internal_matrix>(A, A_stride, piv, piv_stride,
                                                                x + off, rhs_stride, xcol_stride);
            }
            if constexpr (nrhs % columns_per_pass > 0) {
                constexpr int c0 = nrhs - nrhs % columns_per_pass;
                const std::size_t off =
                    internal_rhs ? unsigned(c0) * unsigned(N) : TDLS_STRIDED_WIDE(c0, xcol_stride);
                backward_inplace_multirhs_pass<nrhs % columns_per_pass, internal_rhs, internal_piv,
                                               internal_matrix>(A, A_stride, piv, piv_stride,
                                                                x + off, rhs_stride, xcol_stride);
            }
        }
    }
 
    /* =====================================================================
       SOLVE - convenience wrappers.
       ===================================================================== */
 
    template<bool internal_rhs, bool internal_piv, bool internal_matrix,
             bool oot_diagnostics = true, int A_static_stride = dynamic_stride,
             int piv_static_stride = dynamic_stride, int rhs_static_stride = dynamic_stride>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv, const int piv_stride,
          const T* TDLS_RESTRICT b, T* TDLS_RESTRICT x, const int rhs_stride,
          int& oot_count) noexcept {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return solve<internal_rhs, internal_piv, internal_matrix, oot_diagnostics>(
                A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), b, x,
                detail::stride_or<rhs_static_stride>(rhs_stride), oot_count);
        } else {
            if (!factorize<internal_piv, internal_matrix, oot_diagnostics>(A, A_stride, piv,
                                                                           piv_stride, oot_count))
                return false;
            substitute<internal_rhs, internal_piv, internal_matrix>(A, A_stride, piv, piv_stride, b,
                                                                    x, rhs_stride);
            return true;
        }
    }
 
    template<bool internal_rhs, bool internal_piv, bool internal_matrix,
             int A_static_stride = dynamic_stride, int piv_static_stride = dynamic_stride,
             int rhs_static_stride = dynamic_stride>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv, const int piv_stride,
          const T* TDLS_RESTRICT b, T* TDLS_RESTRICT x, const int rhs_stride) noexcept {
        int unused = 0;
        return solve<internal_rhs, internal_piv, internal_matrix, false, A_static_stride,
                     piv_static_stride, rhs_static_stride>(A, A_stride, piv, piv_stride, b, x,
                                                           rhs_stride, unused);
    }
 
    template<int nrhs, bool internal_rhs, bool internal_piv, bool internal_matrix,
             int pass_width = 0, bool oot_diagnostics = true, int A_static_stride = dynamic_stride,
             int piv_static_stride = dynamic_stride, int rhs_static_stride = dynamic_stride>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve_multirhs(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
                   const int piv_stride, const T* TDLS_RESTRICT b, T* TDLS_RESTRICT x,
                   const int rhs_stride, const std::size_t xcol_stride, int& oot_count) noexcept {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return solve_multirhs<nrhs, internal_rhs, internal_piv, internal_matrix, pass_width,
                                  oot_diagnostics>(
                A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), b, x,
                detail::stride_or<rhs_static_stride>(rhs_stride), xcol_stride, oot_count);
        } else {
            if (!factorize<internal_piv, internal_matrix, oot_diagnostics>(A, A_stride, piv,
                                                                           piv_stride, oot_count))
                return false;
            substitute_multirhs<nrhs, internal_rhs, internal_piv, internal_matrix, pass_width>(
                A, A_stride, piv, piv_stride, b, x, rhs_stride, xcol_stride);
            return true;
        }
    }
 
    template<int nrhs, bool internal_rhs, bool internal_piv, bool internal_matrix,
             int pass_width = 0, int A_static_stride = dynamic_stride,
             int piv_static_stride = dynamic_stride, int rhs_static_stride = dynamic_stride>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve_multirhs(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
                   const int piv_stride, const T* TDLS_RESTRICT b, T* TDLS_RESTRICT x,
                   const int rhs_stride, const std::size_t xcol_stride) noexcept {
        int unused = 0;
        return solve_multirhs<nrhs, internal_rhs, internal_piv, internal_matrix, pass_width, false,
                              A_static_stride, piv_static_stride, rhs_static_stride>(
            A, A_stride, piv, piv_stride, b, x, rhs_stride, xcol_stride, unused);
    }
 
    template<bool internal_rhs, bool internal_piv, bool internal_matrix,
             bool oot_diagnostics = true, int A_static_stride = dynamic_stride,
             int piv_static_stride = dynamic_stride, int rhs_static_stride = dynamic_stride>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve_inplace(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
                  const int piv_stride, T* TDLS_RESTRICT y, const int rhs_stride,
                  int& oot_count) noexcept {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return solve_inplace<internal_rhs, internal_piv, internal_matrix, oot_diagnostics>(
                A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), y,
                detail::stride_or<rhs_static_stride>(rhs_stride), oot_count);
        } else {
            if constexpr (physical) {
                int weak      = 0;
                const bool ok = physical_factorize<internal_piv, internal_matrix>(A, A_stride, piv,
                                                                                  piv_stride, weak);
                if constexpr (oot_diagnostics) oot_count = weak;
                if (!ok) return false;
                physical_multirhs<1, 0, 3, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, y, rhs_stride, 0);
                return true;
            }
            if (!factorize<internal_piv, internal_matrix, oot_diagnostics, true, internal_rhs>(
                    A, A_stride, piv, piv_stride, oot_count, y, rhs_stride))
                return false;
 
            // Backward pass only: the forward one happened inside factorize.
            if constexpr (transposed) {
                // Column-major layout: L^T w = z, then x = P^T w.
                transposed_bwd_only<1, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, y, rhs_stride, 0);
                unpermute<1, internal_rhs, internal_piv>(piv, piv_stride, y, rhs_stride, 0);
            } else {
                bwd_only<1, internal_rhs, internal_piv, internal_matrix>(
                    A, A_stride, piv, piv_stride, y, rhs_stride, 0);
            }
            return true;
        }
    }
 
    template<bool internal_rhs, bool internal_piv, bool internal_matrix,
             int A_static_stride = dynamic_stride, int piv_static_stride = dynamic_stride,
             int rhs_static_stride = dynamic_stride>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve_inplace(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
                  const int piv_stride, T* TDLS_RESTRICT y, const int rhs_stride) noexcept {
        int unused = 0;
        return solve_inplace<internal_rhs, internal_piv, internal_matrix, false, A_static_stride,
                             piv_static_stride, rhs_static_stride>(A, A_stride, piv, piv_stride, y,
                                                                   rhs_stride, unused);
    }
 
    template<int nrhs, bool internal_rhs, bool internal_piv, bool internal_matrix,
             int pass_width = 0, bool oot_diagnostics = true, int A_static_stride = dynamic_stride,
             int piv_static_stride = dynamic_stride, int rhs_static_stride = dynamic_stride>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve_inplace_multirhs(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
                           const int piv_stride, T* TDLS_RESTRICT y, const int rhs_stride,
                           const std::size_t xcol_stride, int& oot_count) noexcept {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return solve_inplace_multirhs<nrhs, internal_rhs, internal_piv, internal_matrix,
                                          pass_width, oot_diagnostics>(
                A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), y,
                detail::stride_or<rhs_static_stride>(rhs_stride), xcol_stride, oot_count);
        } else {
            if (!factorize<internal_piv, internal_matrix, oot_diagnostics>(A, A_stride, piv,
                                                                           piv_stride, oot_count))
                return false;
            substitute_inplace_multirhs<nrhs, internal_rhs, internal_piv, internal_matrix,
                                        pass_width>(A, A_stride, piv, piv_stride, y, rhs_stride,
                                                    xcol_stride);
            return true;
        }
    }
 
    template<int nrhs, bool internal_rhs, bool internal_piv, bool internal_matrix,
             int pass_width = 0, int A_static_stride = dynamic_stride,
             int piv_static_stride = dynamic_stride, int rhs_static_stride = dynamic_stride>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve_inplace_multirhs(T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
                           const int piv_stride, T* TDLS_RESTRICT y, const int rhs_stride,
                           const std::size_t xcol_stride) noexcept {
        int unused = 0;
        return solve_inplace_multirhs<nrhs, internal_rhs, internal_piv, internal_matrix, pass_width,
                                      false, A_static_stride, piv_static_stride, rhs_static_stride>(
            A, A_stride, piv, piv_stride, y, rhs_stride, xcol_stride, unused);
    }
 
    /* =====================================================================
       Note on pivot ownership: the pivot storage is ALWAYS caller-provided
       (internal array or remote scratch), for every entry point: one
       uniform calling convention. A solver-internal pivot would only be
       expressible for the single combination {solve/solve_inplace x internal
       pivot}; externalizing it there costs nothing (the caller's int
       piv[N] inlines to the exact same codegen), so no special case is
       kept. For one-shot solves the contents are simply treated as
       scratch; for factorize + substitute* (consistent tangent), the same
       storage carries the permutation across the calls.
       ===================================================================== */
};
 
 
 
#undef TDLS_LUPP_A
#undef TDLS_LUPP_PIV
#undef TDLS_LUPP_X
#undef TDLS_LUPP_B
#undef TDLS_LUPP_XW
#undef TDLS_LUPP_BW
#undef TDLS_LUPP_Y
#undef TDLS_LUPP_FOR
 
 
 
} // namespace tdls
 
#if defined(TDLS_DETAIL_COMPILER_GCC)
#pragma GCC diagnostic pop
#endif
 
#if defined(TDLS_DETAIL_COMPILER_CLANG)
#pragma clang diagnostic pop
#endif
 
 
 
#endif // TDLS_SOLVERS_TILED_LUPP_SOLVER_STATIC_HPP