solver_static.hpp

Go to the documentation of this file.

#ifndef TDLS_SOLVERS_COOPERATIVE_LUPP_SOLVER_STATIC_HPP
#define TDLS_SOLVERS_COOPERATIVE_LUPP_SOLVER_STATIC_HPP
 
 
 
/* The elimination and substitution kernels of this file are derived from
   MAGMA 2.10.0, magmablas/zgesv_batched_small.cu, function
   zgesv_batched_small_device (authors Azzam Haidar and Ahmad Abdelfattah).
   Those portions are distributed under the MAGMA license, reproduced
   below as required by it:
 
   Copyright (c) 2009-2023, The University of Tennessee
   All rights reserved.
 
     Redistribution and use in source and binary forms, with or without
     modification, are permitted provided that the following conditions
     are met:
 
     * Redistributions of source code must retain the above copyright
       notice, this list of conditions and the following disclaimer.
     * Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
     * Neither the name of the University of Tennessee, Knoxville nor the
       names of its contributors may be used to endorse or promote products
       derived from this software without specific prior written permission.
 
     This software is provided by the copyright holders and contributors
     ``as is'' and any express or implied warranties, including, but not
     limited to, the implied warranties of merchantability and fitness for
     a particular purpose are disclaimed. In no event shall the copyright
     holders or contributors be liable for any direct, indirect, incidental,
     special, exemplary, or consequential damages (including, but not
     limited to, procurement of substitute goods or services; loss of use,
     data, or profits; or business interruption) however caused and on any
     theory of liability, whether in contract, strict liability, or tort
     (including negligence or otherwise) arising in any way out of the use
     of this software, even if advised of the possibility of such damage.
*/
 
 
 
 
 
 
#include <type_traits>
#include <utility>
 
#include <tdls/core/group.hpp>
#include <tdls/core/macros.hpp>
#include <tdls/core/math.hpp>
#include <tdls/core/target.hpp>
#include <tdls/solvers/cooperative_lupp/config.hpp>
 
 
 
// Older gcc releases cannot prove that the rows of a thread are loaded
// before being read: the load of a phantom-capable row slot is guarded by
// the row index and its later uses by rowid, which starts equal to it and
// stays at or above N for a phantom row. The resulting
// -Wmaybe-uninitialized reports are spurious. Under UBSan instrumentation
// at -O2 and above, gcc additionally loses the value ranges of the fully
// unrolled loop indices and reports impossible subscripts through
// -Warray-bounds, exactly as in the TiledLUpp headers; 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.
#if defined(TDLS_DETAIL_COMPILER_GCC)
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
#pragma GCC diagnostic ignored "-Warray-bounds"
#endif
 
 
 
// clang reports a forced unrolling that the optimizer could not perform
// through -Wpass-failed. The unrolling requested by TDLS_UNROLL_FORCE is
// a performance hint: a failed hint does not affect correctness. The
// suppression is scoped to this header and to that warning only.
#if defined(TDLS_DETAIL_COMPILER_CLANG)
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wpass-failed"
#endif
 
 
 
// nvc++ reports through no_device_stack, under the managed memory of
// -stdpar=gpu, every lambda that captures local objects by reference:
// such a lambda handed to a parallel algorithm would read the host stack.
// The lambdas of this header are called on the spot, on the stack of the
// thread that runs the solver: the report is a false positive. The
// suppression is scoped to this header and to that warning only.
#if defined(TDLS_DETAIL_COMPILER_NVCXX)
#pragma diag_push
#pragma diag_suppress no_device_stack
#endif
 
namespace tdls {
 
 
 
namespace detail {
 
template<typename Fn, int... Ks>
TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU constexpr void
unroll_K_impl(Fn& fn, std::integer_sequence<int, Ks...>) noexcept {
    (fn(std::integral_constant<int, Ks>{}), ...);
}
 
template<int R, typename Fn>
TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU constexpr void unroll_K(Fn&& fn) noexcept {
    unroll_K_impl(fn, std::make_integer_sequence<int, R>{});
}
 
} // namespace detail
 
 
 
/* Addressing macros. They expand inside member functions where the
   parameter names (tx, A, A_stride, piv, piv_stride, b, x, y,
   rhs_stride), the residency template booleans and the class constants
   are in scope. K is a row slot of the calling thread, its row is
   tx + K * threads_per_system. #undef'd at the end of this header. */
 
#define TDLS_COOP_LUPP_A(K, c)                                                                     \
    A[internal_matrix ? unsigned((c) * rows_per_thread + (K))                                      \
                      : TDLS_STRIDED((c) * N + tx + (K) * threads_per_system, A_stride)]
#define TDLS_COOP_LUPP_PIV(K)                                                                      \
    piv[internal_piv ? unsigned(K) : TDLS_STRIDED(tx + (K) * threads_per_system, piv_stride)]
#define TDLS_COOP_LUPP_B(K)                                                                        \
    b[internal_rhs ? unsigned(K) : TDLS_STRIDED(tx + (K) * threads_per_system, rhs_stride)]
#define TDLS_COOP_LUPP_X(K)                                                                        \
    x[internal_rhs ? unsigned(K) : TDLS_STRIDED(tx + (K) * threads_per_system, rhs_stride)]
#define TDLS_COOP_LUPP_Y(K)                                                                        \
    y[internal_rhs ? unsigned(K) : TDLS_STRIDED(tx + (K) * threads_per_system, rhs_stride)]
 
 
 
template<typename T, int N, CooperativeLUppConfig<T> Config = CooperativeLUppConfig<T>{}>
struct CooperativeLUppSolverStatic {
 
