solver_dynamic.hpp

Go to the documentation of this file.

#ifndef TDLS_SOLVERS_COOPERATIVE_LUPP_SOLVER_DYNAMIC_HPP
#define TDLS_SOLVERS_COOPERATIVE_LUPP_SOLVER_DYNAMIC_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 <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>
 
 
 
// gcc's flow analysis cannot prove that the right-hand-side entries of a
// thread are loaded before being read: the load stops at the first
// phantom row and every later use stops at the same row, a property
// validated bitwise against the compile-time solver. The resulting
// -Wmaybe-uninitialized reports are spurious. Under UBSan instrumentation
// at -O2 and above, gcc additionally loses the value ranges of the 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,
// exactly as in solver_static.hpp.
#if defined(TDLS_DETAIL_COMPILER_GCC)
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
#pragma GCC diagnostic ignored "-Warray-bounds"
#endif
 
namespace tdls {
 
 
 
/* Addressing macros. They expand inside member functions where n, the
   parameter names (A, A_stride, piv, piv_stride, b, x, y, rhs_stride)
   and the Config value are in scope. #undef'd at the end of this
   header. */
 
#define TDLS_COOP_LUPP_DYN_A(r, c) A[TDLS_STRIDED((c) * n + (r), A_stride)]
#define TDLS_COOP_LUPP_DYN_PIV(r) piv[TDLS_STRIDED(r, piv_stride)]
#define TDLS_COOP_LUPP_DYN_B(r) b[TDLS_STRIDED(r, rhs_stride)]
#define TDLS_COOP_LUPP_DYN_X(r) x[TDLS_STRIDED(r, rhs_stride)]
#define TDLS_COOP_LUPP_DYN_Y(r) y[TDLS_STRIDED(r, rhs_stride)]
 
 
 
template<typename T, CooperativeLUppConfig<T> Config = CooperativeLUppConfig<T>{}>
struct CooperativeLUppSolverDynamic {
 
