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 <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(__GNUC__) && !defined(__clang__)
#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_LAYOUT_INDEX(r, c, n) * unsigned(A_stride)]
#define TDLS_LUPP_DYN_PIV(i) piv[unsigned(i) * unsigned(piv_stride)]
#define TDLS_LUPP_DYN_X(i) x[unsigned(i) * unsigned(rhs_stride)]
#define TDLS_LUPP_DYN_B(i) b[unsigned(i) * unsigned(rhs_stride)]
#define TDLS_LUPP_DYN_XW(w, i)                                                                     \
    x[unsigned(i) * unsigned(rhs_stride) + unsigned(w) * unsigned(xcol_stride)]
#define TDLS_LUPP_DYN_BW(w, i)                                                                     \
    b[unsigned(i) * unsigned(rhs_stride) + unsigned(w) * unsigned(xcol_stride)]
#define TDLS_LUPP_DYN_Y(i) y[unsigned(i) * unsigned(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_threshold = Config.oot_threshold;
    static constexpr T singular_floor = Config.singular_floor;
 
    static_assert(
        Config.oot_threshold.is_finite() && Config.singular_floor.is_finite(),
        "TiledLUppSolverDynamic: oot_threshold and singular_floor must be finite (and fit a "
        "63-bit mantissa)");
    static_assert(singular_floor > T(0), "TiledLUppSolverDynamic: singular_floor must be positive");
    static_assert(singular_floor <= oot_threshold,
                  "TiledLUppSolverDynamic: singular_floor must not exceed oot_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_FORCEINLINE static constexpr int num_tiles(const int n) noexcept {
        return (n + tile_size - 1) / tile_size;
    }
 
    TDLS_HOST_DEVICE TDLS_FORCEINLINE 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_FORCEINLINE 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_FORCEINLINE 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_FORCEINLINE 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_FORCEINLINE 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_FORCEINLINE 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_FORCEINLINE 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_FORCEINLINE 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_FORCEINLINE 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_FORCEINLINE 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_FORCEINLINE 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_FORCEINLINE 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_FORCEINLINE 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;
            }
        }
 
        int piv_row; // winning global (logical) row
 
        if (best >= oot_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_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) {
                // The fused RHS follows the rows through pivoting.
                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_FORCEINLINE 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_FORCEINLINE 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_FORCEINLINE 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_FORCEINLINE 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.
        if constexpr (fuse_rhs) {
            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_FORCEINLINE 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: unit-lower-solve it while the tile is in registers.
        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
        for (int tj = k + 1; tj < nt; ++tj)
            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_FORCEINLINE 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_FORCEINLINE 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.
        if constexpr (fuse_rhs) {
            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_FORCEINLINE 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);
    }
 
    template<bool oot_diagnostics, bool fuse_rhs>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE 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: unit-lower-solve it while the tile is in registers.
        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
        for (int tj = k + 1; tj < nt; ++tj)
            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;
    }
 
    /* =====================================================================
       FACTORIZE - public entry point.
       ===================================================================== */
 
    template<bool oot_diagnostics = true, bool fuse_rhs = false>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr bool
    factorize(const 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 {
 
        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_FORCEINLINE static constexpr bool
    factorize(const int n, T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
              const int piv_stride) noexcept {
        int unused = 0;
        return factorize<false>(n, A, A_stride, piv, piv_stride, unused);
    }
 
    /* =====================================================================
       SUBSTITUTION - schedule-independent.
       ===================================================================== */
 
    TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 int 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_FORCEINLINE 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 int 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_FORCEINLINE 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 int 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_FORCEINLINE 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 int 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_FORCEINLINE 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 int 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_FORCEINLINE 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 int 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_FORCEINLINE 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;
            }
        }
    }
 
    /* =====================================================================
       SUBSTITUTION - public entry points.
       ===================================================================== */
 
    TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr void
    substitute(const 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 {
        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_FORCEINLINE static constexpr void
    substitute_canonical(const 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 {
        substitute_canonical_multirhs_pass(n, 1, A, A_stride, piv, piv_stride, col, x, rhs_stride,
                                           0);
    }
 
    TDLS_HOST_DEVICE TDLS_FORCEINLINE 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 int xcol_stride) noexcept {
        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_FORCEINLINE static constexpr void
    substitute_canonical_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, const int col0, T* TDLS_RESTRICT x,
                                  const int rhs_stride, const int xcol_stride) noexcept {
        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 unsigned off = unsigned(c0) * unsigned(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_FORCEINLINE 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 int xcol_stride) noexcept {
        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_FORCEINLINE static constexpr void
    substitute_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,
                        const T* TDLS_RESTRICT b, T* TDLS_RESTRICT x, const int rhs_stride,
                        const int xcol_stride) noexcept {
        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 unsigned off = unsigned(c0) * unsigned(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_FORCEINLINE static constexpr void
    substitute_inplace(const 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 {
        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_FORCEINLINE 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 int xcol_stride) noexcept {
        for (int w = 0; w < nrhs; ++w)
            permute_inplace(n, piv, piv_stride, x + unsigned(w) * unsigned(xcol_stride),
                            rhs_stride);
        fwd_bwd(n, nrhs, A, A_stride, piv, piv_stride, x, rhs_stride, xcol_stride);
    }
 
    template<int pass_width = 0>
    TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr void
    substitute_inplace_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 int xcol_stride) noexcept {
        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 unsigned off = unsigned(c0) * unsigned(xcol_stride);
                substitute_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_FORCEINLINE static constexpr bool
    solve(const 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 {
        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_FORCEINLINE static constexpr bool
    solve(const 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 {
        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_FORCEINLINE static constexpr bool
    solve_multirhs(const int n, const 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 int xcol_stride,
                   int& oot_count) noexcept {
        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_FORCEINLINE static constexpr bool
    solve_multirhs(const int n, const 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 int xcol_stride) noexcept {
        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_FORCEINLINE static constexpr bool
    solve_inplace(const 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 {
        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.
        bwd_only(n, 1, A, A_stride, piv, piv_stride, y, rhs_stride, 0);
        return true;
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr bool
    solve_inplace(const 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 {
        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_FORCEINLINE static constexpr bool
    solve_inplace_multirhs(const int n, const 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 int xcol_stride, int& oot_count) noexcept {
        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_FORCEINLINE static constexpr bool
    solve_inplace_multirhs(const int n, const 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 int xcol_stride) noexcept {
        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(__GNUC__) && !defined(__clang__)
#pragma GCC diagnostic pop
#endif
 
 
 
#endif // TDLS_SOLVERS_TILED_LUPP_SOLVER_DYNAMIC_HPP