solver_dynamic.hpp

Go to the documentation of this file.

#ifndef TDLS_SOLVERS_TILED_LUPP_SOLVER_DYNAMIC_HPP
#define TDLS_SOLVERS_TILED_LUPP_SOLVER_DYNAMIC_HPP
 
 
 
 
 
 
#include <tdls/core/math.hpp>
#include <tdls/core/target.hpp>
#include <type_traits>
 
#include <tdls/solvers/tiled_lupp/config.hpp>
 
 
 
// gcc's flow analysis cannot prove that the replay buffer of the
// out-of-tile search is written before being read (the loop structure
// guarantees it, a property validated bitwise against the static
// solver), and emits spurious -Wmaybe-uninitialized warnings at -O2.
//
// 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_static.hpp.
#if defined(TDLS_DETAIL_COMPILER_GCC)
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
#pragma GCC diagnostic ignored "-Warray-bounds"
#endif
 
namespace tdls {
 
 
 
/* Addressing macros. They expand inside member functions where n, the
   parameter names (A, A_stride, piv, piv_stride, x, b, y, rhs_stride,
   xcol_stride) and the Config value are in scope. #undef'd at the end
   of this header. */
 
#define TDLS_LUPP_DYN_A(r, c) A[TDLS_STRIDED((r) * n + (c), A_stride)]
#define TDLS_LUPP_DYN_PIV(i) piv[TDLS_STRIDED(i, piv_stride)]
#define TDLS_LUPP_DYN_X(i) x[TDLS_STRIDED(i, rhs_stride)]
#define TDLS_LUPP_DYN_B(i) b[TDLS_STRIDED(i, rhs_stride)]
#define TDLS_LUPP_DYN_XW(w, i) x[TDLS_STRIDED(i, rhs_stride) + TDLS_STRIDED_WIDE(w, xcol_stride)]
#define TDLS_LUPP_DYN_BW(w, i) b[TDLS_STRIDED(i, rhs_stride) + TDLS_STRIDED_WIDE(w, xcol_stride)]
#define TDLS_LUPP_DYN_Y(i) y[TDLS_STRIDED(i, rhs_stride)]
 
 
 
template<typename T, TiledLUppConfig<T> Config = TiledLUppConfig<T>{}>
struct TiledLUppSolverDynamic {
 
    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(Config.row_interchange != RowInterchange::Auto,
                  "TiledLUppSolverDynamic: 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(),
        "TiledLUppSolverDynamic: oot_pivot_threshold and singular_floor must be finite (and fit a "
        "63-bit mantissa)");
    static_assert(Config.relative_pivot_threshold.is_finite(),
                  "TiledLUppSolverDynamic: relative_pivot_threshold must be finite (and fit a "
                  "63-bit mantissa)");
    static_assert(relative_pivot_threshold > T(0) && relative_pivot_threshold <= T(1),
                  "TiledLUppSolverDynamic: relative_pivot_threshold must be in (0, 1]");
    static_assert(singular_floor > T(0), "TiledLUppSolverDynamic: singular_floor must be positive");
    static_assert(singular_floor <= oot_pivot_threshold,
                  "TiledLUppSolverDynamic: singular_floor must not exceed oot_pivot_threshold (the "
                  "floor applies to the out-of-tile recovery path)");
    static_assert(tile_size >= 1, "TiledLUppSolverDynamic: tile size must be >= 1");
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr int num_tiles(const int n) noexcept {
        return (n + tile_size - 1) / tile_size;
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr int tile_size_at(const int t0,
                                                                             const int n) noexcept {
        return (n - t0 < tile_size) ? (n - t0) : tile_size;
    }
 
    /* =====================================================================
       Register-tile micro-kernels, runtime extents. A tile is a
       tile_size x tile_size array with row stride tile_size; the runtime
       extents bound the active part.
       ===================================================================== */
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    ops_swap_rows(T* TDLS_RESTRICT t, const int k, const int r, const int ke) noexcept {
        for (int j = 0; j < ke; ++j) {
            const T tmp          = t[k * tile_size + j];
            t[k * tile_size + j] = t[r * tile_size + j];
            t[r * tile_size + j] = tmp;
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    ops_eliminate_column(T* TDLS_RESTRICT t, const int k, const int re, const int ce) noexcept {
        const T inv_pivot    = T(1) / t[k * tile_size + k];
        t[k * tile_size + k] = inv_pivot;
        for (int i = k + 1; i < re; ++i) {
            t[i * tile_size + k] *= inv_pivot;
            for (int j = k + 1; j < ce; ++j)
                t[i * tile_size + j] -= t[i * tile_size + k] * t[k * tile_size + j];
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    ops_trsm_left_unit(const T* TDLS_RESTRICT lu, T* TDLS_RESTRICT B, const int kd,
                       const int ce) noexcept {
        for (int k = 0; k < kd; ++k) {
            for (int i = k + 1; i < kd; ++i) {
                const T L_ik = lu[i * tile_size + k];
                for (int j = 0; j < ce; ++j)
                    B[i * tile_size + j] -= L_ik * B[k * tile_size + j];
            }
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    ops_trsm_right(const T* TDLS_RESTRICT lu, T* TDLS_RESTRICT B, const int kd,
                   const int re) noexcept {
        for (int k = 0; k < kd; ++k) {
            const T U_kk_inv = lu[k * tile_size + k];
            for (int i = 0; i < re; ++i)
                B[i * tile_size + k] *= U_kk_inv;
            for (int j = k + 1; j < kd; ++j) {
                const T U_kj = lu[k * tile_size + j];
                for (int i = 0; i < re; ++i)
                    B[i * tile_size + j] -= B[i * tile_size + k] * U_kj;
            }
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    ops_gemm_sub(T* TDLS_RESTRICT Ct, const T* TDLS_RESTRICT At, const T* TDLS_RESTRICT Bt,
                 const int re, const int ce, const int kd) noexcept {
        for (int i = 0; i < re; ++i) {
            for (int j = 0; j < ce; ++j) {
                T sum = T(0);
                for (int k = 0; k < kd; ++k)
                    sum += At[i * tile_size + k] * Bt[k * tile_size + j];
                Ct[i * tile_size + j] -= sum;
            }
        }
    }
 
    /* =====================================================================
       Remote <-> register tile movement, runtime extents.
       ===================================================================== */
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_tile(const int n, const T* TDLS_RESTRICT A, const int A_stride,
              const int* TDLS_RESTRICT prow, const int col0, T* TDLS_RESTRICT t, const int re,
              const int ce) noexcept {
        for (int i = 0; i < re; ++i) {
            for (int j = 0; j < ce; ++j)
                t[i * tile_size + j] = TDLS_LUPP_DYN_A(prow[i], col0 + j);
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    store_tile(const int n, T* TDLS_RESTRICT A, const int A_stride, const int* TDLS_RESTRICT prow,
               const int col0, const T* TDLS_RESTRICT t, const int re, const int ce) noexcept {
        for (int i = 0; i < re; ++i) {
            for (int j = 0; j < ce; ++j)
                TDLS_LUPP_DYN_A(prow[i], col0 + j) = t[i * tile_size + j];
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_tile_piv(const int n, 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, const int re, const int ce) noexcept {
        for (int i = 0; i < re; ++i) {
            const int phys = TDLS_LUPP_DYN_PIV(row0 + i);
            for (int j = 0; j < ce; ++j)
                t[i * tile_size + j] = TDLS_LUPP_DYN_A(phys, col0 + j);
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    store_tile_piv(const int n, 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, const int re, const int ce) noexcept {
        for (int i = 0; i < re; ++i) {
            const int phys = TDLS_LUPP_DYN_PIV(row0 + i);
            for (int j = 0; j < ce; ++j)
                TDLS_LUPP_DYN_A(phys, col0 + j) = t[i * tile_size + j];
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_tile_piv_lower(const int n, 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, const int re) noexcept {
        for (int i = 1; i < re; ++i) {
            const int phys = TDLS_LUPP_DYN_PIV(row0 + i);
            for (int j = 0; j < i; ++j)
                t[i * tile_size + j] = TDLS_LUPP_DYN_A(phys, col0 + j);
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_tile_piv_upper(const int n, 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, const int re) noexcept {
        for (int i = 0; i < re; ++i) {
            const int phys = TDLS_LUPP_DYN_PIV(row0 + i);
            for (int j = i; j < re; ++j)
                t[i * tile_size + j] = TDLS_LUPP_DYN_A(phys, col0 + j);
        }
    }
 
    /* =====================================================================
       Diagonal-tile factorization with out-of-tile pivoting. Same
       algorithm as the static TiledLUpp solver, runtime column/extent parameters,
       direct (unpredicated) fused-RHS swap.
       ===================================================================== */
 
    template<bool oot_diagnostics, bool fuse_rhs>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    factor_diag_column(const int n, 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, const int ke, T* TDLS_RESTRICT y,
                       const int rhs_stride) noexcept {
 
        const int gc = k0 + c; // global column
 
        // In-tile pivot search (rows c..ke of the register tile)
        int best_r = c;
        T best     = detail::abs(tile[c * tile_size + c]);
        for (int r = c + 1; r < ke; ++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 (ke < 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).
            if constexpr (oot_diagnostics) ++oot_count;
            T gbest       = best;
            int gbest_row = k0 + best_r;
 
            for (int row = k0 + ke; row < n; ++row) {
                const int phys = TDLS_LUPP_DYN_PIV(row);
 
                T corrected = TDLS_LUPP_DYN_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_DYN_A(phys, bj0 + p) *
                                         TDLS_LUPP_DYN_A(TDLS_LUPP_DYN_PIV(bj0 + p), gc);
                }
 
                if (c > 0) {
                    T L_row[tile_size];
                    for (int t = 0; t < c; ++t) {
                        T a_t = TDLS_LUPP_DYN_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_DYN_A(phys, bj0 + p) *
                                           TDLS_LUPP_DYN_A(TDLS_LUPP_DYN_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.
                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 (self-swap is a
            // bitwise no-op, and the branch would diverge on GPU).
            const int tmp              = TDLS_LUPP_DYN_PIV(gc);
            TDLS_LUPP_DYN_PIV(gc)      = TDLS_LUPP_DYN_PIV(piv_row);
            TDLS_LUPP_DYN_PIV(piv_row) = tmp;
 
            if constexpr (fuse_rhs && !transposed) {
                // The fused RHS follows the rows through pivoting (row-major
                // layout; under the column-major one y is indexed by the
                // columns of S, which pivoting never moves).
                const T ty               = TDLS_LUPP_DYN_Y(gc);
                TDLS_LUPP_DYN_Y(gc)      = TDLS_LUPP_DYN_Y(piv_row);
                TDLS_LUPP_DYN_Y(piv_row) = ty;
            }
 
            if (piv_row < k0 + ke) {
                ops_swap_rows(tile, c, piv_row - k0, ke);
            } 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_DYN_PIV(gc);
                for (int j = 0; j < ke; ++j) {
                    tile[c * tile_size + j] = TDLS_LUPP_DYN_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_DYN_A(phys, bj0 + p) *
                                    TDLS_LUPP_DYN_A(TDLS_LUPP_DYN_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 < ke; ++j) {
                        for (int t = 0; t < c; ++t)
                            tile[c * tile_size + j] -= L_row[t] * tile[t * tile_size + j];
                    }
                }
            }
        }
 
        ops_eliminate_column(tile, c, ke, ke);
        return true;
    }
 
    template<bool oot_diagnostics, bool fuse_rhs>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    factor_diag_tile(const int n, 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 ke, T* TDLS_RESTRICT y = nullptr,
                     const int rhs_stride = 1) noexcept {
        for (int c = 0; c < ke; ++c) {
            if (!factor_diag_column<oot_diagnostics, fuse_rhs>(
                    n, A, A_stride, piv, piv_stride, k0, tile, oot_count, c, ke, y, rhs_stride))
                return false;
        }
        return true;
    }
 
    /* =====================================================================
       RIGHT-LOOKING schedule.
       ===================================================================== */
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    rl_trsm_right_one(const int n, T* TDLS_RESTRICT A, const int A_stride,
                      const int* TDLS_RESTRICT pk, const T* TDLS_RESTRICT tile, const int j0,
                      const int ke, const int je) noexcept {
        T Akj[tile_size * tile_size];
        load_tile(n, A, A_stride, pk, j0, Akj, ke, je);
        ops_trsm_left_unit(tile, Akj, ke, je);
        store_tile(n, A, A_stride, pk, j0, Akj, ke, je);
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    rl_trsm_right_push_one(const int n, T* TDLS_RESTRICT A, const int A_stride,
                           const int* TDLS_RESTRICT pk, const T* TDLS_RESTRICT tile, const int j0,
                           const int ke, const int je, const int k0, T* TDLS_RESTRICT y,
                           const int rhs_stride) noexcept {
        T Akj[tile_size * tile_size];
        load_tile(n, A, A_stride, pk, j0, Akj, ke, je);
        ops_trsm_left_unit(tile, Akj, ke, je);
        store_tile(n, A, A_stride, pk, j0, Akj, ke, je);
        transposed_push_one(Akj, y, rhs_stride, 0, 1, k0, j0, ke, je);
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    rl_schur_one(const int n, 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,
                 const int ke, const int ie, const int je) noexcept {
        T Aij[tile_size * tile_size];
        load_tile(n, A, A_stride, pi, j0, Aij, ie, je);
 
        for (int p = 0; p < ke; ++p) {
            T Akj_row[tile_size];
            for (int j = 0; j < je; ++j)
                Akj_row[j] = TDLS_LUPP_DYN_A(pk[p], j0 + j);
            for (int i = 0; i < ie; ++i) {
                const T L_ip = Aik[i * tile_size + p];
                for (int j = 0; j < je; ++j)
                    Aij[i * tile_size + j] -= L_ip * Akj_row[j];
            }
        }
 
        store_tile(n, A, A_stride, pi, j0, Aij, ie, je);
    }
 
    template<bool fuse_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    rl_update_row_one(const int n, 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, const int ke, const int ie, T* TDLS_RESTRICT y = nullptr,
                      const int rhs_stride = 1) noexcept {
        const int k0 = k * tile_size;
        const int nt = num_tiles(n);
 
        int pi[tile_size];
        for (int i = 0; i < ie; ++i)
            pi[i] = TDLS_LUPP_DYN_PIV(i0 + i);
 
        // TRSM down: Aik := Aik * U^-1
        T Aik[tile_size * tile_size];
        load_tile(n, A, A_stride, pi, k0, Aik, ie, ke);
        ops_trsm_right(tile, Aik, ke, ie);
        store_tile(n, A, A_stride, pi, k0, Aik, ie, ke);
 
        // Fused forward substitution: push the solved y_k segment into
        // this row block while its L panel sits in registers (row-major
        // layout; the column-major one pushes from the U row panel).
        if constexpr (fuse_rhs && !transposed) {
            for (int r = 0; r < ie; ++r) {
                T sum = T(0);
                for (int j = 0; j < ke; ++j)
                    sum += Aik[r * tile_size + j] * TDLS_LUPP_DYN_Y(k0 + j);
                TDLS_LUPP_DYN_Y(i0 + r) -= sum;
            }
        }
 
        // Schur sweep over the trailing columns
        for (int tj = k + 1; tj < nt; ++tj)
            rl_schur_one(n, A, A_stride, pk, pi, Aik, tj * tile_size, ke, ie,
                         tile_size_at(tj * tile_size, n));
    }
 
    template<bool oot_diagnostics, bool fuse_rhs>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    rl_step(const int n, 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;
        const int nt = num_tiles(n);
        const int ke = tile_size_at(k0, n);
 
        T tile[tile_size * tile_size];
        load_tile_piv(n, A, A_stride, piv, piv_stride, k0, k0, tile, ke, ke);
 
        if (!factor_diag_tile<oot_diagnostics, fuse_rhs>(n, A, A_stride, piv, piv_stride, k0, tile,
                                                         oot_count, ke, y, rhs_stride))
            return false;
 
        // Physical rows of the tile after the swaps of this step
        int pk[tile_size];
        for (int i = 0; i < ke; ++i)
            pk[i] = TDLS_LUPP_DYN_PIV(k0 + i);
 
        store_tile(n, A, A_stride, pk, k0, tile, ke, ke);
 
        // 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 && transposed) {
            transposed_solve_tile_upper(tile, y, rhs_stride, 0, 1, k0, ke);
        } else if constexpr (fuse_rhs) {
            for (int kk = 0; kk < ke; ++kk) {
                for (int i = kk + 1; i < ke; ++i)
                    TDLS_LUPP_DYN_Y(k0 + i) -= tile[i * tile_size + kk] * TDLS_LUPP_DYN_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 tj = k + 1; tj < nt; ++tj) {
            if constexpr (fuse_rhs && transposed)
                rl_trsm_right_push_one(n, A, A_stride, pk, tile, tj * tile_size, ke,
                                       tile_size_at(tj * tile_size, n), k0, y, rhs_stride);
            else
                rl_trsm_right_one(n, A, A_stride, pk, tile, tj * tile_size, ke,
                                  tile_size_at(tj * tile_size, n));
        }
 
        // TRSM down + Schur over the rows below
        for (int ti = k + 1; ti < nt; ++ti)
            rl_update_row_one<fuse_rhs>(n, A, A_stride, piv, piv_stride, pk, tile, k,
                                        ti * tile_size, ke, tile_size_at(ti * tile_size, n), y,
                                        rhs_stride);
        return true;
    }
 
    /* =====================================================================
       LEFT-LOOKING schedule.
       ===================================================================== */
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    ll_correct_tile(const int n, 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, const int re,
                    const int ce) noexcept {
        for (int bj0 = 0; bj0 < k0; bj0 += tile_size) {
            T Lt[tile_size * tile_size];
            load_tile_piv(n, A, A_stride, piv, piv_stride, row0, bj0, Lt, re, tile_size);
            T Ut[tile_size * tile_size];
            load_tile_piv(n, A, A_stride, piv, piv_stride, bj0, col0, Ut, tile_size, ce);
            ops_gemm_sub(t, Lt, Ut, re, ce, tile_size);
        }
    }
 
    template<bool fuse_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void ll_update_below_one(
        const int n, 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, const int ke,
        const int ie, T* TDLS_RESTRICT y = nullptr, const int rhs_stride = 1) noexcept {
        T B[tile_size * tile_size];
        load_tile_piv(n, A, A_stride, piv, piv_stride, i0, k0, B, ie, ke);
        ll_correct_tile(n, A, A_stride, piv, piv_stride, i0, k0, k0, B, ie, ke);
        ops_trsm_right(tile, B, ke, ie);
        store_tile_piv(n, A, A_stride, piv, piv_stride, i0, k0, B, ie, ke);
 
        // Fused forward substitution: B is the final L(i,k) panel, so push
        // the solved y_k segment into this row block (row-major layout; the
        // column-major one pushes from the U row panel).
        if constexpr (fuse_rhs && !transposed) {
            for (int r = 0; r < ie; ++r) {
                T sum = T(0);
                for (int j = 0; j < ke; ++j)
                    sum += B[r * tile_size + j] * TDLS_LUPP_DYN_Y(k0 + j);
                TDLS_LUPP_DYN_Y(i0 + r) -= sum;
            }
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    ll_update_right_one(const int n, 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, const int ke,
                        const int je) noexcept {
        T B[tile_size * tile_size];
        load_tile_piv(n, A, A_stride, piv, piv_stride, k0, j0, B, ke, je);
        ll_correct_tile(n, A, A_stride, piv, piv_stride, k0, j0, k0, B, ke, je);
        ops_trsm_left_unit(tile, B, ke, je);
        store_tile_piv(n, A, A_stride, piv, piv_stride, k0, j0, B, ke, je);
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    ll_update_right_push_one(const int n, 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, const int ke,
                             const int je, T* TDLS_RESTRICT y, const int rhs_stride) noexcept {
        T B[tile_size * tile_size];
        load_tile_piv(n, A, A_stride, piv, piv_stride, k0, j0, B, ke, je);
        ll_correct_tile(n, A, A_stride, piv, piv_stride, k0, j0, k0, B, ke, je);
        ops_trsm_left_unit(tile, B, ke, je);
        store_tile_piv(n, A, A_stride, piv, piv_stride, k0, j0, B, ke, je);
        transposed_push_one(B, y, rhs_stride, 0, 1, k0, j0, ke, je);
    }
 
    template<bool oot_diagnostics, bool fuse_rhs>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    ll_step(const int n, 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;
        const int nt = num_tiles(n);
        const int ke = tile_size_at(k0, n);
 
        T tile[tile_size * tile_size];
        load_tile_piv(n, A, A_stride, piv, piv_stride, k0, k0, tile, ke, ke);
        ll_correct_tile(n, A, A_stride, piv, piv_stride, k0, k0, k0, tile, ke, ke);
 
        if (!factor_diag_tile<oot_diagnostics, fuse_rhs>(n, A, A_stride, piv, piv_stride, k0, tile,
                                                         oot_count, ke, y, rhs_stride))
            return false;
 
        store_tile_piv(n, A, A_stride, piv, piv_stride, k0, k0, tile, ke, ke);
 
        // 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 && transposed) {
            transposed_solve_tile_upper(tile, y, rhs_stride, 0, 1, k0, ke);
        } else if constexpr (fuse_rhs) {
            for (int kk = 0; kk < ke; ++kk) {
                for (int i = kk + 1; i < ke; ++i)
                    TDLS_LUPP_DYN_Y(k0 + i) -= tile[i * tile_size + kk] * TDLS_LUPP_DYN_Y(k0 + kk);
            }
        }
 
        // L panel below the diagonal
        for (int ti = k + 1; ti < nt; ++ti)
            ll_update_below_one<fuse_rhs>(n, A, A_stride, piv, piv_stride, tile, k0, ti * tile_size,
                                          ke, tile_size_at(ti * tile_size, n), 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 tj = k + 1; tj < nt; ++tj) {
            if constexpr (fuse_rhs && transposed)
                ll_update_right_push_one(n, A, A_stride, piv, piv_stride, tile, k0, tj * tile_size,
                                         ke, tile_size_at(tj * tile_size, n), y, rhs_stride);
            else
                ll_update_right_one(n, A, A_stride, piv, piv_stride, tile, k0, tj * tile_size, ke,
                                    tile_size_at(tj * tile_size, n));
        }
        return true;
    }
 
    /* =====================================================================
       Physical row interchanges (Config.row_interchange == Physical): the
       whole matrix is one tile, the rows exchanged as in LAPACK; the
       operations of TiledLUppSolverStatic in the same order, so that both
       solvers agree bit for bit.
       ===================================================================== */
 
    template<bool check>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr int
    physical_factorize_impl(const int n, T* TDLS_RESTRICT A, const int A_stride,
                            int* TDLS_RESTRICT piv, const int piv_stride,
                            [[maybe_unused]] const T u, int& weak) noexcept {
        bool ok       = true;
        bool singular = false;
        weak          = 0;
        for (int j = 0; j < n; ++j) {
            T best = detail::abs(TDLS_LUPP_DYN_A(j, j));
            int r  = j;
            for (int i = j + 1; i < n; ++i) {
                const T a = detail::abs(TDLS_LUPP_DYN_A(i, j));
                if (a > best) {
                    best = a;
                    r    = i;
                }
            }
            if constexpr (relative_pivot_threshold < T(1)) {
                const T in_place = detail::abs(TDLS_LUPP_DYN_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_DYN_PIV(j) = r;
            if (r != j) {
                for (int c = 0; c < n; ++c) {
                    const T t             = TDLS_LUPP_DYN_A(j, c);
                    TDLS_LUPP_DYN_A(j, c) = TDLS_LUPP_DYN_A(r, c);
                    TDLS_LUPP_DYN_A(r, c) = t;
                }
            }
            const T inv           = pivot_ok ? T(1) / TDLS_LUPP_DYN_A(j, j) : T(0);
            TDLS_LUPP_DYN_A(j, j) = inv;
            for (int i = j + 1; i < n; ++i) {
                const T l             = TDLS_LUPP_DYN_A(i, j) * inv;
                TDLS_LUPP_DYN_A(i, j) = l;
                if constexpr (check) ok = ok && u * detail::abs(l) <= T(1);
                for (int c = j + 1; c < n; ++c)
                    TDLS_LUPP_DYN_A(i, c) -= l * TDLS_LUPP_DYN_A(j, c);
            }
        }
        return singular ? 2 : ok ? 0 : 1;
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    physical_factorize(const int n, T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
                       const int piv_stride, int& weak) noexcept {
        return physical_factorize_impl<false>(n, A, A_stride, piv, piv_stride, T(0), weak) != 2;
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr int
    physical_factorize_checked(const int n, T* TDLS_RESTRICT A, const int A_stride,
                               int* TDLS_RESTRICT piv, const int piv_stride, const T u) noexcept {
        int weak = 0;
        return physical_factorize_impl<true>(n, A, A_stride, piv, piv_stride, u, weak);
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    physical_permute(const int n, const int w, const int* TDLS_RESTRICT piv, const int piv_stride,
                     T* TDLS_RESTRICT x, const int rhs_stride, const std::size_t xcol_stride,
                     const bool reverse) noexcept {
        for (int t = 0; t < n; ++t) {
            const int j = reverse ? n - 1 - t : t;
            const int r = TDLS_LUPP_DYN_PIV(j);
            if (r != j)
                for (int c = 0; c < w; ++c) {
                    const T v              = TDLS_LUPP_DYN_XW(c, j);
                    TDLS_LUPP_DYN_XW(c, j) = TDLS_LUPP_DYN_XW(c, r);
                    TDLS_LUPP_DYN_XW(c, r) = v;
                }
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    physical_forward(const int n, const int w, 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) {
            physical_permute(n, w, piv, piv_stride, x, rhs_stride, xcol_stride, false);
            for (int j = 0; j < n; ++j)
                for (int i = j + 1; i < n; ++i) {
                    const T l = TDLS_LUPP_DYN_A(i, j);
                    for (int c = 0; c < w; ++c)
                        TDLS_LUPP_DYN_XW(c, i) -= l * TDLS_LUPP_DYN_XW(c, j);
                }
        } else {
            for (int j = 0; j < n; ++j) {
                for (int c = 0; c < w; ++c)
                    TDLS_LUPP_DYN_XW(c, j) *= TDLS_LUPP_DYN_A(j, j);
                for (int q = j + 1; q < n; ++q) {
                    const T uq = TDLS_LUPP_DYN_A(j, q);
                    for (int c = 0; c < w; ++c)
                        TDLS_LUPP_DYN_XW(c, q) -= uq * TDLS_LUPP_DYN_XW(c, j);
                }
            }
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    physical_backward(const int n, const int w, 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 j = n - 1; j >= 0; --j)
                for (int c = 0; c < w; ++c) {
                    T v = TDLS_LUPP_DYN_XW(c, j);
                    for (int q = j + 1; q < n; ++q)
                        v -= TDLS_LUPP_DYN_A(j, q) * TDLS_LUPP_DYN_XW(c, q);
                    TDLS_LUPP_DYN_XW(c, j) = v * TDLS_LUPP_DYN_A(j, j);
                }
        } else {
            for (int j = n - 1; j >= 0; --j)
                for (int q = j + 1; q < n; ++q) {
                    const T l = TDLS_LUPP_DYN_A(q, j);
                    for (int c = 0; c < w; ++c)
                        TDLS_LUPP_DYN_XW(c, j) -= l * TDLS_LUPP_DYN_XW(c, q);
                }
            physical_permute(n, w, piv, piv_stride, x, rhs_stride, xcol_stride, true);
        }
    }
 
    template<int pass_width, int stages>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    physical_multirhs(const int n, const int nrhs, 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 {
        const int cpp = (pass_width <= 0 || pass_width >= nrhs) ? nrhs : pass_width;
        for (int c0 = 0; c0 < nrhs; c0 += cpp) {
            const int w         = nrhs - c0 < cpp ? nrhs - c0 : cpp;
            T* TDLS_RESTRICT xp = x + TDLS_STRIDED_WIDE(c0, xcol_stride);
            if constexpr ((stages & 1) != 0)
                physical_forward(n, w, A, A_stride, piv, piv_stride, xp, rhs_stride, xcol_stride);
            if constexpr ((stages & 2) != 0)
                physical_backward(n, w, A, A_stride, piv, piv_stride, xp, rhs_stride, xcol_stride);
        }
    }
 
    /* =====================================================================
       FACTORIZE - public entry point.
       ===================================================================== */
 
    template<bool oot_diagnostics = true, bool fuse_rhs = false>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    factorize(int n, 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 {
        n = detail::runtime_dimension(n);
        if constexpr (physical) {
            int weak      = 0;
            const bool ok = physical_factorize(n, A, A_stride, piv, piv_stride, weak);
            if constexpr (oot_diagnostics) oot_count = weak;
            if (!ok) return false;
            if constexpr (fuse_rhs)
                physical_forward(n, 1, A, A_stride, piv, piv_stride, y, rhs_stride, 0);
            return true;
        }
 
        if constexpr (oot_diagnostics) oot_count = 0;
 
        for (int i = 0; i < n; ++i)
            TDLS_LUPP_DYN_PIV(i) = i;
 
        const int nt = num_tiles(n);
        for (int k = 0; k < nt; ++k) {
            if constexpr (schedule == Schedule::RightLooking) {
                if (!rl_step<oot_diagnostics, fuse_rhs>(n, A, A_stride, piv, piv_stride, k,
                                                        oot_count, y, rhs_stride))
                    return false;
            } else {
                if (!ll_step<oot_diagnostics, fuse_rhs>(n, A, A_stride, piv, piv_stride, k,
                                                        oot_count, y, rhs_stride))
                    return false;
            }
        }
 
        return true;
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    factorize(int n, T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
              const int piv_stride) noexcept {
        n          = detail::runtime_dimension(n);
        int unused = 0;
        return factorize<false>(n, A, A_stride, piv, piv_stride, unused);
    }
 
    /* =====================================================================
       SUBSTITUTION - schedule-independent.
       ===================================================================== */
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    fwd_push_one(const int n, const int nrhs, 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,
                 const int ke, const int me) noexcept {
        T Lmk[tile_size * tile_size];
        load_tile_piv(n, A, A_stride, piv, piv_stride, m0, k0, Lmk, me, ke);
        for (int i = 0; i < me; ++i) {
            for (int w = 0; w < nrhs; ++w) {
                T sum = T(0);
                for (int j = 0; j < ke; ++j)
                    sum += Lmk[i * tile_size + j] * TDLS_LUPP_DYN_XW(w, k0 + j);
                TDLS_LUPP_DYN_XW(w, m0 + i) -= sum;
            }
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    fwd_step(const int n, const int nrhs, 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;
        const int nt = num_tiles(n);
        const int ke = tile_size_at(k0, n);
 
        T Lkk[tile_size * tile_size];
        load_tile_piv_lower(n, A, A_stride, piv, piv_stride, k0, k0, Lkk, ke);
 
        // In-tile unit-lower solve
        for (int kk = 0; kk < ke; ++kk) {
            for (int i = kk + 1; i < ke; ++i) {
                for (int w = 0; w < nrhs; ++w)
                    TDLS_LUPP_DYN_XW(w, k0 + i) -=
                        Lkk[i * tile_size + kk] * TDLS_LUPP_DYN_XW(w, k0 + kk);
            }
        }
 
        // Push into the tiles below
        for (int m = k + 1; m < nt; ++m)
            fwd_push_one(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride, k0,
                         m * tile_size, ke, tile_size_at(m * tile_size, n));
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    bwd_pull_one(const int n, const int nrhs, 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,
                 const int ke, const int me) noexcept {
        T Ukm[tile_size * tile_size];
        load_tile_piv(n, A, A_stride, piv, piv_stride, k0, m0, Ukm, ke, me);
        for (int i = 0; i < ke; ++i) {
            for (int w = 0; w < nrhs; ++w) {
                T sum = T(0);
                for (int j = 0; j < me; ++j)
                    sum += Ukm[i * tile_size + j] * TDLS_LUPP_DYN_XW(w, m0 + j);
                TDLS_LUPP_DYN_XW(w, k0 + i) -= sum;
            }
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    bwd_step(const int n, const int nrhs, 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;
        const int nt = num_tiles(n);
        const int ke = tile_size_at(k0, n);
 
        // Pull the trailing contributions
        for (int m = k + 1; m < nt; ++m)
            bwd_pull_one(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride, k0,
                         m * tile_size, ke, tile_size_at(m * tile_size, n));
 
        // In-tile upper solve
        T Ukk[tile_size * tile_size];
        load_tile_piv_upper(n, A, A_stride, piv, piv_stride, k0, k0, Ukk, ke);
 
        for (int kk = ke - 1; kk >= 0; --kk) {
            for (int w = 0; w < nrhs; ++w)
                TDLS_LUPP_DYN_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 < nrhs; ++w)
                    TDLS_LUPP_DYN_XW(w, k0 + i) -=
                        Ukk[i * tile_size + kk] * TDLS_LUPP_DYN_XW(w, k0 + kk);
            }
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    bwd_only(const int n, const int nrhs, 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 = num_tiles(n) - 1; k >= 0; --k)
            bwd_step(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride, k);
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    fwd_bwd(const int n, const int nrhs, 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 {
        const int nt = num_tiles(n);
        for (int k = 0; k < nt; ++k)
            fwd_step(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride, k);
 
        bwd_only(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    permute_inplace(const int n, const int* TDLS_RESTRICT piv, const int piv_stride,
                    T* TDLS_RESTRICT x, const int rhs_stride) noexcept {
        if (n <= 64) {
            // Bitmask cycle decomposition: one visited bit per row, the
            // whole state in a single 64-bit register.
            unsigned long long visited = 0ull;
            for (int s = 0; s < n; ++s) {
                if ((visited >> s) & 1ull) continue;
                const T tmp = TDLS_LUPP_DYN_X(s);
                int cur     = s;
                int nxt     = TDLS_LUPP_DYN_PIV(cur);
                while (nxt != s) {
                    TDLS_LUPP_DYN_X(cur) = TDLS_LUPP_DYN_X(nxt);
                    visited |= 1ull << cur;
                    cur = nxt;
                    nxt = TDLS_LUPP_DYN_PIV(cur);
                }
                TDLS_LUPP_DYN_X(cur) = tmp;
                visited |= 1ull << 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_DYN_PIV(s);
                while (probe > s)
                    probe = TDLS_LUPP_DYN_PIV(probe);
                if (probe != s) continue;
 
                const T tmp = TDLS_LUPP_DYN_X(s);
                int cur     = s;
                int nxt     = TDLS_LUPP_DYN_PIV(cur);
                while (nxt != s) {
                    TDLS_LUPP_DYN_X(cur) = TDLS_LUPP_DYN_X(nxt);
                    cur                  = nxt;
                    nxt                  = TDLS_LUPP_DYN_PIV(cur);
                }
                TDLS_LUPP_DYN_X(cur) = tmp;
            }
        }
    }
 
    /* =====================================================================
       TRANSPOSED SUBSTITUTION - column-major layout, as in the
       compile-time solver: A = U^T L^T P, so A x = b is U^T z = b, then
       L^T w = z, then x = P^T w, with the same tile kernels for the fused
       forward pass of solve_inplace and for the separate passes.
       ===================================================================== */
 
    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 nrhs, const int k0,
                                const int ke) noexcept {
        for (int kk = 0; kk < ke; ++kk) {
            for (int w = 0; w < nrhs; ++w)
                TDLS_LUPP_DYN_XW(w, k0 + kk) *= t[kk * tile_size + kk]; // diag holds 1/pivot
            for (int i = kk + 1; i < ke; ++i) {
                for (int w = 0; w < nrhs; ++w)
                    TDLS_LUPP_DYN_XW(w, k0 + i) -=
                        t[kk * tile_size + i] * TDLS_LUPP_DYN_XW(w, k0 + kk);
            }
        }
    }
 
    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 nrhs, const int k0, const int m0,
                        const int ke, const int me) noexcept {
        for (int i = 0; i < me; ++i) {
            for (int w = 0; w < nrhs; ++w) {
                T sum = T(0);
                for (int j = 0; j < ke; ++j)
                    sum += u[j * tile_size + i] * TDLS_LUPP_DYN_XW(w, k0 + j);
                TDLS_LUPP_DYN_XW(w, m0 + i) -= sum;
            }
        }
    }
 
    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 nrhs, const int k0, const int m0,
                        const int ke, const int me) noexcept {
        for (int i = 0; i < ke; ++i) {
            for (int w = 0; w < nrhs; ++w) {
                T sum = T(0);
                for (int j = 0; j < me; ++j)
                    sum += l[j * tile_size + i] * TDLS_LUPP_DYN_XW(w, m0 + j);
                TDLS_LUPP_DYN_XW(w, k0 + i) -= sum;
            }
        }
    }
 
    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 nrhs, const int k0, const int ke) noexcept {
        for (int kk = ke - 1; kk > 0; --kk) {
            for (int i = 0; i < kk; ++i) {
                for (int w = 0; w < nrhs; ++w)
                    TDLS_LUPP_DYN_XW(w, k0 + i) -=
                        l[kk * tile_size + i] * TDLS_LUPP_DYN_XW(w, k0 + kk);
            }
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_fwd_step(const int n, const int nrhs, 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;
        const int nt = num_tiles(n);
        const int ke = tile_size_at(k0, n);
        T Ukk[tile_size * tile_size];
        load_tile_piv_upper(n, A, A_stride, piv, piv_stride, k0, k0, Ukk, ke);
        transposed_solve_tile_upper(Ukk, x, rhs_stride, xcol_stride, nrhs, k0, ke);
        for (int m = k + 1; m < nt; ++m) {
            const int me = tile_size_at(m * tile_size, n);
            T Ukm[tile_size * tile_size];
            load_tile_piv(n, A, A_stride, piv, piv_stride, k0, m * tile_size, Ukm, ke, me);
            transposed_push_one(Ukm, x, rhs_stride, xcol_stride, nrhs, k0, m * tile_size, ke, me);
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_bwd_step(const int n, const int nrhs, 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;
        const int nt = num_tiles(n);
        const int ke = tile_size_at(k0, n);
        for (int m = k + 1; m < nt; ++m) {
            const int me = tile_size_at(m * tile_size, n);
            T Lmk[tile_size * tile_size];
            load_tile_piv(n, A, A_stride, piv, piv_stride, m * tile_size, k0, Lmk, me, ke);
            transposed_pull_one(Lmk, x, rhs_stride, xcol_stride, nrhs, k0, m * tile_size, ke, me);
        }
        T Lkk[tile_size * tile_size];
        load_tile_piv_lower(n, A, A_stride, piv, piv_stride, k0, k0, Lkk, ke);
        transposed_solve_tile_lower_unit(Lkk, x, rhs_stride, xcol_stride, nrhs, k0, ke);
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_bwd_only(const int n, const int nrhs, 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 = num_tiles(n) - 1; k >= 0; --k)
            transposed_bwd_step(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride,
                                k);
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_fwd_bwd(const int n, const int nrhs, 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 {
        const int nt = num_tiles(n);
        for (int k = 0; k < nt; ++k)
            transposed_fwd_step(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride,
                                k);
        transposed_bwd_only(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    unpermute_inplace(const int n, const int* TDLS_RESTRICT piv, const int piv_stride,
                      T* TDLS_RESTRICT x, const int rhs_stride) noexcept {
        if (n <= 64) {
            unsigned long long visited = 0ull;
            for (int s = 0; s < n; ++s) {
                if ((visited >> s) & 1ull) continue;
                T carry = TDLS_LUPP_DYN_X(s);
                int dst = TDLS_LUPP_DYN_PIV(s);
                while (dst != s) {
                    const T displaced    = TDLS_LUPP_DYN_X(dst);
                    TDLS_LUPP_DYN_X(dst) = carry;
                    carry                = displaced;
                    visited |= 1ull << dst;
                    dst = TDLS_LUPP_DYN_PIV(dst);
                }
                TDLS_LUPP_DYN_X(s) = carry;
                visited |= 1ull << s;
            }
        } else {
            for (int s = 0; s < n; ++s) {
                int probe = TDLS_LUPP_DYN_PIV(s);
                while (probe > s)
                    probe = TDLS_LUPP_DYN_PIV(probe);
                if (probe != s) continue;
 
                T carry = TDLS_LUPP_DYN_X(s);
                int dst = TDLS_LUPP_DYN_PIV(s);
                while (dst != s) {
                    const T displaced    = TDLS_LUPP_DYN_X(dst);
                    TDLS_LUPP_DYN_X(dst) = carry;
                    carry                = displaced;
                    dst                  = TDLS_LUPP_DYN_PIV(dst);
                }
                TDLS_LUPP_DYN_X(s) = carry;
            }
        }
    }
 
    /* =====================================================================
       SUBSTITUTION - public entry points.
       ===================================================================== */
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute(int n, 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 {
        n = detail::runtime_dimension(n);
        if constexpr (physical) {
            for (int i = 0; i < n; ++i)
                TDLS_LUPP_DYN_X(i) = TDLS_LUPP_DYN_B(i);
            physical_multirhs<0, 3>(n, 1, A, A_stride, piv, piv_stride, x, rhs_stride, 0);
            return;
        }
        if constexpr (transposed) {
            // Column-major layout: b consumed unpermuted, x = P^T w at the end.
            for (int i = 0; i < n; ++i)
                TDLS_LUPP_DYN_X(i) = TDLS_LUPP_DYN_B(i);
            transposed_fwd_bwd(n, 1, A, A_stride, piv, piv_stride, x, rhs_stride, 0);
            unpermute_inplace(n, piv, piv_stride, x, rhs_stride);
        } else {
            for (int i = 0; i < n; ++i)
                TDLS_LUPP_DYN_X(i) = TDLS_LUPP_DYN_B(TDLS_LUPP_DYN_PIV(i));
            fwd_bwd(n, 1, A, A_stride, piv, piv_stride, x, rhs_stride, 0);
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_canonical(int n, 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 {
        n = detail::runtime_dimension(n);
        if constexpr (physical) {
            for (int i = 0; i < n; ++i)
                TDLS_LUPP_DYN_X(i) = i == col ? T(1) : T(0);
            physical_multirhs<0, 3>(n, 1, A, A_stride, piv, piv_stride, x, rhs_stride, 0);
            return;
        }
        substitute_canonical_multirhs_pass(n, 1, A, A_stride, piv, piv_stride, col, x, rhs_stride,
                                           0);
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void substitute_canonical_multirhs_pass(
        const int n, const int nrhs, 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, x = P^T w at the end.
            for (int i = 0; i < n; ++i)
                for (int w = 0; w < nrhs; ++w)
                    TDLS_LUPP_DYN_XW(w, i) = (i == col0 + w) ? T(1) : T(0);
            transposed_fwd_bwd(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
            for (int w = 0; w < nrhs; ++w)
                unpermute_inplace(n, piv, piv_stride, x + TDLS_STRIDED_WIDE(w, xcol_stride),
                                  rhs_stride);
        } else {
            for (int i = 0; i < n; ++i) {
                const int p = TDLS_LUPP_DYN_PIV(i);
                for (int w = 0; w < nrhs; ++w)
                    TDLS_LUPP_DYN_XW(w, i) = (p == col0 + w) ? T(1) : T(0);
            }
            fwd_bwd(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
        }
    }
 
    template<int pass_width = 0>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_canonical_multirhs(int n, int nrhs, 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 {
        n    = detail::runtime_dimension(n);
        nrhs = detail::runtime_dimension(nrhs);
        if constexpr (physical) {
            for (int c = 0; c < nrhs; ++c)
                for (int i = 0; i < n; ++i)
                    TDLS_LUPP_DYN_XW(c, i) = i == col0 + c ? T(1) : T(0);
            physical_multirhs<pass_width, 3>(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride,
                                             xcol_stride);
            return;
        }
        static_assert(pass_width >= 0, "tdls: pass_width must not be negative");
        if constexpr (pass_width == 0) {
            substitute_canonical_multirhs_pass(n, nrhs, A, A_stride, piv, piv_stride, col0, x,
                                               rhs_stride, xcol_stride);
        } else {
            // A pass_width of nrhs or more runs the loop once, the whole
            // block being its single pass.
            for (int c0 = 0; c0 < nrhs; c0 += pass_width) {
                const int cw          = (pass_width > nrhs - c0) ? nrhs - c0 : pass_width;
                const std::size_t off = TDLS_STRIDED_WIDE(c0, xcol_stride);
                substitute_canonical_multirhs_pass(n, cw, A, A_stride, piv, piv_stride, col0 + c0,
                                                   x + off, rhs_stride, xcol_stride);
            }
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_multirhs_pass(const int n, const int nrhs, 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.
            for (int i = 0; i < n; ++i)
                for (int w = 0; w < nrhs; ++w)
                    TDLS_LUPP_DYN_XW(w, i) = TDLS_LUPP_DYN_BW(w, i);
            transposed_fwd_bwd(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
            for (int w = 0; w < nrhs; ++w)
                unpermute_inplace(n, piv, piv_stride, x + TDLS_STRIDED_WIDE(w, xcol_stride),
                                  rhs_stride);
        } else {
            for (int i = 0; i < n; ++i) {
                const int p = TDLS_LUPP_DYN_PIV(i);
                for (int w = 0; w < nrhs; ++w)
                    TDLS_LUPP_DYN_XW(w, i) = TDLS_LUPP_DYN_BW(w, p);
            }
            fwd_bwd(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
        }
    }
 
    template<int pass_width = 0>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_multirhs(int n, int nrhs, 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 {
        n    = detail::runtime_dimension(n);
        nrhs = detail::runtime_dimension(nrhs);
        if constexpr (physical) {
            for (int c = 0; c < nrhs; ++c)
                for (int i = 0; i < n; ++i)
                    TDLS_LUPP_DYN_XW(c, i) = TDLS_LUPP_DYN_BW(c, i);
            physical_multirhs<pass_width, 3>(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride,
                                             xcol_stride);
            return;
        }
        static_assert(pass_width >= 0, "tdls: pass_width must not be negative");
        if constexpr (pass_width == 0) {
            substitute_multirhs_pass(n, nrhs, A, A_stride, piv, piv_stride, b, x, rhs_stride,
                                     xcol_stride);
        } else {
            // A pass_width of nrhs or more runs the loop once, the whole
            // block being its single pass.
            for (int c0 = 0; c0 < nrhs; c0 += pass_width) {
                const int cw          = (pass_width > nrhs - c0) ? nrhs - c0 : pass_width;
                const std::size_t off = TDLS_STRIDED_WIDE(c0, xcol_stride);
                substitute_multirhs_pass(n, cw, A, A_stride, piv, piv_stride, b + off, x + off,
                                         rhs_stride, xcol_stride);
            }
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_inplace(int n, 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 {
        n = detail::runtime_dimension(n);
        if constexpr (physical) {
            physical_multirhs<0, 3>(n, 1, 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(n, 1, A, A_stride, piv, piv_stride, x, rhs_stride, 0);
            unpermute_inplace(n, piv, piv_stride, x, rhs_stride);
        } else {
            permute_inplace(n, piv, piv_stride, x, rhs_stride);
            fwd_bwd(n, 1, A, A_stride, piv, piv_stride, x, rhs_stride, 0);
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_inplace_multirhs_pass(const int n, const int nrhs, 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(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
            for (int w = 0; w < nrhs; ++w)
                unpermute_inplace(n, piv, piv_stride, x + TDLS_STRIDED_WIDE(w, xcol_stride),
                                  rhs_stride);
        } else {
            for (int w = 0; w < nrhs; ++w)
                permute_inplace(n, piv, piv_stride, x + TDLS_STRIDED_WIDE(w, xcol_stride),
                                rhs_stride);
            fwd_bwd(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    forward_inplace_multirhs_pass(const int n, const int nrhs, 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 {
        const int nt = num_tiles(n);
        if constexpr (transposed) {
            for (int k = 0; k < nt; ++k)
                transposed_fwd_step(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride,
                                    xcol_stride, k);
        } else {
            for (int w = 0; w < nrhs; ++w)
                permute_inplace(n, piv, piv_stride, x + TDLS_STRIDED_WIDE(w, xcol_stride),
                                rhs_stride);
            for (int k = 0; k < nt; ++k)
                fwd_step(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride, k);
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    backward_inplace_multirhs_pass(const int n, const int nrhs, 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(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
            for (int w = 0; w < nrhs; ++w)
                unpermute_inplace(n, piv, piv_stride, x + TDLS_STRIDED_WIDE(w, xcol_stride),
                                  rhs_stride);
        } else {
            bwd_only(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
        }
    }
 
    template<int pass_width = 0>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_inplace_multirhs(int n, int nrhs, 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 {
        n    = detail::runtime_dimension(n);
        nrhs = detail::runtime_dimension(nrhs);
        if constexpr (physical) {
            physical_multirhs<pass_width, 3>(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride,
                                             xcol_stride);
            return;
        }
        static_assert(pass_width >= 0, "tdls: pass_width must not be negative");
        if constexpr (pass_width == 0) {
            substitute_inplace_multirhs_pass(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride,
                                             xcol_stride);
        } else {
            // A pass_width of nrhs or more runs the loop once, the whole
            // block being its single pass.
            for (int c0 = 0; c0 < nrhs; c0 += pass_width) {
                const int cw          = (pass_width > nrhs - c0) ? nrhs - c0 : pass_width;
                const std::size_t off = TDLS_STRIDED_WIDE(c0, xcol_stride);
                substitute_inplace_multirhs_pass(n, cw, A, A_stride, piv, piv_stride, x + off,
                                                 rhs_stride, xcol_stride);
            }
        }
    }
 
    template<int pass_width = 0>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    forward_inplace_multirhs(int n, int nrhs, 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 {
        n    = detail::runtime_dimension(n);
        nrhs = detail::runtime_dimension(nrhs);
        if constexpr (physical) {
            physical_multirhs<pass_width, 1>(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride,
                                             xcol_stride);
            return;
        }
        static_assert(pass_width >= 0, "tdls: pass_width must not be negative");
        if constexpr (pass_width == 0) {
            forward_inplace_multirhs_pass(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride,
                                          xcol_stride);
        } else {
            for (int c0 = 0; c0 < nrhs; c0 += pass_width) {
                const int cw          = (pass_width > nrhs - c0) ? nrhs - c0 : pass_width;
                const std::size_t off = TDLS_STRIDED_WIDE(c0, xcol_stride);
                forward_inplace_multirhs_pass(n, cw, A, A_stride, piv, piv_stride, x + off,
                                              rhs_stride, xcol_stride);
            }
        }
    }
 
    template<int pass_width = 0>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    backward_inplace_multirhs(int n, int nrhs, 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 {
        n    = detail::runtime_dimension(n);
        nrhs = detail::runtime_dimension(nrhs);
        if constexpr (physical) {
            physical_multirhs<pass_width, 2>(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride,
                                             xcol_stride);
            return;
        }
        static_assert(pass_width >= 0, "tdls: pass_width must not be negative");
        if constexpr (pass_width == 0) {
            backward_inplace_multirhs_pass(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride,
                                           xcol_stride);
        } else {
            for (int c0 = 0; c0 < nrhs; c0 += pass_width) {
                const int cw          = (pass_width > nrhs - c0) ? nrhs - c0 : pass_width;
                const std::size_t off = TDLS_STRIDED_WIDE(c0, xcol_stride);
                backward_inplace_multirhs_pass(n, cw, A, A_stride, piv, piv_stride, x + off,
                                               rhs_stride, xcol_stride);
            }
        }
    }
 
    /* =====================================================================
       SOLVE - convenience wrappers.
       ===================================================================== */
 
    template<bool oot_diagnostics = true>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve(int n, 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 {
        n = detail::runtime_dimension(n);
        if (!factorize<oot_diagnostics>(n, A, A_stride, piv, piv_stride, oot_count)) return false;
        substitute(n, A, A_stride, piv, piv_stride, b, x, rhs_stride);
        return true;
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve(int n, 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 {
        n          = detail::runtime_dimension(n);
        int unused = 0;
        return solve<false>(n, A, A_stride, piv, piv_stride, b, x, rhs_stride, unused);
    }
 
    template<int pass_width = 0, bool oot_diagnostics = true>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve_multirhs(int n, int nrhs, 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 {
        n    = detail::runtime_dimension(n);
        nrhs = detail::runtime_dimension(nrhs);
        if (!factorize<oot_diagnostics>(n, A, A_stride, piv, piv_stride, oot_count)) return false;
        substitute_multirhs<pass_width>(n, nrhs, A, A_stride, piv, piv_stride, b, x, rhs_stride,
                                        xcol_stride);
        return true;
    }
 
    template<int pass_width = 0>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve_multirhs(int n, int nrhs, 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 {
        n          = detail::runtime_dimension(n);
        nrhs       = detail::runtime_dimension(nrhs);
        int unused = 0;
        return solve_multirhs<pass_width, false>(n, nrhs, A, A_stride, piv, piv_stride, b, x,
                                                 rhs_stride, xcol_stride, unused);
    }
 
    template<bool oot_diagnostics = true>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve_inplace(int n, 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 {
        n = detail::runtime_dimension(n);
        if constexpr (physical) {
            int weak      = 0;
            const bool ok = physical_factorize(n, A, A_stride, piv, piv_stride, weak);
            if constexpr (oot_diagnostics) oot_count = weak;
            if (!ok) return false;
            physical_multirhs<0, 3>(n, 1, A, A_stride, piv, piv_stride, y, rhs_stride, 0);
            return true;
        }
        if (!factorize<oot_diagnostics, true>(n, 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(n, 1, A, A_stride, piv, piv_stride, y, rhs_stride, 0);
            unpermute_inplace(n, piv, piv_stride, y, rhs_stride);
        } else {
            bwd_only(n, 1, A, A_stride, piv, piv_stride, y, rhs_stride, 0);
        }
        return true;
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve_inplace(int n, 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 {
        n          = detail::runtime_dimension(n);
        int unused = 0;
        return solve_inplace<false>(n, A, A_stride, piv, piv_stride, y, rhs_stride, unused);
    }
 
    template<int pass_width = 0, bool oot_diagnostics = true>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve_inplace_multirhs(int n, int nrhs, 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 {
        n    = detail::runtime_dimension(n);
        nrhs = detail::runtime_dimension(nrhs);
        if (!factorize<oot_diagnostics>(n, A, A_stride, piv, piv_stride, oot_count)) return false;
        substitute_inplace_multirhs<pass_width>(n, nrhs, A, A_stride, piv, piv_stride, y,
                                                rhs_stride, xcol_stride);
        return true;
    }
 
    template<int pass_width = 0>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve_inplace_multirhs(int n, int nrhs, 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 {
        n          = detail::runtime_dimension(n);
        nrhs       = detail::runtime_dimension(nrhs);
        int unused = 0;
        return solve_inplace_multirhs<pass_width, false>(n, nrhs, A, A_stride, piv, piv_stride, y,
                                                         rhs_stride, xcol_stride, unused);
    }
};
 
 
 
#undef TDLS_LUPP_DYN_A
#undef TDLS_LUPP_DYN_PIV
#undef TDLS_LUPP_DYN_X
#undef TDLS_LUPP_DYN_B
#undef TDLS_LUPP_DYN_XW
#undef TDLS_LUPP_DYN_BW
#undef TDLS_LUPP_DYN_Y
 
 
 
} // namespace tdls
 
#if defined(TDLS_DETAIL_COMPILER_GCC)
#pragma GCC diagnostic pop
#endif
 
 
 
#endif // TDLS_SOLVERS_TILED_LUPP_SOLVER_DYNAMIC_HPP