    static_assert(N >= 1, "CooperativeLUppSolverStatic: N must be >= 1");
    static_assert(
        Config.row_interchange != RowInterchange::Auto,
        "CooperativeLUppSolverStatic: row_interchange must be Logical or Physical (Auto is the "
        "default of the blocks of ArrowheadLUpp only)");
    static_assert(Config.rows_per_thread >= 1,
                  "CooperativeLUppSolverStatic: rows_per_thread must be >= 1");
    static_assert(Config.relative_pivot_threshold.is_finite() && Config.singular_floor.is_finite(),
                  "CooperativeLUppSolverStatic: relative_pivot_threshold and singular_floor must "
                  "be finite (and fit a 63-bit mantissa)");
 
    static constexpr int rows_per_thread = Config.rows_per_thread < 1   ? 1
                                           : Config.rows_per_thread < N ? Config.rows_per_thread
                                                                        : N;
    static constexpr int threads_per_system = (N + rows_per_thread - 1) / rows_per_thread;
    static constexpr int workspace_size =
        detail::barrier_elements + (Config.row_interchange == RowInterchange::Logical
                                        ? 3 * N
                                        : 2 * N + 2 * threads_per_system + 5);
    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::RowMajor;
 
    static_assert(relative_pivot_threshold > T(0) && relative_pivot_threshold <= T(1),
                  "CooperativeLUppSolverStatic: relative_pivot_threshold must be in (0, 1]");
    static_assert(singular_floor > T(0),
                  "CooperativeLUppSolverStatic: singular_floor must be positive");
 
