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 <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(__GNUC__) && !defined(__clang__)
#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(__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 ? TDLS_LAYOUT_INDEX(r, c, N)                                                 \
                      : TDLS_LAYOUT_INDEX(r, c, N) * unsigned(A_stride)]
#define TDLS_LUPP_PIV(i) piv[internal_piv ? unsigned(i) : unsigned(i) * unsigned(piv_stride)]
#define TDLS_LUPP_X(i) x[internal_rhs ? unsigned(i) : unsigned(i) * unsigned(rhs_stride)]
#define TDLS_LUPP_B(i) b[internal_rhs ? unsigned(i) : unsigned(i) * unsigned(rhs_stride)]
#define TDLS_LUPP_XW(w, i)                                                                         \
    x[internal_rhs ? unsigned((w) * N + (i))                                                       \
                   : unsigned(i) * unsigned(rhs_stride) + unsigned(w) * unsigned(xcol_stride)]
#define TDLS_LUPP_BW(w, i)                                                                         \
    b[internal_rhs ? unsigned((w) * N + (i))                                                       \
                   : unsigned(i) * unsigned(rhs_stride) + unsigned(w) * unsigned(xcol_stride)]
#define TDLS_LUPP_Y(i) y[internal_rhs ? unsigned(i) : unsigned(i) * unsigned(rhs_stride)]
 
 
 
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_threshold = Config.oot_threshold;
    static constexpr T singular_floor = Config.singular_floor;
 
    static_assert(N >= 1, "TiledLUppSolverStatic: N must be >= 1");
    static_assert(
        Config.oot_threshold.is_finite() && Config.singular_floor.is_finite(),
        "TiledLUppSolverStatic: oot_threshold and singular_floor must be finite (and fit a "
        "63-bit mantissa)");
    static_assert(singular_floor > T(0), "TiledLUppSolverStatic: singular_floor must be positive");
    static_assert(singular_floor <= oot_threshold,
                  "TiledLUppSolverStatic: singular_floor must not exceed oot_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_inner>;
 
    /* =====================================================================
       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_FORCEINLINE 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_inner) {
            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_FORCEINLINE 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_inner) {
            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_FORCEINLINE 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_inner) {
            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_FORCEINLINE 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_inner) {
            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_FORCEINLINE 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_inner) {
            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_FORCEINLINE 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_inner) {
            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_FORCEINLINE 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_inner) {
            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;
                }
            }
        }
 
        int piv_row; // winning global (logical) row
 
        if (best >= oot_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_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 (lanes of a warp pick
            // different pivots almost every column).
            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) {
                // 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.
                if constexpr (internal_rhs) {
                    if constexpr (Config.unroll_inner) {
                        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_inner) {
                    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_FORCEINLINE 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_inner) {
            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 (they fold to nothing when the
       permutation itself is internal).
       ===================================================================== */
 
    template<int k_extent, int j_extent, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 i_extent, int j_extent, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCEINLINE 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_inner) {
            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_FORCEINLINE 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_inner) {
            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).
        if constexpr (fuse_rhs) {
            if constexpr (Config.unroll_inner) {
                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_FORCEINLINE 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_inner) {
            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: unit-lower-solve it while the tile is in registers.
        if constexpr (fuse_rhs) {
            if constexpr (Config.unroll_inner) {
                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
        for (int j = k + 1; j < full_tiles; ++j)
            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)
            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_FORCEINLINE 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_FORCEINLINE 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.
        if constexpr (fuse_rhs) {
            if constexpr (Config.unroll_inner) {
                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_FORCEINLINE 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, bool internal_piv, bool internal_matrix, bool oot_diagnostics,
             bool fuse_rhs = false, bool internal_rhs = true>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE 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: unit-lower-solve it while the tile is in registers.
        if constexpr (fuse_rhs) {
            if constexpr (Config.unroll_inner) {
                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
        for (int j = k + 1; j < full_tiles; ++j)
            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)
            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;
    }
 
    /* =====================================================================
       FACTORIZE - public entry point.
       ===================================================================== */
 
    template<bool internal_piv, bool internal_matrix, bool oot_diagnostics = true,
             bool fuse_rhs = false, bool internal_rhs = true>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 (oot_diagnostics) oot_count = 0;
 
        if constexpr (Config.unroll_inner) {
            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>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE 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>(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_FORCEINLINE 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 int 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_inner) {
            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_FORCEINLINE 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 int 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_inner) {
            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_FORCEINLINE 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 int 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_inner) {
            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_FORCEINLINE 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 int 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_inner) {
            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_FORCEINLINE 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 int 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_FORCEINLINE 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 int 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);
    }
 
    /* =====================================================================
       SUBSTITUTION - public entry points.
       ===================================================================== */
 
    template<bool internal_rhs, bool internal_piv, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 (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_inner) {
                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_inner) {
                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>
    TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 {
        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_FORCEINLINE 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 int xcol_stride) noexcept {
        if constexpr (Config.unroll_inner) {
            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>
    TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 int xcol_stride) noexcept {
        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 unsigned off =
                internal_rhs ? unsigned(c0) * unsigned(N) : unsigned(c0) * unsigned(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 unsigned off =
                internal_rhs ? unsigned(c0) * unsigned(N) : unsigned(c0) * unsigned(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_FORCEINLINE 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 int xcol_stride) noexcept {
        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_inner) {
                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_inner) {
                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>
    TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 int xcol_stride) noexcept {
        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 unsigned off =
                internal_rhs ? unsigned(c0) * unsigned(N) : unsigned(c0) * unsigned(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 unsigned off =
                internal_rhs ? unsigned(c0) * unsigned(N) : unsigned(c0) * unsigned(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>
    TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 (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_inner) {
                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, bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 int 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_inner) {
                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];
            }
        }
        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>
    TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 int xcol_stride) noexcept {
        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 unsigned off =
                internal_rhs ? unsigned(c0) * unsigned(N) : unsigned(c0) * unsigned(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 unsigned off =
                internal_rhs ? unsigned(c0) * unsigned(N) : unsigned(c0) * unsigned(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);
        }
    }
 
    /* =====================================================================
       SOLVE - convenience wrappers.
       ===================================================================== */
 
    template<bool internal_rhs, bool internal_piv, bool internal_matrix,
             bool oot_diagnostics = true>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 (!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>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE 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, 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>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 int xcol_stride, int& oot_count) noexcept {
        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>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 int xcol_stride) noexcept {
        int unused = 0;
        return solve_multirhs<nrhs, internal_rhs, internal_piv, internal_matrix, pass_width, false>(
            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>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 (!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.
        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>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE 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, 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>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 int xcol_stride, int& oot_count) noexcept {
        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>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 int xcol_stride) noexcept {
        int unused = 0;
        return solve_inplace_multirhs<nrhs, internal_rhs, internal_piv, internal_matrix, pass_width,
                                      false>(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
 
 
 
} // namespace tdls
 
#if defined(__GNUC__) && !defined(__clang__)
#pragma GCC diagnostic pop
#endif
 
#if defined(__clang__)
#pragma clang diagnostic pop
#endif
 
 
 
#endif // TDLS_SOLVERS_TILED_LUPP_SOLVER_STATIC_HPP