    static_assert(
        Config.row_interchange != RowInterchange::Auto,
        "CooperativeLUppSolverDynamic: row_interchange must be Logical or Physical (Auto is the "
        "default of the blocks of ArrowheadLUpp only)");
    static_assert(Config.rows_per_thread >= 1,
                  "CooperativeLUppSolverDynamic: rows_per_thread must be >= 1");
    static_assert(Config.relative_pivot_threshold.is_finite() && Config.singular_floor.is_finite(),
                  "CooperativeLUppSolverDynamic: 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;
    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),
                  "CooperativeLUppSolverDynamic: relative_pivot_threshold must be in (0, 1]");
    static_assert(singular_floor > T(0),
                  "CooperativeLUppSolverDynamic: singular_floor must be positive");
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr int
    threads_per_system(const int n) noexcept {
        return (n + rows_per_thread - 1) / rows_per_thread;
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr int
    workspace_size(const int n) noexcept {
        return detail::barrier_elements + (Config.row_interchange == RowInterchange::Logical
                                               ? 3 * n
                                               : 2 * n + 2 * threads_per_system(n) + 5);
    }
 
    /* =====================================================================
       Pivot entries and right-hand-side entries.
       Slot K of thread tx holds row tx + K * threads, a phantom row when
       it reaches n. The rows themselves stay in the matrix; only their
       pivot entries and their right-hand-side entries are kept by the
       thread. Each thread reads and writes only the entries of its own
       rows, so the operands never carry data between threads: the
       workspace does.
       ===================================================================== */
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    init_piv(const int tx, const int threads, int (&rpiv)[rows_per_thread]) noexcept {
        for (int K = 0; K < rows_per_thread; ++K)
            rpiv[K] = tx + K * threads;
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_piv(const int n, const int tx, const int threads, const int* TDLS_RESTRICT piv,
             const int piv_stride, int (&rpiv)[rows_per_thread]) noexcept {
        for (int K = 0; K < rows_per_thread; ++K) {
            const int r = tx + K * threads;
            rpiv[K]     = r < n ? TDLS_COOP_LUPP_DYN_PIV(r) : r;
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    store_piv(const int n, const int tx, const int threads, int* TDLS_RESTRICT piv,
              const int piv_stride, const int (&rpiv)[rows_per_thread]) noexcept {
        for (int K = 0; K < rows_per_thread; ++K) {
            const int r = tx + K * threads;
            if (r >= n) break;
            TDLS_COOP_LUPP_DYN_PIV(r) = rpiv[K];
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_rhs(const int n, const int tx, const int threads, const T* TDLS_RESTRICT b,
             const int rhs_stride, T (&rB)[rows_per_thread]) noexcept {
        for (int K = 0; K < rows_per_thread; ++K) {
            const int r = tx + K * threads;
            if (r >= n) break;
            rB[K] = TDLS_COOP_LUPP_DYN_B(r);
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    store_solution(const int n, const int tx, const int threads, T* TDLS_RESTRICT x,
                   const int rhs_stride, const T (&rB)[rows_per_thread]) noexcept {
        for (int K = 0; K < rows_per_thread; ++K) {
            const int r = tx + K * threads;
            if (r >= n) break;
            TDLS_COOP_LUPP_DYN_X(r) = rB[K];
        }
    }
 
    /* =====================================================================
       FACTORIZATION, right-looking, one column per step, as in the
       compile-time solver, in place in the matrix. Logical row
       interchanges: 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: 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.
       ===================================================================== */
 
    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(const int n, const int i, const int tx, const int threads,
                             T* TDLS_RESTRICT A, const int A_stride,
                             [[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.
        for (int K = 0; K < rows_per_thread; ++K) {
            const int r = tx + K * threads;
            if (r >= n) break;
            dsx[rowid[K]] = detail::abs(TDLS_COOP_LUPP_DYN_A(r, i));
        }
        sync();
 
        // Every thread scans the same magnitudes and finds the same pivot.
        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.
        for (int K = 0; K < rows_per_thread; ++K) {
            const int r = tx + K * threads;
            if (r >= n) break;
            if (rowid[K] == max_id) {
                rowid[K] = i;
                for (int j = i; j < n; j++)
                    sx[j] = update * TDLS_COOP_LUPP_DYN_A(r, j);
                if constexpr (fuse_rhs && !transposed) sB[0] = rB[K];
            } else if (rowid[K] == i) {
                rowid[K] = max_id;
            }
            // S^T x = b: the owner of entry i publishes b'(i), final at
            // this step.
            if constexpr (fuse_rhs && transposed) {
                if (r == i) sB[0] = rB[K];
            }
        }
        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.
        for (int K = 0; K < rows_per_thread; ++K) {
            const int r = tx + K * threads;
            if (r >= n) break;
            if (rowid[K] > i) {
                TDLS_COOP_LUPP_DYN_A(r, i) *= reg;
                for (int j = i + 1; j < n; j++)
                    TDLS_COOP_LUPP_DYN_A(r, j) -= TDLS_COOP_LUPP_DYN_A(r, i) * sx[j];
                if constexpr (fuse_rhs && !transposed) rB[K] -= TDLS_COOP_LUPP_DYN_A(r, i) * sB[0];
            }
            // the owner of entry i keeps z_i in its place; each entry c > i
            // loses U(i, c) z_i
            if constexpr (fuse_rhs && transposed) {
                if (r == i) rB[K] = z_i;
                if (r > i) rB[K] -= sx[r] * z_i;
            }
        }
        sync();
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    save_row(const int n, const T* TDLS_RESTRICT A, const int A_stride, const int r,
             T* record) noexcept {
        for (int c = 0; c < n; c++)
            record[c] = TDLS_COOP_LUPP_DYN_A(r, c);
    }
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    load_row(const int n, T* TDLS_RESTRICT A, const int A_stride, const int r,
             const T* record) noexcept {
        for (int c = 0; c < n; c++)
            TDLS_COOP_LUPP_DYN_A(r, c) = record[c];
    }
 
    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 n, const int i, const int tx, const int threads,
                              T* TDLS_RESTRICT A, const int A_stride,
                              [[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; under
        // S^T x = b, the owner of position i is the owner of entry i: its
        // record carries b'(i), and it keeps z_i in its place.
        constexpr bool carries_rhs = fuse_rhs && !transposed;
        // Workspace layout of the compile-time solver.
        T* const candidate = work;
        T* const where     = candidate + threads;
        T* const in_place  = where + threads;
        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 and publishes its magnitude, as in the compile-time
        // solver.
        T best       = T(-1);
        int best_pos = n;
        T local      = T(0);
        for (int K = 0; K < rows_per_thread; ++K) {
            const int pos = tx + K * threads;
            if (pos >= n) break;
            if (pos < i) continue;
            const T v = detail::abs(TDLS_COOP_LUPP_DYN_A(pos, 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) in_place[0] = local;
        sync();
 
        // Every thread finds the same column maximum and the same pivot, as
        // in the compile-time solver.
        T rx_abs_max = candidate[i % threads];
        for (int t = 0; t < threads; 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; 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;
 
        // As in the compile-time solver: the owner of position i publishes
        // its row whole and, when the pivot is not in place, the owner of
        // position max_id publishes the pivot row. Each record holds the n
        // entries of its row, then its right-hand-side entry (fuse_rhs
        // only) and its original index. The rows are addressed by position,
        // and the entries of the slots move through loops kept small, which
        // the compiler unrolls: a slot indexed at run time would demote rB
        // and origin to the stack.
        if (tx == i % threads) save_row(n, A, A_stride, i, leaving);
        if (max_id != i && tx == max_id % threads) save_row(n, A, A_stride, max_id, entering);
        for (int K = 0; K < rows_per_thread; ++K) {
            const int pos = tx + K * threads;
            if (pos != i && (pos != max_id || max_id == i)) continue;
            T* const record = pos == i ? leaving : entering;
            if constexpr (fuse_rhs) record[n] = rB[K];
            record[n + 1] = static_cast<T>(origin[K]);
        }
        sync();
 
        if (max_id != i) {
            if (tx == i % threads) load_row(n, A, A_stride, i, entering);
            if (tx == max_id % threads) load_row(n, A, A_stride, max_id, leaving);
            for (int K = 0; K < rows_per_thread; ++K) {
                const int pos = tx + K * threads;
                if (pos != max_id && pos != i) continue;
                const T* const record = pos == i ? entering : leaving;
                if constexpr (carries_rhs) rB[K] = record[n];
                origin[K] = static_cast<int>(record[n + 1]);
            }
        }
 
        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.
        for (int K = 0; K < rows_per_thread; ++K) {
            const int pos = tx + K * threads;
            if (pos >= n) break;
            if (pos > i) {
                TDLS_COOP_LUPP_DYN_A(pos, i) *= reg;
                for (int j = i + 1; j < n; j++)
                    TDLS_COOP_LUPP_DYN_A(pos, j) -= TDLS_COOP_LUPP_DYN_A(pos, i) * prow[j];
                if constexpr (carries_rhs) rB[K] -= TDLS_COOP_LUPP_DYN_A(pos, i) * prow[n];
            }
            if constexpr (fuse_rhs && transposed) {
                if (pos == i) rB[K] = z_i;
                if (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 n, const int tx, const int threads, T* TDLS_RESTRICT A, const int A_stride,
              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;
            for (int i = 0; i < n; i++)
                logical_elimination_step<fuse_rhs>(n, i, tx, threads, A, A_stride, rB, rpiv, sB, sx,
                                                   dsx, linfo, sync);
        } else {
            for (int i = 0; i < n; i++)
                physical_elimination_step<fuse_rhs>(n, i, tx, threads, A, A_stride, 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, as in the
       compile-time solver: the forward pass
       replays the elimination of the right-hand side from the factored
       rows, the backward pass publishes column i of U, 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.
       ===================================================================== */
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    forward_pass(const int n, const int tx, const int threads, const T* TDLS_RESTRICT A,
                 const int A_stride, T (&rB)[rows_per_thread], const int (&rowid)[rows_per_thread],
                 T* work, Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        T* const sB = work;
        for (int i = 0; i < n; i++) {
            for (int K = 0; K < rows_per_thread; ++K) {
                if (tx + K * threads >= n) break;
                if (rowid[K] == i) sB[0] = rB[K];
            }
            sync();
 
            for (int K = 0; K < rows_per_thread; ++K) {
                const int r = tx + K * threads;
                if (r >= n) break;
                if (rowid[K] > i) rB[K] -= TDLS_COOP_LUPP_DYN_A(r, i) * sB[0];
            }
            sync();
        }
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    backward_pass(const int n, const int tx, const int threads, const T* TDLS_RESTRICT A,
                  const int A_stride, 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;
        for (int K = 0; K < rows_per_thread; ++K) {
            if (tx + K * threads >= n) break;
            sB[rowid[K]] = rB[K];
        }
        sync();
 
        for (int i = n - 1; i >= 0; i--) {
            for (int K = 0; K < rows_per_thread; ++K) {
                const int r = tx + K * threads;
                if (r >= n) break;
                sx[rowid[K]] = TDLS_COOP_LUPP_DYN_A(r, i);
            }
            sync();
 
            const T reg = sB[i] / sx[i];
 
            for (int K = 0; K < rows_per_thread; ++K) {
                const int idx = tx + K * threads;
                if (idx >= n) break;
                if (idx < i) sB[idx] = sB[idx] - reg * sx[idx];
                rB[K] = (idx == i) ? reg : rB[K];
            }
            sync();
        }
    }
 
    /* =====================================================================
       SUBSTITUTION of S^T x = b, the row-major layout: the passes of the
       compile-time solver, U^T z = b then L^T w = z on the positions held
       at the entries of the same index, one barrier per step, then
       x = P^T w, with the rows of S in memory. The positions of the rows
       are published once; then each thread loads the entries of S it needs
       itself, before the barrier of the step, so that the latency of the
       memory sits off the chain of the steps. The arithmetic is that of the
       compile-time solver. Workspace: the row of S at each position
       (logical row interchanges), then a double-buffered broadcast value,
       n + 2 elements; the first n elements stage x = P^T w at the end,
       once no step reads them any more.
       ===================================================================== */
 
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr int
    row_at([[maybe_unused]] const T* rows, const int k) noexcept {
        if constexpr (Config.row_interchange == RowInterchange::Logical)
            return static_cast<int>(rows[k]);
        else
            return 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 n, const int tx, const int threads, const T* TDLS_RESTRICT A,
                      const int A_stride, T (&rB)[rows_per_thread], const T* rows, T* buf,
                      Sync& sync) noexcept(detail::nothrow_sync<Sync>) {
        if constexpr (!backward_only) {
            for (int i = 0; i < n; i++) {
                const int ri = row_at(rows, i);
                const T u_ii = TDLS_COOP_LUPP_DYN_A(ri, i);
                T u[rows_per_thread];
                for (int K = 0; K < rows_per_thread; ++K) {
                    const int c = tx + K * threads;
                    if (c >= n) break;
                    u[K] = c > i ? TDLS_COOP_LUPP_DYN_A(ri, c) : T(0);
                    if (c == i) buf[i & 1] = rB[K];
                }
                sync();
 
                const T z_i = buf[i & 1] * (T(1) / u_ii);
                for (int K = 0; K < rows_per_thread; ++K) {
                    const int c = tx + K * threads;
                    if (c >= n) break;
                    if (c == i) rB[K] = z_i;
                    if (c > i) rB[K] -= u[K] * z_i;
                }
            }
        }
        for (int k = n - 1; k >= 1; k--) {
            const int rk = row_at(rows, k);
            T l[rows_per_thread];
            for (int K = 0; K < rows_per_thread; ++K) {
                const int c = tx + K * threads;
                if (c >= n) break;
                l[K] = c < k ? TDLS_COOP_LUPP_DYN_A(rk, c) : T(0);
                if (c == k) buf[(k + 1) & 1] = rB[K];
            }
            sync();
 
            const T w_k = buf[(k + 1) & 1];
            for (int K = 0; K < rows_per_thread; ++K) {
                const int c = tx + K * threads;
                if (c >= n) break;
                if (c < k) rB[K] -= l[K] * w_k;
            }
        }
    }
 
    /* =====================================================================
       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 n, const int tx, const int threads, const T* TDLS_RESTRICT A,
                 const int A_stride, 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) {
            // the row of S at each position
            if constexpr (logical) {
                for (int K = 0; K < rows_per_thread; ++K) {
                    const int r = tx + K * threads;
                    if (r >= n) break;
                    work[rpiv[K]] = static_cast<T>(r);
                }
                sync();
            }
            transposed_passes<backward_only>(n, tx, threads, A, A_stride, rB, work, work + n, sync);
            // x = P^T w: 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.
            T* const staged = work;
            for (int K = 0; K < rows_per_thread; ++K) {
                const int c = tx + K * threads;
                if (c >= n) break;
                staged[logical ? c : rpiv[K]] = rB[K];
            }
            sync();
            for (int K = 0; K < rows_per_thread; ++K) {
                const int c = tx + K * threads;
                if (c >= n) break;
                rB[K] = staged[logical ? rpiv[K] : c];
            }
        } else if constexpr (logical) {
            if constexpr (!backward_only)
                forward_pass(n, tx, threads, A, A_stride, rB, rpiv, work, sync);
            backward_pass(n, tx, threads, A, A_stride, rB, rpiv, work, sync);
        } else {
            int position[rows_per_thread];
            init_piv(tx, threads, position);
            if constexpr (!backward_only) {
                // Gather in the pivoted order: each slot receives the entry
                // of the original row it holds. The entries are staged apart
                // from the broadcast slot of the pass, which a faster thread
                // may already write.
                T* const staged = work + n;
                for (int K = 0; K < rows_per_thread; ++K) {
                    const int r = tx + K * threads;
                    if (r >= n) break;
                    staged[r] = rB[K];
                }
                sync();
                for (int K = 0; K < rows_per_thread; ++K) {
                    if (tx + K * threads >= n) break;
                    rB[K] = staged[rpiv[K]];
                }
                forward_pass(n, tx, threads, A, A_stride, rB, position, work, sync);
            }
            backward_pass(n, tx, threads, A, A_stride, rB, position, work, sync);
        }
    }
 
    /* =====================================================================
       Entry points.
       ===================================================================== */
 
    TDLS_EXEC_CHECK_DISABLE
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr auto make_sync(const int n,
                                                                           T* work) noexcept {
        AutoSync deduce;
        return detail::resolve_sync<0>(deduce, work, threads_per_system(n));
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync = tdls::AutoSync>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    factorize(const int n, 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>) {
        const int threads = threads_per_system(n);
        auto&& group      = detail::resolve_sync<0>(sync, work, threads);
        T* const exchange = work + detail::barrier_elements;
        T rB[rows_per_thread];
        int rpiv[rows_per_thread];
        init_piv(tx, threads, rpiv);
        const int linfo = eliminate<false>(n, tx, threads, A, A_stride, rB, rpiv, exchange, group);
        store_piv(n, tx, threads, piv, piv_stride, rpiv);
        group(); // every output visible to the whole group on return
        return linfo == 0;
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync = tdls::AutoSync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute(const int n, 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>) {
        const int threads = threads_per_system(n);
        auto&& group      = detail::resolve_sync<0>(sync, work, threads);
        T* const exchange = work + detail::barrier_elements;
        T rB[rows_per_thread];
        int rpiv[rows_per_thread];
        load_piv(n, tx, threads, piv, piv_stride, rpiv);
        load_rhs(n, tx, threads, b, rhs_stride, rB);
        substitution<false>(n, tx, threads, A, A_stride, rB, rpiv, exchange, group);
        store_solution(n, tx, threads, x, rhs_stride, rB);
        group(); // every output visible to the whole group on return
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync = tdls::AutoSync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_canonical(const int n, 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>) {
        const int threads = threads_per_system(n);
        auto&& group      = detail::resolve_sync<0>(sync, work, threads);
        T* const exchange = work + detail::barrier_elements;
        T rB[rows_per_thread];
        int rpiv[rows_per_thread];
        load_piv(n, tx, threads, piv, piv_stride, rpiv);
        for (int K = 0; K < rows_per_thread; ++K)
            rB[K] = (tx + K * threads == col) ? T(1) : T(0);
        substitution<false>(n, tx, threads, A, A_stride, rB, rpiv, exchange, group);
        store_solution(n, tx, threads, x, rhs_stride, rB);
        group(); // every output visible to the whole group on return
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync = tdls::AutoSync>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
    substitute_inplace(const int n, 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>) {
        const int threads = threads_per_system(n);
        auto&& group      = detail::resolve_sync<0>(sync, work, threads);
        T* const exchange = work + detail::barrier_elements;
        T rB[rows_per_thread];
        int rpiv[rows_per_thread];
        load_piv(n, tx, threads, piv, piv_stride, rpiv);
        load_rhs(n, tx, threads, x, rhs_stride, rB);
        substitution<false>(n, tx, threads, A, A_stride, rB, rpiv, exchange, group);
        store_solution(n, tx, threads, x, rhs_stride, rB);
        group(); // every output visible to the whole group on return
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync = tdls::AutoSync>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve(const int n, 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>) {
        const int threads = threads_per_system(n);
        auto&& group      = detail::resolve_sync<0>(sync, work, threads);
        T* const exchange = work + detail::barrier_elements;
        T rB[rows_per_thread];
        int rpiv[rows_per_thread];
        init_piv(tx, threads, rpiv);
        const int linfo = eliminate<false>(n, tx, threads, A, A_stride, rB, rpiv, exchange, group);
        store_piv(n, tx, threads, piv, piv_stride, rpiv);
        load_rhs(n, tx, threads, b, rhs_stride, rB);
        substitution<false>(n, tx, threads, A, A_stride, rB, rpiv, exchange, group);
        store_solution(n, tx, threads, x, rhs_stride, rB);
        group(); // every output visible to the whole group on return
        return linfo == 0;
    }
 
    TDLS_EXEC_CHECK_DISABLE
    template<typename Sync = tdls::AutoSync>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
    solve_inplace(const int n, 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>) {
        const int threads = threads_per_system(n);
        auto&& group      = detail::resolve_sync<0>(sync, work, threads);
        T* const exchange = work + detail::barrier_elements;
        T rB[rows_per_thread];
        int rpiv[rows_per_thread];
        init_piv(tx, threads, rpiv);
        load_rhs(n, tx, threads, y, rhs_stride, rB);
        const int linfo = eliminate<true>(n, tx, threads, A, A_stride, rB, rpiv, exchange, group);
        store_piv(n, tx, threads, piv, piv_stride, rpiv);
        substitution<true>(n, tx, threads, A, A_stride, rB, rpiv, exchange, group);
        for (int K = 0; K < rows_per_thread; ++K) {
            const int r = tx + K * threads;
            if (r >= n) break;
            TDLS_COOP_LUPP_DYN_Y(r) = rB[K];
        }
        group(); // every output visible to the whole group on return
        return linfo == 0;
    }
};
 
 
 
#undef TDLS_COOP_LUPP_DYN_A
#undef TDLS_COOP_LUPP_DYN_PIV
#undef TDLS_COOP_LUPP_DYN_B
#undef TDLS_COOP_LUPP_DYN_X
#undef TDLS_COOP_LUPP_DYN_Y
 
 
 
} // namespace tdls
 
#if defined(TDLS_DETAIL_COMPILER_GCC)
#pragma GCC diagnostic pop
#endif
 
 
 
#endif // TDLS_SOLVERS_COOPERATIVE_LUPP_SOLVER_DYNAMIC_HPP