    /* =====================================================================
       Row slots.
       Slot K of thread tx holds row tx + K * threads_per_system. When N is
       not a multiple of rows_per_thread, the last slots hold phantom rows
       (index >= N) for some or all threads. The classification is a
       compile-time function of K: a slot without any real row is skipped
       entirely, a slot full for every thread needs no guard, and only a
       mixed slot tests its row at run time.
       ===================================================================== */
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    slot_has_rows(const int K) noexcept {
        return K * threads_per_system < N;
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    slot_is_full(const int K) noexcept {
        return (K + 1) * threads_per_system <= N;
    }
 
    template<typename Fn>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    at_position(const int tx, const int pos, Fn&& fn) noexcept {
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (K == pos / threads_per_system && tx == pos % threads_per_system) fn(K_tag);
            }
        });
    }
 
    /* =====================================================================
       Operand movement.
       Each thread reads and writes only the entries of its own rows, so
       the operands themselves never carry data between threads: the
       workspace does.
       ===================================================================== */
 
    template<bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_rows(const int tx, const T* TDLS_RESTRICT A, const int A_stride,
              T (&rA)[rows_per_thread][N]) noexcept {
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (slot_is_full(K) || tx + K * threads_per_system < N) {
                    if constexpr (Config.unroll_loops) {
                        TDLS_UNROLL_FORCE
                        for (int c = 0; c < N; ++c)
                            rA[K][c] = TDLS_COOP_LUPP_A(K, c);
                    } else {
                        for (int c = 0; c < N; ++c)
                            rA[K][c] = TDLS_COOP_LUPP_A(K, c);
                    }
                }
            }
        });
    }
 
    template<bool internal_matrix>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    store_rows(const int tx, T* TDLS_RESTRICT A, const int A_stride,
               const T (&rA)[rows_per_thread][N]) noexcept {
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (slot_is_full(K) || tx + K * threads_per_system < N) {
                    if constexpr (Config.unroll_loops) {
                        TDLS_UNROLL_FORCE
                        for (int c = 0; c < N; ++c)
                            TDLS_COOP_LUPP_A(K, c) = rA[K][c];
                    } else {
                        for (int c = 0; c < N; ++c)
                            TDLS_COOP_LUPP_A(K, c) = rA[K][c];
                    }
                }
            }
        });
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    init_piv(const int tx, int (&rpiv)[rows_per_thread]) noexcept {
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            rpiv[K]         = tx + K * threads_per_system;
        });
    }
 
    template<bool internal_piv>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_piv(const int tx, const int* TDLS_RESTRICT piv, const int piv_stride,
             int (&rpiv)[rows_per_thread]) noexcept {
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if (slot_is_full(K) || tx + K * threads_per_system < N)
                rpiv[K] = TDLS_COOP_LUPP_PIV(K);
            else
                rpiv[K] = tx + K * threads_per_system;
        });
    }
 
    template<bool internal_piv>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    store_piv(const int tx, int* TDLS_RESTRICT piv, const int piv_stride,
              const int (&rpiv)[rows_per_thread]) noexcept {
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (slot_is_full(K) || tx + K * threads_per_system < N)
                    TDLS_COOP_LUPP_PIV(K) = rpiv[K];
            }
        });
    }
 
    template<bool internal_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_rhs(const int tx, const T* TDLS_RESTRICT b, const int rhs_stride,
             T (&rB)[rows_per_thread]) noexcept {
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (slot_is_full(K) || tx + K * threads_per_system < N) rB[K] = TDLS_COOP_LUPP_B(K);
            }
        });
    }
 
    template<bool internal_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    store_solution(const int tx, T* TDLS_RESTRICT x, const int rhs_stride,
                   const T (&rB)[rows_per_thread]) noexcept {
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (slot_is_full(K) || tx + K * threads_per_system < N) TDLS_COOP_LUPP_X(K) = rB[K];
            }
        });
    }
 
    /* =====================================================================
       FACTORIZATION, right-looking, one column per step, as in MAGMA.
       Logical row interchanges, three barrier intervals per column: the
       magnitudes of the column are published at the positions of their
       rows, every thread finds the same pivot, the owner of the pivot row
       publishes it, and every thread eliminates its own rows below the
       pivot. Physical row interchanges, two barrier intervals per column:
       each thread publishes its best candidate, every thread finds the
       same pivot, the owner of position i publishes its row, the owner of
       the pivot row as well when the pivot is not in place, and every
       thread eliminates its rows below the pivot after the exchange. The
       pivot, the operations and their order are the same in both schemes.
       ===================================================================== */
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    keeps_pivot(const T in_place, const T column_max) noexcept {
        const T bound = relative_pivot_threshold * column_max;
        return in_place >= (bound > singular_floor ? bound : singular_floor);
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<bool fuse_rhs, typename Sync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    logical_elimination_step([[maybe_unused]] const int tx, const int i,
                             T (&rA)[rows_per_thread][N], [[maybe_unused]] T (&rB)[rows_per_thread],
                             int (&rowid)[rows_per_thread], T* sB, T* sx, T* dsx, int& linfo,
                             Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        // izamax: the magnitudes of column i, at the positions of their rows.
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (slot_is_full(K) || rowid[K] < N) dsx[rowid[K]] = detail::abs(rA[K][i]);
            }
        });
        sync();
 
        // Every thread scans the same magnitudes and finds the same pivot.
        // The scan reads the workspace, not the registers: no pragma.
        T rx_abs_max = dsx[i];
        int max_id   = i;
        for (int j = i + 1; j < N; j++) {
            if (dsx[j] > rx_abs_max) {
                max_id     = j;
                rx_abs_max = dsx[j];
            }
        }
        // Below the default threshold 1, the row in place keeps the pivot
        // when it is large enough.
        if constexpr (relative_pivot_threshold < T(1)) {
            if (keeps_pivot(dsx[i], rx_abs_max)) max_id = i;
        }
        const bool zero_pivot = (rx_abs_max < singular_floor);
        linfo                 = (zero_pivot && linfo == 0) ? (i + 1) : linfo;
        const T update        = zero_pivot ? T(0) : T(1);
 
        // The owner of the pivot row publishes it; positions i and max_id
        // are swapped. A phantom row never matches: its position is >= N.
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (rowid[K] == max_id) {
                    rowid[K] = i;
                    if constexpr (Config.unroll_loops) {
                        TDLS_UNROLL_FORCE
                        for (int j = i; j < N; j++)
                            sx[j] = update * rA[K][j];
                    } else {
                        for (int j = i; j < N; j++)
                            sx[j] = update * rA[K][j];
                    }
                    if constexpr (fuse_rhs && !transposed) sB[0] = rB[K];
                } else if (rowid[K] == i) {
                    rowid[K] = max_id;
                }
            }
        });
        if constexpr (fuse_rhs && transposed)
            at_position(tx, i, [&](auto K_tag) { sB[0] = rB[decltype(K_tag)::value]; });
        sync();
 
        const T reg            = zero_pivot ? T(1) : T(1) / sx[i];
        [[maybe_unused]] T z_i = T(0);
        if constexpr (fuse_rhs && transposed) z_i = sB[0] * reg;
 
        // scal and ger: every thread eliminates its rows below the pivot.
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if ((slot_is_full(K) || rowid[K] < N) && rowid[K] > i) {
                    rA[K][i] *= reg;
                    if constexpr (Config.unroll_loops) {
                        TDLS_UNROLL_FORCE
                        for (int j = i + 1; j < N; j++)
                            rA[K][j] -= rA[K][i] * sx[j];
                    } else {
                        for (int j = i + 1; j < N; j++)
                            rA[K][j] -= rA[K][i] * sx[j];
                    }
                    if constexpr (fuse_rhs && !transposed) rB[K] -= rA[K][i] * sB[0];
                }
                if constexpr (fuse_rhs && transposed) {
                    const int c = tx + K * threads_per_system;
                    if (c == i) rB[K] = z_i;
                    if ((slot_is_full(K) || c < N) && c > i) rB[K] -= sx[c] * z_i;
                }
            }
        });
        sync();
    }
 
    template<bool fuse_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    save_row(const T (&row)[N], [[maybe_unused]] const T& rhs, const int origin,
             T* record) noexcept {
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int c = 0; c < N; c++)
                record[c] = row[c];
        } else {
            for (int c = 0; c < N; c++)
                record[c] = row[c];
        }
        if constexpr (fuse_rhs) record[N] = rhs;
        record[N + 1] = static_cast<T>(origin);
    }
 
    template<bool fuse_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_row(T (&row)[N], [[maybe_unused]] T& rhs, int& origin, const T* record) noexcept {
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int c = 0; c < N; c++)
                row[c] = record[c];
        } else {
            for (int c = 0; c < N; c++)
                row[c] = record[c];
        }
        if constexpr (fuse_rhs) rhs = record[N];
        origin = static_cast<int>(record[N + 1]);
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool below(const int tx, const int K,
                                                                       const int i) noexcept {
        return (K + 1) * threads_per_system > i + 1 &&
               (slot_is_full(K) || tx + K * threads_per_system < N);
    }
 
    template<bool fuse_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    gather_row(const int tx, const int i, const int pos, const T (&rA)[rows_per_thread][N],
               [[maybe_unused]] const T (&rB)[rows_per_thread],
               const int (&origin)[rows_per_thread], T* record) noexcept {
        const int slot   = pos / threads_per_system;
        const auto entry = [&](const int c) {
            T v = T(0);
            detail::unroll_K<rows_per_thread>([&](auto K_tag) {
                constexpr int K = decltype(K_tag)::value;
                if constexpr (slot_has_rows(K)) {
                    if (below(tx, K, i)) v = K == slot ? rA[K][c] : v;
                }
            });
            return v;
        };
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int c = 0; c < N; c++)
                record[c] = entry(c);
        } else {
            for (int c = 0; c < N; c++)
                record[c] = entry(c);
        }
        T rhs     = T(0);
        int index = 0;
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (below(tx, K, i)) {
                    if constexpr (fuse_rhs) rhs = K == slot ? rB[K] : rhs;
                    index = K == slot ? origin[K] : index;
                }
            }
        });
        if constexpr (fuse_rhs) record[N] = rhs;
        record[N + 1] = static_cast<T>(index);
    }
 
    template<bool fuse_rhs>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    scatter_row(const int tx, const int i, const int pos, T (&rA)[rows_per_thread][N],
                [[maybe_unused]] T (&rB)[rows_per_thread], int (&origin)[rows_per_thread],
                const T* record) noexcept {
        const int slot   = pos / threads_per_system;
        const auto entry = [&](const int c) {
            const T v = record[c];
            detail::unroll_K<rows_per_thread>([&](auto K_tag) {
                constexpr int K = decltype(K_tag)::value;
                if constexpr (slot_has_rows(K)) {
                    if (below(tx, K, i)) rA[K][c] = K == slot ? v : rA[K][c];
                }
            });
        };
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int c = 0; c < N; c++)
                entry(c);
        } else {
            for (int c = 0; c < N; c++)
                entry(c);
        }
        [[maybe_unused]] const T rhs = fuse_rhs ? record[N] : T(0);
        const int index              = static_cast<int>(record[N + 1]);
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (below(tx, K, i)) {
                    if constexpr (fuse_rhs) rB[K] = K == slot ? rhs : rB[K];
                    origin[K] = K == slot ? index : origin[K];
                }
            }
        });
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<bool fuse_rhs, typename Sync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    physical_elimination_step(const int tx, const int i, T (&rA)[rows_per_thread][N],
                              [[maybe_unused]] T (&rB)[rows_per_thread],
                              int (&origin)[rows_per_thread], T* work, int& linfo,
                              Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        // The right-hand side moves with the rows only under S x = b.
        constexpr bool carries_rhs = fuse_rhs && !transposed;
        // Workspace: the candidates and their positions, the magnitude of
        // the row in place, then the exchange records of the rows leaving
        // and entering position i. Positions travel as values of T, as in
        // save_row.
        T* const candidate = work;
        T* const where     = candidate + threads_per_system;
        T* const in_place  = where + threads_per_system;
        T* const leaving   = in_place + 1;
        T* const entering  = leaving + N + 2;
 
        // Candidate of the thread. The owner of position i starts from the
        // row in place, so that a NaN there stops the search, as in the
        // logical scan, and publishes its magnitude.
        T best       = T(-1);
        int best_pos = N;
        T local      = T(0);
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                const int pos = tx + K * threads_per_system;
                if ((slot_is_full(K) || pos < N) && pos >= i) {
                    const T v = detail::abs(rA[K][i]);
                    if (pos == i || v > best) {
                        best     = v;
                        best_pos = pos;
                    }
                    if (pos == i) local = v;
                }
            }
        });
        candidate[tx] = best;
        where[tx]     = static_cast<T>(best_pos);
        if (tx == i % threads_per_system) in_place[0] = local;
        sync();
 
        // Every thread finds the same column maximum. The row in place keeps
        // the pivot when it reaches it, as in the logical scan, which keeps
        // the first maximum; otherwise the pivot is the first position
        // holding it. Only that rare case reads the positions. The loops
        // read the workspace, not the registers: no pragma.
        T rx_abs_max = candidate[i % threads_per_system];
        for (int t = 0; t < threads_per_system; t++)
            rx_abs_max = candidate[t] > rx_abs_max ? candidate[t] : rx_abs_max;
        bool keep = !(in_place[0] < rx_abs_max);
        if constexpr (relative_pivot_threshold < T(1)) keep = keeps_pivot(in_place[0], rx_abs_max);
        int max_id = i;
        if (!keep) {
            max_id = N;
            for (int t = 0; t < threads_per_system; t++) {
                const int pos = static_cast<int>(where[t]);
                if (candidate[t] == rx_abs_max && pos < max_id) max_id = pos;
            }
        }
        const bool zero_pivot = (rx_abs_max < singular_floor);
        linfo                 = (zero_pivot && linfo == 0) ? (i + 1) : linfo;
 
        // The owner of position i publishes its row whole, from a slot
        // known at compile time in an unrolled sweep. When the pivot is not
        // in place, the owner of position max_id publishes the pivot row as
        // well: it enters position i, and the row in place leaves it for
        // position max_id. The elimination reads the pivot row and its
        // right-hand-side entry from the record of the row that ends at
        // position i, under S x = b; under S^T x = b, the entry b'(i) of
        // the owner of position i, from the record of the row in place. A
        // zero pivot is not scaled away as in the logical scheme: the
        // factorization of a singular matrix is unspecified, and reg = 1
        // keeps it finite.
        at_position(tx, i, [&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            save_row<fuse_rhs>(rA[K], rB[K], origin[K], leaving);
        });
        if (max_id != i && tx == max_id % threads_per_system)
            gather_row<carries_rhs>(tx, i, max_id, rA, rB, origin, entering);
        sync();
 
        if (max_id != i) {
            at_position(tx, i, [&](auto K_tag) {
                constexpr int K = decltype(K_tag)::value;
                load_row<carries_rhs>(rA[K], rB[K], origin[K], entering);
            });
            if (tx == max_id % threads_per_system)
                scatter_row<carries_rhs>(tx, i, max_id, rA, rB, origin, leaving);
        }
 
        const T* const prow    = max_id == i ? leaving : entering;
        const T reg            = zero_pivot ? T(1) : T(1) / prow[i];
        [[maybe_unused]] T z_i = T(0);
        if constexpr (fuse_rhs && transposed) z_i = leaving[N] * reg;
 
        // scal and ger: every thread eliminates its rows below the pivot.
        // The next step writes only the candidates, which no thread reads
        // any more: no barrier closes the step.
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                const int pos = tx + K * threads_per_system;
                if ((slot_is_full(K) || pos < N) && pos > i) {
                    rA[K][i] *= reg;
                    if constexpr (Config.unroll_loops) {
                        TDLS_UNROLL_FORCE
                        for (int j = i + 1; j < N; j++)
                            rA[K][j] -= rA[K][i] * prow[j];
                    } else {
                        for (int j = i + 1; j < N; j++)
                            rA[K][j] -= rA[K][i] * prow[j];
                    }
                    if constexpr (carries_rhs) rB[K] -= rA[K][i] * prow[N];
                }
                if constexpr (fuse_rhs && transposed) {
                    if (pos == i) rB[K] = z_i;
                    if ((slot_is_full(K) || pos < N) && pos > i) rB[K] -= prow[pos] * z_i;
                }
            }
        });
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<bool fuse_rhs, typename Sync>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr int
    eliminate(const int tx, T (&rA)[rows_per_thread][N], T (&rB)[rows_per_thread],
              int (&rpiv)[rows_per_thread], T* work,
              Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        int linfo = 0;
        if constexpr (Config.row_interchange == RowInterchange::Logical) {
            T* const sB  = work;
            T* const sx  = work + N;
            T* const dsx = work + 2 * N;
            if constexpr (Config.unroll_loops) {
                TDLS_UNROLL_FORCE
                for (int i = 0; i < N; i++)
                    logical_elimination_step<fuse_rhs>(tx, i, rA, rB, rpiv, sB, sx, dsx, linfo,
                                                       sync);
            } else {
                for (int i = 0; i < N; i++)
                    logical_elimination_step<fuse_rhs>(tx, i, rA, rB, rpiv, sB, sx, dsx, linfo,
                                                       sync);
            }
        } else {
            if constexpr (Config.unroll_loops) {
                TDLS_UNROLL_FORCE
                for (int i = 0; i < N; i++)
                    physical_elimination_step<fuse_rhs>(tx, i, rA, rB, rpiv, work, linfo, sync);
            } else {
                for (int i = 0; i < N; i++)
                    physical_elimination_step<fuse_rhs>(tx, i, rA, rB, rpiv, work, linfo, sync);
            }
            sync(); // no barrier closes the steps: this one frees the workspace
        }
        return linfo;
    }
 
    /* =====================================================================
       SUBSTITUTION of S x = b, the column-major layout.
       Forward: the owner of the row at position i publishes its
       right-hand-side entry, every thread updates its rows below. This is
       the elimination of the right-hand side that solve_inplace folds into
       the factorization, replayed from the factored rows. Backward: the
       right-hand side is staged at the positions of its rows, then column
       i of U is published, every thread computes x_i, the owner of entry i
       keeps it and every thread updates its entries above. Under physical
       row interchanges, the right-hand side is first gathered in the
       pivoted order, and both passes run on the positions of the rows,
       which are static: the arithmetic is the logical one by construction,
       down to the contractions a compiler applies, and in an unrolled
       sweep the owner of each position is known at compile time.
       ===================================================================== */
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    forward_step(const int i, const T (&rA)[rows_per_thread][N], T (&rB)[rows_per_thread],
                 const int (&rowid)[rows_per_thread], T* sB,
                 Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (rowid[K] == i) sB[0] = rB[K];
            }
        });
        sync();
 
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if ((slot_is_full(K) || rowid[K] < N) && rowid[K] > i) rB[K] -= rA[K][i] * sB[0];
            }
        });
        sync();
    }
 
    template<typename Sync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    forward_pass(const T (&rA)[rows_per_thread][N], T (&rB)[rows_per_thread],
                 const int (&rowid)[rows_per_thread], T* work,
                 Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        T* const sB = work;
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = 0; i < N; i++)
                forward_step(i, rA, rB, rowid, sB, sync);
        } else {
            for (int i = 0; i < N; i++)
                forward_step(i, rA, rB, rowid, sB, sync);
        }
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    backward_step(const int i, const int tx, const T (&rA)[rows_per_thread][N],
                  T (&rB)[rows_per_thread], const int (&rowid)[rows_per_thread], T* sB, T* sx,
                  Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (slot_is_full(K) || rowid[K] < N) sx[rowid[K]] = rA[K][i];
            }
        });
        sync();
 
        const T reg = sB[i] / sx[i];
 
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                const int idx = tx + K * threads_per_system;
                if (slot_is_full(K) || idx < N) {
                    if (idx < i) sB[idx] = sB[idx] - reg * sx[idx];
                    rB[K] = (idx == i) ? reg : rB[K];
                }
            }
        });
        sync();
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    backward_pass(const int tx, const T (&rA)[rows_per_thread][N], T (&rB)[rows_per_thread],
                  const int (&rowid)[rows_per_thread], T* work,
                  Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        T* const sB = work;
        T* const sx = work + N;
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (slot_is_full(K) || rowid[K] < N) sB[rowid[K]] = rB[K];
            }
        });
        sync();
 
        if constexpr (Config.unroll_loops) {
            TDLS_UNROLL_FORCE
            for (int i = N - 1; i >= 0; i--)
                backward_step(i, tx, rA, rB, rowid, sB, sx, sync);
        } else {
            for (int i = N - 1; i >= 0; i--)
                backward_step(i, tx, rA, rB, rowid, sB, sx, sync);
        }
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    gather_rhs(const int tx, T (&rB)[rows_per_thread], const int (&origin)[rows_per_thread],
               T* staged, Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (slot_is_full(K) || tx + K * threads_per_system < N)
                    staged[tx + K * threads_per_system] = rB[K];
            }
        });
        sync();
 
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (slot_is_full(K) || tx + K * threads_per_system < N) rB[K] = staged[origin[K]];
            }
        });
    }
 
    /* =====================================================================
       SUBSTITUTION of S^T x = b, the row-major layout.
       With P S = L U, S^T = U^T L^T P. Both passes run on the positions of
       the rows, held at the entries of the same index: z_p and w_p live in
       the slot of entry p, which b'(p) leaves when it is consumed, so no
       second vector is needed. Forward, U^T z = b: the row at position i
       publishes its row of U, the owner of entry i its b'(i), every thread
       computes z_i = b'(i) / U(i, i) through the reciprocal of the pivot,
       the owner of entry i keeps z_i, and each entry c > i loses
       U(i, c) z_i. This is the forward pass that solve_inplace folds into
       the factorization, replayed from the factored rows. Backward,
       L^T w = z: the row at position k publishes its multipliers L(k, c),
       c < k, the owner of entry k its w_k, and each entry c < k loses
       L(k, c) w_k; L has a unit diagonal, so no division sits on the chain
       of the steps. One barrier per step: each step writes its own half of
       a double buffer, which no thread still reads (see segment). Then
       x = P^T w, one exchange through the workspace: under logical row
       interchanges, row r of S takes w at its position; under physical
       ones, w_p goes to the original index of the row at position p. These
       are the mirror of the staging of the right-hand side that the
       substitutions of S x = b make at the start.
       ===================================================================== */
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr T* segment(T* work,
                                                                       const int step) noexcept {
        return work + (step & 1) * (N + 1);
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_forward_step(const int tx, const int i, const T (&rA)[rows_per_thread][N],
                            T (&rB)[rows_per_thread], const int (&rowid)[rows_per_thread], T* work,
                            Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        T* const sx = segment(work, i);
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (rowid[K] == i) {
                    // Unrolled, the copy runs over all N entries under a
                    // guard: a trip count set by the step lets the AMD
                    // compiler unroll it before the step loop, and the
                    // rows then fall to scratch memory.
                    if constexpr (Config.unroll_loops) {
                        TDLS_UNROLL_FORCE
                        for (int j = 0; j < N; j++)
                            if (j >= i) sx[j] = rA[K][j];
                    } else {
                        for (int j = i; j < N; j++)
                            sx[j] = rA[K][j];
                    }
                }
            }
        });
        at_position(tx, i, [&](auto K_tag) { sx[N] = rB[decltype(K_tag)::value]; });
        sync();
 
        const T z_i = sx[N] * (T(1) / sx[i]);
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                const int c = tx + K * threads_per_system;
                if (c == i) rB[K] = z_i;
                if ((slot_is_full(K) || c < N) && c > i) rB[K] -= sx[c] * z_i;
            }
        });
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_backward_step(const int tx, const int k, const T (&rA)[rows_per_thread][N],
                             T (&rB)[rows_per_thread], const int (&rowid)[rows_per_thread], T* work,
                             Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        T* const sx = segment(work, k + 1);
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (rowid[K] == k) {
                    // Unrolled over all N entries under a guard, as in the
                    // forward step.
                    if constexpr (Config.unroll_loops) {
                        TDLS_UNROLL_FORCE
                        for (int c = 0; c < N; c++)
                            if (c < k) sx[c] = rA[K][c];
                    } else {
                        for (int c = 0; c < k; c++)
                            sx[c] = rA[K][c];
                    }
                }
            }
        });
        at_position(tx, k, [&](auto K_tag) { sx[N] = rB[decltype(K_tag)::value]; });
        sync();
 
        const T w_k = sx[N];
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                const int c = tx + K * threads_per_system;
                if (c < k) rB[K] -= sx[c] * w_k;
            }
        });
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<bool backward_only, typename Sync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    transposed_passes(const int tx, const T (&rA)[rows_per_thread][N], T (&rB)[rows_per_thread],
                      const int (&rowid)[rows_per_thread], T* work,
                      Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        if constexpr (Config.unroll_loops) {
            if constexpr (!backward_only) {
                TDLS_UNROLL_FORCE
                for (int i = 0; i < N; i++)
                    transposed_forward_step(tx, i, rA, rB, rowid, work, sync);
            }
            TDLS_UNROLL_FORCE
            for (int k = N - 1; k >= 1; k--)
                transposed_backward_step(tx, k, rA, rB, rowid, work, sync);
        } else {
            if constexpr (!backward_only) {
                for (int i = 0; i < N; i++)
                    transposed_forward_step(tx, i, rA, rB, rowid, work, sync);
            }
            for (int k = N - 1; k >= 1; k--)
                transposed_backward_step(tx, k, rA, rB, rowid, work, sync);
        }
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    unpermute_solution(const int tx, T (&rB)[rows_per_thread], const int (&rpiv)[rows_per_thread],
                       T* work, Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        constexpr bool logical = Config.row_interchange == RowInterchange::Logical;
        T* const staged        = segment(work, 1);
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (slot_is_full(K) || tx + K * threads_per_system < N)
                    staged[logical ? tx + K * threads_per_system : rpiv[K]] = rB[K];
            }
        });
        sync();
 
        detail::unroll_K<rows_per_thread>([&](auto K_tag) {
            constexpr int K = decltype(K_tag)::value;
            if constexpr (slot_has_rows(K)) {
                if (slot_is_full(K) || tx + K * threads_per_system < N)
                    rB[K] = staged[logical ? rpiv[K] : tx + K * threads_per_system];
            }
        });
    }
 
    /* =====================================================================
       SUBSTITUTION, dispatch on the layout.
       ===================================================================== */
 
    TDLS_EXEC_CHECK_DISABLE
    template<bool backward_only, typename Sync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitution(const int tx, const T (&rA)[rows_per_thread][N], T (&rB)[rows_per_thread],
                 const int (&rpiv)[rows_per_thread], T* work,
                 Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        constexpr bool logical = Config.row_interchange == RowInterchange::Logical;
        if constexpr (transposed) {
            if constexpr (logical) {
                transposed_passes<backward_only>(tx, rA, rB, rpiv, work, sync);
            } else {
                int position[rows_per_thread];
                init_piv(tx, position);
                transposed_passes<backward_only>(tx, rA, rB, position, work, sync);
            }
            unpermute_solution(tx, rB, rpiv, work, sync);
        } else if constexpr (logical) {
            if constexpr (!backward_only) forward_pass(rA, rB, rpiv, work, sync);
            backward_pass(tx, rA, rB, rpiv, work, sync);
        } else {
            int position[rows_per_thread];
            init_piv(tx, position);
            if constexpr (!backward_only) {
                // Staged apart from the broadcast slot of the pass, which a
                // faster thread may already write.
                gather_rhs(tx, rB, rpiv, work + N, sync);
                forward_pass(rA, rB, position, work, sync);
            }
            backward_pass(tx, rA, rB, position, work, sync);
        }
    }
 
    /* =====================================================================
       Entry points.
       ===================================================================== */
 
    TDLS_EXEC_CHECK_DISABLE
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr auto make_sync(T* work) noexcept {
        AutoSync deduce;
        return detail::resolve_sync<threads_per_system>(deduce, work, threads_per_system);
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<bool internal_piv, bool internal_matrix, int A_static_stride = dynamic_stride,
             int piv_static_stride = dynamic_stride, typename Sync = tdls::AutoSync>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    factorize(const int tx, T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
              const int piv_stride, T* work,
              Sync&& sync = Sync{}) noexcept(detail::nothrow_sync<Sync>) {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride>) {
            return factorize<internal_piv, internal_matrix>(
                tx, A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), work, std::forward<Sync>(sync));
        } else {
            auto&& group = detail::resolve_sync<threads_per_system>(sync, work, threads_per_system);
            T* const exchange = work + detail::barrier_elements;
            T rA[rows_per_thread][N];
            T rB[rows_per_thread];
            int rpiv[rows_per_thread];
            load_rows<internal_matrix>(tx, A, A_stride, rA);
            init_piv(tx, rpiv);
            const int linfo = eliminate<false>(tx, rA, rB, rpiv, exchange, group);
            store_rows<internal_matrix>(tx, A, A_stride, rA);
            store_piv<internal_piv>(tx, piv, piv_stride, rpiv);
            group(); // every output visible to the whole group on return
            return linfo == 0;
        }
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<bool internal_rhs, bool internal_piv, bool internal_matrix,
             int A_static_stride = dynamic_stride, int piv_static_stride = dynamic_stride,
             int rhs_static_stride = dynamic_stride, typename Sync = tdls::AutoSync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute(const int tx, 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, T* work,
               Sync&& sync = Sync{}) noexcept(detail::nothrow_sync<Sync>) {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return substitute<internal_rhs, internal_piv, internal_matrix>(
                tx, A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), b, x,
                detail::stride_or<rhs_static_stride>(rhs_stride), work, std::forward<Sync>(sync));
        } else {
            auto&& group = detail::resolve_sync<threads_per_system>(sync, work, threads_per_system);
            T* const exchange = work + detail::barrier_elements;
            T rA[rows_per_thread][N];
            T rB[rows_per_thread];
            int rpiv[rows_per_thread];
            load_rows<internal_matrix>(tx, A, A_stride, rA);
            load_piv<internal_piv>(tx, piv, piv_stride, rpiv);
            load_rhs<internal_rhs>(tx, b, rhs_stride, rB);
            substitution<false>(tx, rA, rB, rpiv, exchange, group);
            store_solution<internal_rhs>(tx, x, rhs_stride, rB);
            group(); // every output visible to the whole group on return
        }
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<bool internal_rhs, bool internal_piv, bool internal_matrix,
             int A_static_stride = dynamic_stride, int piv_static_stride = dynamic_stride,
             int rhs_static_stride = dynamic_stride, typename Sync = tdls::AutoSync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_canonical(const int tx, 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, T* work,
                         Sync&& sync = Sync{}) noexcept(detail::nothrow_sync<Sync>) {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return substitute_canonical<internal_rhs, internal_piv, internal_matrix>(
                tx, A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), col, x,
                detail::stride_or<rhs_static_stride>(rhs_stride), work, std::forward<Sync>(sync));
        } else {
            auto&& group = detail::resolve_sync<threads_per_system>(sync, work, threads_per_system);
            T* const exchange = work + detail::barrier_elements;
            T rA[rows_per_thread][N];
            T rB[rows_per_thread];
            int rpiv[rows_per_thread];
            load_rows<internal_matrix>(tx, A, A_stride, rA);
            load_piv<internal_piv>(tx, piv, piv_stride, rpiv);
            detail::unroll_K<rows_per_thread>([&](auto K_tag) {
                constexpr int K = decltype(K_tag)::value;
                rB[K]           = (tx + K * threads_per_system == col) ? T(1) : T(0);
            });
            substitution<false>(tx, rA, rB, rpiv, exchange, group);
            store_solution<internal_rhs>(tx, x, rhs_stride, rB);
            group(); // every output visible to the whole group on return
        }
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<bool internal_rhs, bool internal_piv, bool internal_matrix,
             int A_static_stride = dynamic_stride, int piv_static_stride = dynamic_stride,
             int rhs_static_stride = dynamic_stride, typename Sync = tdls::AutoSync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_inplace(const int tx, 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, T* work,
                       Sync&& sync = Sync{}) noexcept(detail::nothrow_sync<Sync>) {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return substitute_inplace<internal_rhs, internal_piv, internal_matrix>(
                tx, A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), x,
                detail::stride_or<rhs_static_stride>(rhs_stride), work, std::forward<Sync>(sync));
        } else {
            auto&& group = detail::resolve_sync<threads_per_system>(sync, work, threads_per_system);
            T* const exchange = work + detail::barrier_elements;
            T rA[rows_per_thread][N];
            T rB[rows_per_thread];
            int rpiv[rows_per_thread];
            load_rows<internal_matrix>(tx, A, A_stride, rA);
            load_piv<internal_piv>(tx, piv, piv_stride, rpiv);
            load_rhs<internal_rhs>(tx, x, rhs_stride, rB);
            substitution<false>(tx, rA, rB, rpiv, exchange, group);
            store_solution<internal_rhs>(tx, x, rhs_stride, rB);
            group(); // every output visible to the whole group on return
        }
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<bool internal_rhs, bool internal_piv, bool internal_matrix,
             int A_static_stride = dynamic_stride, int piv_static_stride = dynamic_stride,
             int rhs_static_stride = dynamic_stride, typename Sync = tdls::AutoSync>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve(const int tx, 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,
          T* work, Sync&& sync = Sync{}) noexcept(detail::nothrow_sync<Sync>) {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return solve<internal_rhs, internal_piv, internal_matrix>(
                tx, A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), b, x,
                detail::stride_or<rhs_static_stride>(rhs_stride), work, std::forward<Sync>(sync));
        } else {
            auto&& group = detail::resolve_sync<threads_per_system>(sync, work, threads_per_system);
            T* const exchange = work + detail::barrier_elements;
            T rA[rows_per_thread][N];
            T rB[rows_per_thread];
            int rpiv[rows_per_thread];
            load_rows<internal_matrix>(tx, A, A_stride, rA);
            init_piv(tx, rpiv);
            const int linfo = eliminate<false>(tx, rA, rB, rpiv, exchange, group);
            store_rows<internal_matrix>(tx, A, A_stride, rA);
            store_piv<internal_piv>(tx, piv, piv_stride, rpiv);
            load_rhs<internal_rhs>(tx, b, rhs_stride, rB);
            substitution<false>(tx, rA, rB, rpiv, exchange, group);
            store_solution<internal_rhs>(tx, x, rhs_stride, rB);
            group(); // every output visible to the whole group on return
            return linfo == 0;
        }
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<bool internal_rhs, bool internal_piv, bool internal_matrix,
             int A_static_stride = dynamic_stride, int piv_static_stride = dynamic_stride,
             int rhs_static_stride = dynamic_stride, typename Sync = tdls::AutoSync>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve_inplace(const int tx, T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv,
                  const int piv_stride, T* TDLS_RESTRICT y, const int rhs_stride, T* work,
                  Sync&& sync = Sync{}) noexcept(detail::nothrow_sync<Sync>) {
        if constexpr (detail::any_static_stride<A_static_stride, piv_static_stride,
                                                rhs_static_stride>) {
            return solve_inplace<internal_rhs, internal_piv, internal_matrix>(
                tx, A, detail::stride_or<A_static_stride>(A_stride), piv,
                detail::stride_or<piv_static_stride>(piv_stride), y,
                detail::stride_or<rhs_static_stride>(rhs_stride), work, std::forward<Sync>(sync));
        } else {
            auto&& group = detail::resolve_sync<threads_per_system>(sync, work, threads_per_system);
            T* const exchange = work + detail::barrier_elements;
            T rA[rows_per_thread][N];
            T rB[rows_per_thread];
            int rpiv[rows_per_thread];
            load_rows<internal_matrix>(tx, A, A_stride, rA);
            init_piv(tx, rpiv);
            load_rhs<internal_rhs>(tx, y, rhs_stride, rB);
            const int linfo = eliminate<true>(tx, rA, rB, rpiv, exchange, group);
            store_rows<internal_matrix>(tx, A, A_stride, rA);
            store_piv<internal_piv>(tx, piv, piv_stride, rpiv);
            substitution<true>(tx, rA, rB, rpiv, exchange, group);
            detail::unroll_K<rows_per_thread>([&](auto K_tag) {
                constexpr int K = decltype(K_tag)::value;
                if constexpr (slot_has_rows(K)) {
                    if (slot_is_full(K) || tx + K * threads_per_system < N)
                        TDLS_COOP_LUPP_Y(K) = rB[K];
                }
            });
            group(); // every output visible to the whole group on return
            return linfo == 0;
        }
    }
};
 
 
 
#undef TDLS_COOP_LUPP_A
#undef TDLS_COOP_LUPP_PIV
#undef TDLS_COOP_LUPP_B
#undef TDLS_COOP_LUPP_X
#undef TDLS_COOP_LUPP_Y
 
 
 
} // namespace tdls
 
#if defined(TDLS_DETAIL_COMPILER_NVCXX)
#pragma diag_pop
#endif
 
#if defined(TDLS_DETAIL_COMPILER_GCC)
#pragma GCC diagnostic pop
#endif
 
#if defined(TDLS_DETAIL_COMPILER_CLANG)
#pragma clang diagnostic pop
#endif
 
 
 
#endif // TDLS_SOLVERS_COOPERATIVE_LUPP_SOLVER_STATIC_HPP