solver_dynamic.hpp

Go to the documentation of this file.

#ifndef TDLS_SOLVERS_ARROWHEAD_LUPP_SOLVER_DYNAMIC_HPP
#define TDLS_SOLVERS_ARROWHEAD_LUPP_SOLVER_DYNAMIC_HPP
 
 
 
 
 
 
#include <limits>
 
#include <tdls/core/macros.hpp>
#include <tdls/core/math.hpp>
#include <tdls/solvers/arrowhead_lupp/config.hpp>
#include <tdls/solvers/arrowhead_lupp/storage_map.hpp>
#include <tdls/solvers/tiled_lupp/solver_dynamic.hpp>
 
 
 
namespace tdls {
 
 
 
/* Addressing macros. They expand inside member functions where the
   parameter names (A, A_stride, y, rhs_stride, work, work_stride, piv,
   piv_stride) are in scope. #undef'd at the end of this header. */
 
#define TDLS_ARROW_LUPP_DYN_A(idx) A[TDLS_STRIDED(idx, A_stride)]
#define TDLS_ARROW_LUPP_DYN_Y(i) y[TDLS_STRIDED(i, rhs_stride)]
#define TDLS_ARROW_LUPP_DYN_YW(m, i)                                                               \
    y[TDLS_STRIDED(i, rhs_stride) +                                                                \
      (ycols == 1 ? std::size_t(0) : TDLS_STRIDED_WIDE(m, xcol_stride))]
#define TDLS_ARROW_LUPP_DYN_W(e) work[TDLS_STRIDED(e, work_stride)]
#define TDLS_ARROW_LUPP_DYN_F(r, c) TDLS_ARROW_LUPP_DYN_W((r) * front_cols + (c))
#define TDLS_ARROW_LUPP_DYN_LU(i, j) TDLS_ARROW_LUPP_DYN_W(o.Dl + (i) * S + (j))
#define TDLS_ARROW_LUPP_DYN_RC(c, i) TDLS_ARROW_LUPP_DYN_W(o.W + (c) * S + (i))
 
 
 
template<typename T, ArrowheadLUppConfig<T> Config>
struct ArrowheadLUppSolverDynamic {
    using storage_map = ArrowheadStorageMap<Config.storage, Config.coupling, Config.layout>;
 
  private:
    static constexpr bool col  = Config.layout == MatrixLayout::ColMajor;
    static constexpr bool diag = Config.coupling == ArrowheadCoupling::Diagonal;
    static constexpr bool multiplier_update =
        !col && !diag && Config.condensation == ArrowheadCondensation::MinFlops;
    static constexpr bool in_place = Config.schur_in_place;
 
    static_assert(Config.condensation_pivot_threshold.is_finite() &&
                      T(Config.condensation_pivot_threshold) >= T(0) &&
                      T(Config.condensation_pivot_threshold) <= T(1),
                  "ArrowheadLUppSolverDynamic: condensation_pivot_threshold must be in [0, 1]");
 
    static constexpr T condensation_pivot_threshold = Config.condensation_pivot_threshold;
    static_assert(Config.rhs_pass_width >= -1,
                  "ArrowheadLUppSolverDynamic: rhs_pass_width must be >= -1");
 
    using BlockSolver  = TiledLUppSolverDynamic<T, detail::arrowhead_block_config(Config)>;
    using BorderSolver = TiledLUppSolverDynamic<T, detail::arrowhead_border_config(Config)>;
    using PhysicalBlockSolver =
        TiledLUppSolverDynamic<T, detail::arrowhead_physical_block_config(Config)>;
 
    struct Sizes {
        int B, S, NB, R, R0, front_rows, front_cols, extra_cols, pass_width;
        bool reg;
    };
    template<typename I>
    struct OffsetsOf {
        I extra, Dl, U, W, X, Mb, m, a, acoef, Sacc, racc, Sbb, rb, w, s, refused, block_row, total;
    };
    using Offsets = OffsetsOf<int>;
    template<typename I>
    struct StreamedOffsetsOf {
        I Abk, Dsave, Wsave, total;
    };
    using StreamedOffsets = StreamedOffsetsOf<int>;
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr Sizes
    sizes(const ArrowheadShape& shape, const int nc, const int ny) noexcept {
        Sizes z{};
        z.B                                  = shape.border;
        z.S                                  = shape.block;
        z.NB                                 = shape.blocks;
        z.R                                  = diag ? z.S : z.B;
        z.R0                                 = diag ? shape.coupling_row : 0;
        z.front_rows                         = z.S + z.B;
        z.front_cols                         = z.S + z.B + z.R + nc;
        z.extra_cols                         = z.R > z.S ? z.R - z.S : 0;
        constexpr RowInterchange interchange = Config.block_row_interchange;
        constexpr int width                  = Config.rhs_pass_width;
        z.reg          = detail::arrowhead_register_blocks_here<T>(interchange, z.S);
        const int pass = detail::arrowhead_rhs_pass_width_here(width, z.reg);
        z.pass_width   = pass == 0 || pass > z.B + ny ? z.B + ny : pass;
        return z;
    }
 
    template<typename I = int>
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr OffsetsOf<I>
    offsets(const ArrowheadShape& shape, const int nc) noexcept {
        const Sizes z = sizes(shape, nc, nc);
        OffsetsOf<I> o{};
        I e     = I(z.front_rows) * z.front_cols;
        o.extra = e;
        e += I(z.NB) * z.S * z.extra_cols;
        o.Dl = e;
        e += I(z.S) * z.S;
        o.U = e;
        e += I(z.S) * z.S;
        o.W = e;
        e += (I(z.B) + nc) * z.S;
        o.X = e;
        e += (I(z.B) + nc) * z.S;
        o.Mb = e;
        e += I(z.R) * z.S;
        o.m = e;
        e += z.S;
        o.a = e;
        e += z.S;
        o.acoef = e;
        e += z.S;
        o.Sacc = e;
        e += I(z.R) * z.B;
        o.racc = e;
        e += I(z.R) * nc;
        o.Sbb = e;
        e += I(z.B) * z.B;
        o.rb = e;
        e += I(z.B) * nc;
        o.w = e;
        e += I(z.R) * nc;
        o.s = e;
        e += I(z.S) * nc;
        o.refused = e;
        e += z.NB;
        o.block_row = e;
        e += z.front_rows;
        o.total = e;
        return o;
    }
 
    template<typename I = int>
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr StreamedOffsetsOf<I>
    streamed_offsets(const ArrowheadShape& shape, const int nc) noexcept {
        const Sizes z = sizes(shape, nc, nc);
        StreamedOffsetsOf<I> so{};
        I e    = offsets<I>(shape, nc).total;
        so.Abk = e;
        e += diag ? I(0) : I(z.B) * z.S;
        so.Dsave = e;
        e += I(z.S) * z.S;
        so.Wsave = e;
        e += (I(z.B) + 1) * z.S;
        so.total = e;
        return so;
    }
 
  public:
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool valid(const ArrowheadShape& shape,
                                                               const int nrhs = 1) noexcept {
        return storage_map::valid(shape) && nrhs >= 1 && nrhs <= (1 << 30) &&
               streamed_offsets<long long>(shape, nrhs).total < (1LL << 31);
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr int
    workspace_size(const ArrowheadShape& shape) noexcept {
        return valid(shape) ? offsets(shape, 1).total : 0;
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr int workspace_size(const ArrowheadShape& shape,
                                                                       const int nrhs) noexcept {
        return valid(shape, nrhs) ? offsets(shape, nrhs).total : 0;
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr int
    streamed_workspace_size(const ArrowheadShape& shape) noexcept {
        return valid(shape) ? streamed_offsets(shape, 1).total : 0;
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr int
    streamed_workspace_size(const ArrowheadShape& shape, const int nrhs) noexcept {
        return valid(shape, nrhs) ? streamed_offsets(shape, nrhs).total : 0;
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr int
    pivot_size(const ArrowheadShape& shape) noexcept {
        if (!valid(shape)) return 0;
        return shape.block > shape.border ? shape.block : shape.border;
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
    solve_inplace(const ArrowheadShape& shape, T* TDLS_RESTRICT A, const int A_stride,
                  int* TDLS_RESTRICT piv, const int piv_stride, T* TDLS_RESTRICT work,
                  const int work_stride, T* TDLS_RESTRICT y, const int rhs_stride,
                  int& refused_count) noexcept {
        return solve_core<1, false>(shape, 1, A, A_stride, piv, piv_stride, work, work_stride, y,
                                    rhs_stride, 0, 0, refused_count);
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
    solve_inplace(const ArrowheadShape& shape, T* TDLS_RESTRICT A, const int A_stride,
                  int* TDLS_RESTRICT piv, const int piv_stride, T* TDLS_RESTRICT work,
                  const int work_stride, T* TDLS_RESTRICT y, const int rhs_stride) noexcept {
        int refused_count = 0;
        return solve_inplace(shape, A, A_stride, piv, piv_stride, work, work_stride, y, rhs_stride,
                             refused_count);
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
    solve(const ArrowheadShape& shape, T* TDLS_RESTRICT A, const int A_stride,
          int* TDLS_RESTRICT piv, const int piv_stride, T* TDLS_RESTRICT work,
          const int work_stride, const T* TDLS_RESTRICT b, T* TDLS_RESTRICT x, const int rhs_stride,
          int& refused_count) noexcept {
        if (!valid(shape)) {
            refused_count = 0;
            return false;
        }
        const int n = shape.n();
        for (int i = 0; i < n; ++i)
            x[TDLS_STRIDED(i, rhs_stride)] = b[TDLS_STRIDED(i, rhs_stride)];
        return solve_inplace(shape, A, A_stride, piv, piv_stride, work, work_stride, x, rhs_stride,
                             refused_count);
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
    solve(const ArrowheadShape& shape, T* TDLS_RESTRICT A, const int A_stride,
          int* TDLS_RESTRICT piv, const int piv_stride, T* TDLS_RESTRICT work,
          const int work_stride, const T* TDLS_RESTRICT b, T* TDLS_RESTRICT x,
          const int rhs_stride) noexcept {
        int refused_count = 0;
        return solve(shape, A, A_stride, piv, piv_stride, work, work_stride, b, x, rhs_stride,
                     refused_count);
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
    solve_inplace_multirhs(const ArrowheadShape& shape, const int nrhs, T* TDLS_RESTRICT A,
                           const int A_stride, int* TDLS_RESTRICT piv, const int piv_stride,
                           T* TDLS_RESTRICT work, const int work_stride, T* TDLS_RESTRICT y,
                           const int rhs_stride, const std::size_t xcol_stride,
                           int& refused_count) noexcept {
        return solve_core<0, false>(shape, nrhs, A, A_stride, piv, piv_stride, work, work_stride, y,
                                    rhs_stride, xcol_stride, 0, refused_count);
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
    solve_inplace_multirhs(const ArrowheadShape& shape, const int nrhs, T* TDLS_RESTRICT A,
                           const int A_stride, int* TDLS_RESTRICT piv, const int piv_stride,
                           T* TDLS_RESTRICT work, const int work_stride, T* TDLS_RESTRICT y,
                           const int rhs_stride, const std::size_t xcol_stride) noexcept {
        int refused_count = 0;
        return solve_inplace_multirhs(shape, nrhs, A, A_stride, piv, piv_stride, work, work_stride,
                                      y, rhs_stride, xcol_stride, refused_count);
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
    solve_multirhs(const ArrowheadShape& shape, const int nrhs, T* TDLS_RESTRICT A,
                   const int A_stride, int* TDLS_RESTRICT piv, const int piv_stride,
                   T* TDLS_RESTRICT work, const int work_stride, const T* TDLS_RESTRICT b,
                   T* TDLS_RESTRICT x, const int rhs_stride, const std::size_t xcol_stride,
                   int& refused_count) noexcept {
        if (!valid(shape, nrhs)) {
            refused_count = 0;
            return false;
        }
        const int n = shape.n();
        for (int m = 0; m < nrhs; ++m)
            for (int i = 0; i < n; ++i)
                x[TDLS_STRIDED(i, rhs_stride) + TDLS_STRIDED_WIDE(m, xcol_stride)] =
                    b[TDLS_STRIDED(i, rhs_stride) + TDLS_STRIDED_WIDE(m, xcol_stride)];
        return solve_inplace_multirhs(shape, nrhs, A, A_stride, piv, piv_stride, work, work_stride,
                                      x, rhs_stride, xcol_stride, refused_count);
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
    solve_multirhs(const ArrowheadShape& shape, const int nrhs, T* TDLS_RESTRICT A,
                   const int A_stride, int* TDLS_RESTRICT piv, const int piv_stride,
                   T* TDLS_RESTRICT work, const int work_stride, const T* TDLS_RESTRICT b,
                   T* TDLS_RESTRICT x, const int rhs_stride,
                   const std::size_t xcol_stride) noexcept {
        int refused_count = 0;
        return solve_multirhs(shape, nrhs, A, A_stride, piv, piv_stride, work, work_stride, b, x,
                              rhs_stride, xcol_stride, refused_count);
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool solve_border_canonical_multirhs(
        const ArrowheadShape& shape, const int nrhs, T* TDLS_RESTRICT A, const int A_stride,
        int* TDLS_RESTRICT piv, const int piv_stride, T* TDLS_RESTRICT work, const int work_stride,
        const int col0, T* TDLS_RESTRICT x, const int rhs_stride, const std::size_t xcol_stride,
        int& refused_count) noexcept {
        return solve_core<0, true>(shape, nrhs, A, A_stride, piv, piv_stride, work, work_stride, x,
                                   rhs_stride, xcol_stride, col0, refused_count);
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
    solve_border_canonical_multirhs(const ArrowheadShape& shape, const int nrhs, T* TDLS_RESTRICT A,
                                    const int A_stride, int* TDLS_RESTRICT piv,
                                    const int piv_stride, T* TDLS_RESTRICT work,
                                    const int work_stride, const int col0, T* TDLS_RESTRICT x,
                                    const int rhs_stride, const std::size_t xcol_stride) noexcept {
        int refused_count = 0;
        return solve_border_canonical_multirhs(shape, nrhs, A, A_stride, piv, piv_stride, work,
                                               work_stride, col0, x, rhs_stride, xcol_stride,
                                               refused_count);
    }
 
    struct StreamedBlock {
        T* Dl;      
        T* W;       
        T* Abk;     
        int S;      
        int B;      
        int stride; 
        TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T& d(const int i,
                                                              const int j) const noexcept {
            return Dl[TDLS_STRIDED(col ? j * S + i : i * S + j, stride)];
        }
        TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T& akb(const int i,
                                                                const int c) const noexcept {
            return W[TDLS_STRIDED(c * S + i, stride)];
        }
        TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T& y(const int i) const noexcept {
            return W[TDLS_STRIDED(B * S + i, stride)];
        }
        TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T& coupling(const int i) const noexcept {
            return Abk[TDLS_STRIDED(i, stride)];
        }
        TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T& abk(const int r,
                                                                const int i) const noexcept {
            return Abk[TDLS_STRIDED(r * S + i, stride)];
        }
    };
 
    struct RecomputedBlock : StreamedBlock {
        bool recovery = false;
    };
 
    enum class StreamedOutcome {
        solved,  
        tangent, 
        singular 
    };
 
    template<typename Provider>
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
    solve_inplace_streamed(const ArrowheadShape& shape, Provider&& provider, T* TDLS_RESTRICT A,
                           const int A_stride, int* TDLS_RESTRICT piv, const int piv_stride,
                           T* TDLS_RESTRICT work, const int work_stride, T* TDLS_RESTRICT y,
                           const int rhs_stride, int& refused_count) noexcept {
        return solve_core<1, false>(shape, 1, A, A_stride, piv, piv_stride, work, work_stride, y,
                                    rhs_stride, 0, 0, refused_count, &provider);
    }
 
    template<typename Provider>
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool solve_border_canonical_multirhs_streamed(
        const ArrowheadShape& shape, const int nrhs, Provider&& provider, T* TDLS_RESTRICT A,
        const int A_stride, int* TDLS_RESTRICT piv, const int piv_stride, T* TDLS_RESTRICT work,
        const int work_stride, const int col0, T* TDLS_RESTRICT x, const int rhs_stride,
        const std::size_t xcol_stride, int& refused_count) noexcept {
        return solve_core<0, true>(shape, nrhs, A, A_stride, piv, piv_stride, work, work_stride, x,
                                   rhs_stride, xcol_stride, col0, refused_count, &provider);
    }
 
    template<typename Provider, typename Converged>
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr StreamedOutcome
    solve_inplace_streamed_or_tangent(const ArrowheadShape& shape, const int ntan,
                                      Provider&& provider, Converged&& converged,
                                      T* TDLS_RESTRICT A, const int A_stride,
                                      int* TDLS_RESTRICT piv, const int piv_stride,
                                      T* TDLS_RESTRICT work, const int work_stride,
                                      T* TDLS_RESTRICT y, const int rhs_stride, const int col0,
                                      T* TDLS_RESTRICT x, const int x_stride,
                                      const std::size_t xcol_stride, int& refused_count) noexcept {
        if (!valid(shape, ntan)) {
            refused_count = 0;
            return StreamedOutcome::singular;
        }
        // The pass carries one right-hand side through the blocks (y, with
        // one column); the workspace holds both phases, ntan columns in the
        // front and in the border right-hand sides.
        [[maybe_unused]] constexpr int ycols = 1;
        constexpr int y_xcol_stride          = 0;
        const Sizes z                        = sizes(shape, 1, 1);
        const Offsets o                      = offsets(shape, ntan);
        const StreamedOffsets so             = streamed_offsets(shape, ntan);
        const int B = z.B, R = z.R, R0 = z.R0, NB = z.NB;
        if constexpr (!in_place) {
            for (int r = 0; r < R; ++r) {
                for (int c = 0; c < B; ++c)
                    TDLS_ARROW_LUPP_DYN_W(o.Sacc + r * B + c) =
                        TDLS_ARROW_LUPP_DYN_A(storage_map::abb(shape, R0 + r, c));
                TDLS_ARROW_LUPP_DYN_W(o.racc + r) = TDLS_ARROW_LUPP_DYN_Y(R0 + r);
            }
        }
        refused_count = 0;
        for (int k = 0; k < NB; ++k) {
            StreamedBlock blk = streamed_block(z, o, so, work, work_stride);
            provider(k, blk);
            const bool ok = condense_or_store_streamed<1, false>(shape, z, o, so, A, A_stride, piv,
                                                                 piv_stride, work, work_stride, y,
                                                                 rhs_stride, y_xcol_stride, k, 1);
            TDLS_ARROW_LUPP_DYN_W(o.refused + k) = ok ? T(0) : T(1);
            if (!ok) ++refused_count;
        }
        // The rows as stored and the accumulated rows in separate loops, as
        // in solve_core.
        const int acc0 = in_place ? B : R0;
        const int acc1 = in_place ? B : R0 + R;
        for (int r = 0; r < acc0; ++r)
            for (int c = 0; c < B; ++c)
                TDLS_ARROW_LUPP_DYN_W(o.Sbb + r * B + c) =
                    TDLS_ARROW_LUPP_DYN_A(storage_map::abb(shape, r, c));
        for (int r = acc0; r < acc1; ++r)
            for (int c = 0; c < B; ++c)
                TDLS_ARROW_LUPP_DYN_W(o.Sbb + r * B + c) =
                    TDLS_ARROW_LUPP_DYN_W(o.Sacc + (r - R0) * B + c);
        for (int r = acc1; r < B; ++r)
            for (int c = 0; c < B; ++c)
                TDLS_ARROW_LUPP_DYN_W(o.Sbb + r * B + c) =
                    TDLS_ARROW_LUPP_DYN_A(storage_map::abb(shape, r, c));
        T* const Sbb = work + TDLS_STRIDED(o.Sbb, work_stride);
        T* const rb  = work + TDLS_STRIDED(o.rb, work_stride);
        if (converged()) {
            const Sizes zt = sizes(shape, ntan, 0);
            for (int m = 0; m < ntan; ++m)
                for (int r = 0; r < B; ++r)
                    TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + r) = r == col0 + m ? T(1) : T(0);
            // The refused blocks, stored as given, by the second pass of the
            // canonical solve, which writes neither x nor y.
            if (refused_count > 0 &&
                !second_pass<0, true>(shape, zt, o, A, A_stride, work, work_stride, y, rhs_stride,
                                      y_xcol_stride, ntan))
                return StreamedOutcome::singular;
            const bool ok = ntan == 1 ? BorderSolver::solve_inplace(B, Sbb, work_stride, piv,
                                                                    piv_stride, rb, work_stride)
                                      : BorderSolver::solve_inplace_multirhs(
                                            B, ntan, Sbb, work_stride, piv, piv_stride, rb,
                                            work_stride, TDLS_STRIDED(B, work_stride));
            if (!ok) return StreamedOutcome::singular;
            for (int m = 0; m < ntan; ++m)
                for (int r = 0; r < B; ++r)
                    x[TDLS_STRIDED(r, x_stride) + TDLS_STRIDED_WIDE(m, xcol_stride)] =
                        TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + r);
            return StreamedOutcome::tangent;
        }
        for (int r = 0; r < acc0; ++r)
            TDLS_ARROW_LUPP_DYN_W(o.rb + r) = TDLS_ARROW_LUPP_DYN_Y(r);
        for (int r = acc0; r < acc1; ++r)
            TDLS_ARROW_LUPP_DYN_W(o.rb + r) = TDLS_ARROW_LUPP_DYN_W(o.racc + r - R0);
        for (int r = acc1; r < B; ++r)
            TDLS_ARROW_LUPP_DYN_W(o.rb + r) = TDLS_ARROW_LUPP_DYN_Y(r);
        if (refused_count > 0 && !second_pass<1, false>(shape, z, o, A, A_stride, work, work_stride,
                                                        y, rhs_stride, y_xcol_stride, 1))
            return StreamedOutcome::singular;
        if (!BorderSolver::solve_inplace(B, Sbb, work_stride, piv, piv_stride, rb, work_stride))
            return StreamedOutcome::singular;
        for (int r = 0; r < B; ++r)
            TDLS_ARROW_LUPP_DYN_Y(r) = TDLS_ARROW_LUPP_DYN_W(o.rb + r);
        recover<1>(shape, z, o, A, A_stride, work, work_stride, y, rhs_stride, y_xcol_stride, 1);
        return StreamedOutcome::solved;
    }
 
    enum class RecomputedOutcome {
        solved,   
        tangent,  
        singular, 
        deferred  
    };
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
    valid_recomputed(const ArrowheadShape& shape, const int capacity, const int nrhs = 1) noexcept {
        constexpr int digits =
            std::numeric_limits<T>::digits < 62 ? std::numeric_limits<T>::digits : 62;
        return valid(shape, nrhs) && capacity >= 0 && capacity <= shape.blocks &&
               static_cast<long long>(shape.blocks) <= (1LL << digits) &&
               recomputed_offsets<long long>(shape, capacity, nrhs).total < (1LL << 31);
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr int
    recomputed_workspace_size(const ArrowheadShape& shape, const int capacity,
                              const int nrhs = 1) noexcept {
        return valid_recomputed(shape, capacity, nrhs)
                   ? recomputed_offsets(shape, capacity, nrhs).total
                   : 0;
    }
 
    template<typename Provider>
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr RecomputedOutcome
    solve_inplace_recomputed(const ArrowheadShape& shape, const int capacity, Provider&& provider,
                             const T* TDLS_RESTRICT Abb, const int Abb_stride,
                             int* TDLS_RESTRICT piv, const int piv_stride, T* TDLS_RESTRICT work,
                             const int work_stride, T* TDLS_RESTRICT y, const int rhs_stride,
                             int& refused_count) noexcept {
        return recomputed_core(shape, capacity, 0, provider, Abb, Abb_stride, piv, piv_stride, work,
                               work_stride, y, rhs_stride, 0, nullptr, 1, 0, refused_count);
    }
 
    template<typename Provider, typename Converged>
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr RecomputedOutcome
    solve_inplace_recomputed_or_tangent(
        const ArrowheadShape& shape, const int ntan, const int capacity, Provider&& provider,
        Converged&& converged, const T* TDLS_RESTRICT Abb, const int Abb_stride,
        int* TDLS_RESTRICT piv, const int piv_stride, T* TDLS_RESTRICT work, const int work_stride,
        T* TDLS_RESTRICT y, const int rhs_stride, const int col0, T* TDLS_RESTRICT x,
        const int x_stride, const std::size_t xcol_stride, int& refused_count) noexcept {
        if (ntan < 1) {
            refused_count = 0;
            return RecomputedOutcome::singular;
        }
        return recomputed_core(shape, capacity, ntan, provider, Abb, Abb_stride, piv, piv_stride,
                               work, work_stride, y, rhs_stride, col0, x, x_stride, xcol_stride,
                               refused_count, &converged);
    }
 
  private:
    template<typename I>
    struct RecomputedOffsetsOf {
        I store, ystore, slots, total;
    };
    using RecomputedOffsets = RecomputedOffsetsOf<int>;
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr ArrowheadShape
    store_shape(const ArrowheadShape& shape, const int capacity) noexcept {
        return {shape.border, shape.block, capacity, shape.coupling_row};
    }
 
    template<typename I = int>
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr RecomputedOffsetsOf<I>
    recomputed_offsets(const ArrowheadShape& shape, const int capacity, const int nc) noexcept {
        const ArrowheadShape st = store_shape(shape, capacity);
        RecomputedOffsetsOf<I> ro{};
        I e      = streamed_offsets<I>(st, nc).total;
        ro.store = e;
        e += capacity > 0 ? I(storage_map::size(st)) : I(0);
        ro.ystore = e;
        e += I(shape.border) + I(capacity) * shape.block;
        ro.slots = e;
        e += capacity;
        ro.total = e;
        return ro;
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE static constexpr T border_entry(const T* TDLS_RESTRICT Abb,
                                                                   const int Abb_stride,
                                                                   const int B, const int r,
                                                                   const int c) noexcept {
        return Abb[TDLS_STRIDED(col ? c * B + r : r * B + c, Abb_stride)];
    }
 
    template<typename Provider, typename... Converged>
    TDLS_HOST_DEVICE static constexpr RecomputedOutcome recomputed_core(
        const ArrowheadShape& shape, const int capacity, const int ntan, Provider& provider,
        const T* TDLS_RESTRICT Abb, const int Abb_stride, int* TDLS_RESTRICT piv,
        const int piv_stride, T* TDLS_RESTRICT work, const int work_stride, T* TDLS_RESTRICT y,
        const int rhs_stride, [[maybe_unused]] const int col0, [[maybe_unused]] T* TDLS_RESTRICT x,
        [[maybe_unused]] const int x_stride, [[maybe_unused]] const std::size_t xcol_stride,
        int& refused_count, Converged*... converged) noexcept {
        static_assert(!in_place, "ArrowheadLUppSolverDynamic: the recomputed solve accumulates "
                                 "the Schur complement apart (schur_in_place off)");
        // Columns of the front: the tangent columns, one at least.
        const int fnc = ntan > 1 ? ntan : 1;
        refused_count = 0;
        if (!valid_recomputed(shape, capacity, fnc) || ntan > shape.border)
            return RecomputedOutcome::singular;
        [[maybe_unused]] constexpr int ycols = 1;
        const ArrowheadShape st              = store_shape(shape, capacity);
        const Sizes z                        = sizes(shape, 1, 1);
        const Offsets o                      = offsets(st, fnc);
        const StreamedOffsets so             = streamed_offsets(st, fnc);
        const RecomputedOffsets ro           = recomputed_offsets(shape, capacity, fnc);
        const int B = z.B, S = z.S, R = z.R, R0 = z.R0, NB = z.NB;
        for (int r = 0; r < R; ++r) {
            for (int c = 0; c < B; ++c)
                TDLS_ARROW_LUPP_DYN_W(o.Sacc + r * B + c) =
                    border_entry(Abb, Abb_stride, B, R0 + r, c);
            TDLS_ARROW_LUPP_DYN_W(o.racc + r) = TDLS_ARROW_LUPP_DYN_Y(R0 + r);
        }
        int count = 0;
        for (int k = 0; k < NB; ++k) {
            RecomputedBlock blk{streamed_block(z, o, so, work, work_stride), false};
            provider(k, blk);
            const int verdict = condense_or_keep_recomputed(
                shape, z, o, so, ro, capacity, piv, piv_stride, work, work_stride, k, count);
            refused_count = count;
            if (verdict == 2) return RecomputedOutcome::deferred;
        }
        // The flags of the places of the store, for its second pass.
        for (int j = 0; j < capacity; ++j)
            TDLS_ARROW_LUPP_DYN_W(o.refused + j) = j < count ? T(1) : T(0);
        T* const store  = work + TDLS_STRIDED(ro.store, work_stride);
        T* const ystore = work + TDLS_STRIDED(ro.ystore, work_stride);
        for (int r = 0; r < B; ++r) {
            const bool acc = r >= R0 && r < R0 + R;
            for (int c = 0; c < B; ++c)
                TDLS_ARROW_LUPP_DYN_W(o.Sbb + r * B + c) =
                    acc ? TDLS_ARROW_LUPP_DYN_W(o.Sacc + (r - R0) * B + c)
                        : border_entry(Abb, Abb_stride, B, r, c);
        }
        T* const Sbb = work + TDLS_STRIDED(o.Sbb, work_stride);
        T* const rb  = work + TDLS_STRIDED(o.rb, work_stride);
        if constexpr (sizeof...(Converged) > 0) {
            if ((*streamed_arg(converged...))()) {
                // The tangent: unit columns, the refused blocks of the store
                // by the second pass of the canonical solve.
                const Sizes zt = sizes(st, fnc, 0);
                for (int m = 0; m < fnc; ++m)
                    for (int r = 0; r < B; ++r)
                        TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + r) = r == col0 + m ? T(1) : T(0);
                if (count > 0 && !second_pass<0, true>(st, zt, o, store, work_stride, work,
                                                       work_stride, ystore, work_stride, 0, fnc))
                    return RecomputedOutcome::singular;
                const bool ok = fnc == 1 ? BorderSolver::solve_inplace(B, Sbb, work_stride, piv,
                                                                       piv_stride, rb, work_stride)
                                         : BorderSolver::solve_inplace_multirhs(
                                               B, fnc, Sbb, work_stride, piv, piv_stride, rb,
                                               work_stride, TDLS_STRIDED(B, work_stride));
                if (!ok) return RecomputedOutcome::singular;
                for (int m = 0; m < fnc; ++m)
                    for (int r = 0; r < B; ++r)
                        x[TDLS_STRIDED(r, x_stride) + TDLS_STRIDED_WIDE(m, xcol_stride)] =
                            TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + r);
                return RecomputedOutcome::tangent;
            }
        }
        for (int r = 0; r < B; ++r) {
            const bool acc = r >= R0 && r < R0 + R;
            TDLS_ARROW_LUPP_DYN_W(o.rb + r) =
                acc ? TDLS_ARROW_LUPP_DYN_W(o.racc + r - R0) : TDLS_ARROW_LUPP_DYN_Y(r);
        }
        if (count > 0 && !second_pass<1, false>(st, sizes(st, 1, 1), o, store, work_stride, work,
                                                work_stride, ystore, work_stride, 0, 1))
            return RecomputedOutcome::singular;
        if (!BorderSolver::solve_inplace(B, Sbb, work_stride, piv, piv_stride, rb, work_stride))
            return RecomputedOutcome::singular;
        // The accepted blocks, recomputed, the last one first: the places
        // of the store hold the refused ones in increasing order.
        int j = count - 1;
        for (int k = NB - 1; k >= 0; --k) {
            if (j >= 0 && TDLS_ARROW_LUPP_DYN_W(ro.slots + j) == T(k)) {
                --j;
                continue;
            }
            RecomputedBlock blk{streamed_block(z, o, so, work, work_stride), true};
            provider(k, blk);
            // y_k - A_kb x_b, then the column solved in place of A_kb(:, 0)
            for (int i = 0; i < S; ++i) {
                T v = TDLS_ARROW_LUPP_DYN_W(o.W + B * S + i);
                for (int c = 0; c < B; ++c)
                    v -= TDLS_ARROW_LUPP_DYN_W(o.W + c * S + i) * TDLS_ARROW_LUPP_DYN_W(o.rb + c);
                TDLS_ARROW_LUPP_DYN_W(o.s + i) = v;
            }
            for (int i = 0; i < S; ++i)
                TDLS_ARROW_LUPP_DYN_W(o.W + i) = TDLS_ARROW_LUPP_DYN_W(o.s + i);
            if (z.reg) {
                (void)physical_factor(z, o, work, work_stride, piv, piv_stride);
                physical_solve(z, o, work, work_stride, piv, piv_stride, 1, true);
            } else {
                T* const Dl = work + TDLS_STRIDED(o.Dl, work_stride);
                (void)BlockSolver::factorize(S, Dl, work_stride, piv, piv_stride);
                BlockSolver::substitute_inplace(S, Dl, work_stride, piv, piv_stride,
                                                work + TDLS_STRIDED(o.W, work_stride), work_stride);
            }
            for (int i = 0; i < S; ++i)
                TDLS_ARROW_LUPP_DYN_Y(B + k * S + i) = TDLS_ARROW_LUPP_DYN_W(o.W + i);
        }
        if (count > 0)
            recover_refused(shape, st, z, o, ro, work, work_stride, y, rhs_stride, count);
        for (int r = 0; r < B; ++r)
            TDLS_ARROW_LUPP_DYN_Y(r) = TDLS_ARROW_LUPP_DYN_W(o.rb + r);
        return RecomputedOutcome::solved;
    }
 
    TDLS_HOST_DEVICE static constexpr int
    condense_or_keep_recomputed(const ArrowheadShape& shape, const Sizes& z, const Offsets& o,
                                const StreamedOffsets& so, const RecomputedOffsets& ro,
                                const int capacity, int* TDLS_RESTRICT piv, const int piv_stride,
                                T* TDLS_RESTRICT work, const int work_stride, const int k,
                                int& count) noexcept {
        const int B = z.B, S = z.S;
        const bool forward_check = !z.reg && col && z.pass_width == B + 1;
        for (int e = 0; e < S * S; ++e)
            TDLS_ARROW_LUPP_DYN_W(so.Dsave + e) = TDLS_ARROW_LUPP_DYN_W(o.Dl + e);
        if (forward_check)
            for (int e = 0; e < (B + 1) * S; ++e)
                TDLS_ARROW_LUPP_DYN_W(so.Wsave + e) = TDLS_ARROW_LUPP_DYN_W(o.W + e);
        if (condense<1, false, true>(shape, z, o, nullptr, 1, piv, piv_stride, work, work_stride,
                                     nullptr, 1, 0, k, 1, &so))
            return 0;
        ++count;
        if (count > capacity) return 2;
        // The block as given, in place j of the store (W is intact unless the
        // forward stage ran).
        const int j             = count - 1;
        const ArrowheadShape st = store_shape(shape, capacity);
        const int ow            = forward_check ? so.Wsave : o.W;
        for (int i = 0; i < S; ++i) {
            for (int jj = 0; jj < S; ++jj)
                TDLS_ARROW_LUPP_DYN_W(ro.store + storage_map::d(st, j, i, jj)) =
                    TDLS_ARROW_LUPP_DYN_W(so.Dsave + (col ? jj * S + i : i * S + jj));
            for (int c = 0; c < B; ++c)
                TDLS_ARROW_LUPP_DYN_W(ro.store + storage_map::akb(st, j, i, c)) =
                    TDLS_ARROW_LUPP_DYN_W(ow + c * S + i);
            TDLS_ARROW_LUPP_DYN_W(ro.ystore + B + j * S + i) =
                TDLS_ARROW_LUPP_DYN_W(ow + B * S + i);
        }
        if constexpr (diag) {
            for (int i = 0; i < S; ++i)
                TDLS_ARROW_LUPP_DYN_W(ro.store + storage_map::abk_diag(st, j, i)) =
                    TDLS_ARROW_LUPP_DYN_W(o.acoef + i);
        } else {
            for (int r = 0; r < B; ++r)
                for (int i = 0; i < S; ++i)
                    TDLS_ARROW_LUPP_DYN_W(ro.store + storage_map::abk(st, j, r, i)) =
                        TDLS_ARROW_LUPP_DYN_W(so.Abk + r * S + i);
        }
        TDLS_ARROW_LUPP_DYN_W(ro.slots + j) = T(k);
        return 1;
    }
 
    TDLS_HOST_DEVICE static constexpr void
    recover_refused(const ArrowheadShape& shape, const ArrowheadShape& st, const Sizes& z,
                    const Offsets& o, const RecomputedOffsets& ro, T* TDLS_RESTRICT work,
                    const int work_stride, T* TDLS_RESTRICT y, const int rhs_stride,
                    const int count) noexcept {
        (void)shape;
        const int B = z.B, S = z.S, R = z.R;
        for (int a = 0; a < R; ++a)
            TDLS_ARROW_LUPP_DYN_W(o.w + a) = T(0);
        for (int j = count - 1; j >= 0; --j) {
            const int k = static_cast<int>(TDLS_ARROW_LUPP_DYN_W(ro.slots + j));
            for (int i = 0; i < S; ++i) {
                T v = TDLS_ARROW_LUPP_DYN_W(ro.ystore + B + j * S + i);
                for (int c = 0; c < B; ++c)
                    v -= TDLS_ARROW_LUPP_DYN_W(ro.store + storage_map::akb(st, j, i, c)) *
                         TDLS_ARROW_LUPP_DYN_W(o.rb + c);
                for (int a = 0; a < R; ++a) {
                    const T h =
                        a < S ? TDLS_ARROW_LUPP_DYN_W(ro.store + storage_map::d(st, j, i, a))
                              : TDLS_ARROW_LUPP_DYN_W(o.extra + (j * S + i) * z.extra_cols + a - S);
                    v -= h * TDLS_ARROW_LUPP_DYN_W(o.w + a);
                }
                TDLS_ARROW_LUPP_DYN_W(o.s + i) = v;
            }
            if constexpr (diag) {
                for (int i = 0; i < S; ++i)
                    TDLS_ARROW_LUPP_DYN_W(o.w + i) +=
                        TDLS_ARROW_LUPP_DYN_W(ro.store + storage_map::abk_diag(st, j, i)) *
                        TDLS_ARROW_LUPP_DYN_W(o.s + i);
            } else {
                for (int r = 0; r < B; ++r)
                    for (int i = 0; i < S; ++i)
                        TDLS_ARROW_LUPP_DYN_W(o.w + r) +=
                            TDLS_ARROW_LUPP_DYN_W(ro.store + storage_map::abk(st, j, r, i)) *
                            TDLS_ARROW_LUPP_DYN_W(o.s + i);
            }
            for (int i = 0; i < S; ++i)
                TDLS_ARROW_LUPP_DYN_Y(B + k * S + i) = TDLS_ARROW_LUPP_DYN_W(o.s + i);
        }
    }
 
  private:
    template<typename P>
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL static constexpr P*
    streamed_arg(P* p) noexcept {
        return p;
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL static constexpr StreamedBlock
    streamed_block(const Sizes& z, const Offsets& o, const StreamedOffsets& so,
                   T* TDLS_RESTRICT work, const int work_stride) noexcept {
        StreamedBlock blk{};
        blk.Dl     = work + TDLS_STRIDED(o.Dl, work_stride);
        blk.W      = work + TDLS_STRIDED(o.W, work_stride);
        blk.Abk    = work + TDLS_STRIDED(diag ? o.acoef : so.Abk, work_stride);
        blk.S      = z.S;
        blk.B      = z.B;
        blk.stride = work_stride;
        return blk;
    }
 
    template<int NC, bool canonical, typename... Provider>
    TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL static constexpr bool
    solve_core(const ArrowheadShape& shape, const int nrhs, T* TDLS_RESTRICT A, const int A_stride,
               int* TDLS_RESTRICT piv, const int piv_stride, T* TDLS_RESTRICT work,
               const int work_stride, T* TDLS_RESTRICT y, const int rhs_stride,
               [[maybe_unused]] const std::size_t xcol_stride, [[maybe_unused]] const int col0,
               int& refused_count, Provider*... provider) noexcept {
        const int nc = NC > 0 ? NC : nrhs;
        if (!valid(shape, nc)) {
            refused_count = 0;
            return false;
        }
        const int ny = canonical ? 0 : nc;
        // Columns of y when fixed, for TDLS_ARROW_LUPP_DYN_YW.
        [[maybe_unused]] constexpr int ycols = NC;
        const Sizes z                        = sizes(shape, nc, ny);
        const Offsets o                      = offsets(shape, nc);
        const int B = z.B, R = z.R, R0 = z.R0, NB = z.NB;
        if constexpr (!in_place) {
            for (int r = 0; r < R; ++r) {
                for (int c = 0; c < B; ++c)
                    TDLS_ARROW_LUPP_DYN_W(o.Sacc + r * B + c) =
                        TDLS_ARROW_LUPP_DYN_A(storage_map::abb(shape, R0 + r, c));
                for (int m = 0; m < ny; ++m)
                    TDLS_ARROW_LUPP_DYN_W(o.racc + m * R + r) = TDLS_ARROW_LUPP_DYN_YW(m, R0 + r);
            }
        }
        refused_count = 0;
        if constexpr (sizeof...(Provider) > 0) {
            // Streamed: the provider writes block k in the workspace.
            const StreamedOffsets so = streamed_offsets(shape, nc);
            for (int k = 0; k < NB; ++k) {
                StreamedBlock blk = streamed_block(z, o, so, work, work_stride);
                (*streamed_arg(provider...))(k, blk);
                const bool ok = condense_or_store_streamed<NC, canonical>(
                    shape, z, o, so, A, A_stride, piv, piv_stride, work, work_stride, y, rhs_stride,
                    xcol_stride, k, nc);
                TDLS_ARROW_LUPP_DYN_W(o.refused + k) = ok ? T(0) : T(1);
                if (!ok) ++refused_count;
            }
        } else {
            for (int k = 0; k < NB; ++k) {
                const bool ok =
                    condense<NC, canonical>(shape, z, o, A, A_stride, piv, piv_stride, work,
                                            work_stride, y, rhs_stride, xcol_stride, k, nc);
                TDLS_ARROW_LUPP_DYN_W(o.refused + k) = ok ? T(0) : T(1);
                if (!ok) ++refused_count;
            }
        }
        // The border system: the accumulated rows, the others as stored;
        // column m of its right-hand sides at m * B.
        // The rows as stored and the accumulated rows in separate loops:
        // GCC 11 and 12 vectorize a choice of source per entry wrongly (-O3
        // with AVX2 or AVX-512).
        const int acc0 = in_place ? B : R0;
        const int acc1 = in_place ? B : R0 + R;
        for (int r = 0; r < acc0; ++r) {
            for (int c = 0; c < B; ++c)
                TDLS_ARROW_LUPP_DYN_W(o.Sbb + r * B + c) =
                    TDLS_ARROW_LUPP_DYN_A(storage_map::abb(shape, r, c));
            for (int m = 0; m < nc; ++m) {
                if constexpr (canonical)
                    TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + r) = r == col0 + m ? T(1) : T(0);
                else
                    TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + r) = TDLS_ARROW_LUPP_DYN_YW(m, r);
            }
        }
        for (int r = acc0; r < acc1; ++r) {
            for (int c = 0; c < B; ++c)
                TDLS_ARROW_LUPP_DYN_W(o.Sbb + r * B + c) =
                    TDLS_ARROW_LUPP_DYN_W(o.Sacc + (r - R0) * B + c);
            for (int m = 0; m < nc; ++m) {
                if constexpr (canonical)
                    TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + r) = r == col0 + m ? T(1) : T(0);
                else
                    TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + r) =
                        TDLS_ARROW_LUPP_DYN_W(o.racc + m * R + r - R0);
            }
        }
        for (int r = acc1; r < B; ++r) {
            for (int c = 0; c < B; ++c)
                TDLS_ARROW_LUPP_DYN_W(o.Sbb + r * B + c) =
                    TDLS_ARROW_LUPP_DYN_A(storage_map::abb(shape, r, c));
            for (int m = 0; m < nc; ++m) {
                if constexpr (canonical)
                    TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + r) = r == col0 + m ? T(1) : T(0);
                else
                    TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + r) = TDLS_ARROW_LUPP_DYN_YW(m, r);
            }
        }
        if (refused_count > 0 &&
            !second_pass<NC, canonical>(shape, z, o, A, A_stride, work, work_stride, y, rhs_stride,
                                        xcol_stride, nc))
            return false;
        if (nc == 1) {
            if (!BorderSolver::solve_inplace(B, work + TDLS_STRIDED(o.Sbb, work_stride),
                                             work_stride, piv, piv_stride,
                                             work + TDLS_STRIDED(o.rb, work_stride), work_stride))
                return false;
        } else {
            if (!BorderSolver::solve_inplace_multirhs(
                    B, nc, work + TDLS_STRIDED(o.Sbb, work_stride), work_stride, piv, piv_stride,
                    work + TDLS_STRIDED(o.rb, work_stride), work_stride,
                    TDLS_STRIDED(B, work_stride)))
                return false;
        }
        for (int m = 0; m < nc; ++m)
            for (int r = 0; r < B; ++r) {
                if constexpr (canonical)
                    y[TDLS_STRIDED(r, rhs_stride) + TDLS_STRIDED_WIDE(m, xcol_stride)] =
                        TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + r);
                else
                    TDLS_ARROW_LUPP_DYN_YW(m, r) = TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + r);
            }
        if constexpr (!canonical)
            recover<NC>(shape, z, o, A, A_stride, work, work_stride, y, rhs_stride, xcol_stride,
                        nc);
        return true;
    }
 
    template<int NC, bool canonical, bool recompute = false, typename... Streamed>
    TDLS_HOST_DEVICE static constexpr bool
    condense(const ArrowheadShape& shape, const Sizes& z, const Offsets& o, T* TDLS_RESTRICT A,
             const int A_stride, int* TDLS_RESTRICT piv, const int piv_stride,
             T* TDLS_RESTRICT work, const int work_stride, T* TDLS_RESTRICT y, const int rhs_stride,
             [[maybe_unused]] const std::size_t xcol_stride, const int k,
             [[maybe_unused]] const int nrhs, const Streamed*... so) noexcept {
        constexpr bool streamed = sizeof...(Streamed) > 0;
        const int B = z.B, S = z.S;
        // Right-hand sides carried through the blocks.
        const int ny                         = canonical ? 0 : NC > 0 ? NC : nrhs;
        [[maybe_unused]] constexpr int ycols = NC;
        const T u                            = condensation_pivot_threshold;
        T* const Dl                          = work + TDLS_STRIDED(o.Dl, work_stride);
        T* const W                           = work + TDLS_STRIDED(o.W, work_stride);
        // The block, in the order of its storage (a streamed block is there
        // already).
        if constexpr (!streamed) {
            for (int i = 0; i < S; ++i)
                for (int j = 0; j < S; ++j)
                    TDLS_ARROW_LUPP_DYN_W(o.Dl + (col ? j * S + i : i * S + j)) =
                        TDLS_ARROW_LUPP_DYN_A(storage_map::d(shape, k, i, j));
        }
        // Factored by TiledLUpp, or by physical interchanges
        // (blocks in registers): Dl then holds U in logical row order and the
        // check of L is made.
        bool ok = true;
        if (z.reg) {
            const int verdict = physical_factor(z, o, work, work_stride, piv, piv_stride);
            if (verdict == 2) return false;
            ok = verdict == 0;
        } else {
            if (!BlockSolver::factorize(S, Dl, work_stride, piv, piv_stride)) return false;
        }
        if constexpr (diag && !streamed) {
            for (int i = 0; i < S; ++i)
                TDLS_ARROW_LUPP_DYN_W(o.acoef + i) =
                    TDLS_ARROW_LUPP_DYN_A(storage_map::abk_diag(shape, k, i));
        }
        const bool forward_check = !z.reg && col && z.pass_width == B + ny;
        // The multipliers L, and U in logical row order unless the forward
        // stage of the solve gives the border multipliers.
        if (z.reg) {
            for (int l = 0; l < S; ++l)
                for (int j = l; j < S; ++j)
                    TDLS_ARROW_LUPP_DYN_W(o.U + l * S + j) =
                        TDLS_ARROW_LUPP_DYN_W(o.Dl + l * S + j);
        } else {
            for (int l = 0; l < S; ++l) {
                const int p = piv[TDLS_STRIDED(l, piv_stride)];
                for (int j = 0; j < S; ++j) {
                    if (j < l) {
                        ok = ok && u * detail::abs(TDLS_ARROW_LUPP_DYN_W(o.Dl + p * S + j)) <= T(1);
                    } else if (!forward_check) {
                        TDLS_ARROW_LUPP_DYN_W(o.U + l * S + j) =
                            TDLS_ARROW_LUPP_DYN_W(o.Dl + p * S + j);
                    }
                }
            }
        }
        if (forward_check) {
            if constexpr (!streamed) {
                for (int i = 0; i < S; ++i) {
                    for (int c = 0; c < B; ++c)
                        TDLS_ARROW_LUPP_DYN_W(o.W + c * S + i) =
                            TDLS_ARROW_LUPP_DYN_A(storage_map::akb(shape, k, i, c));
                    for (int m = 0; m < ny; ++m)
                        TDLS_ARROW_LUPP_DYN_W(o.W + (B + m) * S + i) =
                            TDLS_ARROW_LUPP_DYN_YW(m, B + k * S + i);
                }
            }
            BlockSolver::forward_inplace_multirhs(S, B + ny, Dl, work_stride, piv, piv_stride, W,
                                                  work_stride, TDLS_STRIDED(S, work_stride));
            // Column c of U^-T A_kb is row c of the border multipliers.
            for (int e = 0; e < B * S; ++e)
                ok = ok && u * detail::abs(TDLS_ARROW_LUPP_DYN_W(o.W + e)) <= T(1);
            if (!ok) return false;
            BlockSolver::backward_inplace_multirhs(S, B + ny, Dl, work_stride, piv, piv_stride, W,
                                                   work_stride, TDLS_STRIDED(S, work_stride));
            condense_pass<NC, canonical, recompute>(shape, z, o, A, A_stride, piv, piv_stride, work,
                                                    work_stride, y, rhs_stride, xcol_stride, k, 0,
                                                    B + ny, true, so...);
        } else {
            ok = border_multipliers(shape, z, o, A, A_stride, work, work_stride, k, so...) && ok;
            if (!ok) return false;
            // U, recovery data: not for the canonical solve.
            if constexpr (multiplier_update && !canonical && !recompute) {
                for (int l = 0; l < S; ++l)
                    for (int j = l; j < S; ++j)
                        TDLS_ARROW_LUPP_DYN_A(storage_map::d(shape, k, l, j)) =
                            TDLS_ARROW_LUPP_DYN_W(o.U + l * S + j);
            }
            for (int c0 = 0; c0 < B + ny; c0 += z.pass_width) {
                const int w = B + ny - c0 < z.pass_width ? B + ny - c0 : z.pass_width;
                condense_pass<NC, canonical, recompute>(shape, z, o, A, A_stride, piv, piv_stride,
                                                        work, work_stride, y, rhs_stride,
                                                        xcol_stride, k, c0, w, false, so...);
            }
        }
        return true;
    }
 
    template<int NC, bool canonical>
    TDLS_HOST_DEVICE static constexpr bool condense_or_store_streamed(
        const ArrowheadShape& shape, const Sizes& z, const Offsets& o, const StreamedOffsets& so,
        T* TDLS_RESTRICT A, const int A_stride, int* TDLS_RESTRICT piv, const int piv_stride,
        T* TDLS_RESTRICT work, const int work_stride, T* TDLS_RESTRICT y, const int rhs_stride,
        [[maybe_unused]] const std::size_t xcol_stride, const int k, const int nrhs) noexcept {
        const int B = z.B, S = z.S;
        const int ny                         = canonical ? 0 : NC > 0 ? NC : nrhs;
        [[maybe_unused]] constexpr int ycols = NC;
        const bool forward_check             = !z.reg && col && z.pass_width == B + ny;
        for (int e = 0; e < S * S; ++e)
            TDLS_ARROW_LUPP_DYN_W(so.Dsave + e) = TDLS_ARROW_LUPP_DYN_W(o.Dl + e);
        if (forward_check)
            for (int e = 0; e < (B + ny) * S; ++e)
                TDLS_ARROW_LUPP_DYN_W(so.Wsave + e) = TDLS_ARROW_LUPP_DYN_W(o.W + e);
        if (condense<NC, canonical>(shape, z, o, A, A_stride, piv, piv_stride, work, work_stride, y,
                                    rhs_stride, xcol_stride, k, nrhs, &so))
            return true;
        // Refused: the block as given, where the stored path leaves a
        // refused block (W is intact unless the forward stage ran).
        const int ow = forward_check ? so.Wsave : o.W;
        for (int i = 0; i < S; ++i) {
            for (int j = 0; j < S; ++j)
                TDLS_ARROW_LUPP_DYN_A(storage_map::d(shape, k, i, j)) =
                    TDLS_ARROW_LUPP_DYN_W(so.Dsave + (col ? j * S + i : i * S + j));
            for (int c = 0; c < B; ++c)
                TDLS_ARROW_LUPP_DYN_A(storage_map::akb(shape, k, i, c)) =
                    TDLS_ARROW_LUPP_DYN_W(ow + c * S + i);
            for (int m = 0; m < ny; ++m)
                TDLS_ARROW_LUPP_DYN_YW(m, B + k * S + i) = TDLS_ARROW_LUPP_DYN_W(ow + B * S + i);
        }
        if constexpr (diag) {
            for (int i = 0; i < S; ++i)
                TDLS_ARROW_LUPP_DYN_A(storage_map::abk_diag(shape, k, i)) =
                    TDLS_ARROW_LUPP_DYN_W(o.acoef + i);
        } else {
            for (int r = 0; r < B; ++r)
                for (int i = 0; i < S; ++i)
                    TDLS_ARROW_LUPP_DYN_A(storage_map::abk(shape, k, r, i)) =
                        TDLS_ARROW_LUPP_DYN_W(so.Abk + r * S + i);
        }
        return false;
    }
 
    template<typename... Streamed>
    TDLS_HOST_DEVICE static constexpr bool
    border_multipliers(const ArrowheadShape& shape, const Sizes& z, const Offsets& o,
                       const T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT work,
                       const int work_stride, const int k, const Streamed*... so) noexcept {
        constexpr bool streamed = sizeof...(Streamed) > 0;
        const int B = z.B, S = z.S;
        bool ok = true;
        if constexpr (col) {
            for (int c = 0; c < B; ++c) {
                for (int i = 0; i < S; ++i) {
                    if constexpr (streamed)
                        TDLS_ARROW_LUPP_DYN_W(o.m + i) = TDLS_ARROW_LUPP_DYN_W(o.W + c * S + i);
                    else
                        TDLS_ARROW_LUPP_DYN_W(o.m + i) =
                            TDLS_ARROW_LUPP_DYN_A(storage_map::akb(shape, k, i, c));
                }
                ok = multipliers_pass(z, o, work, work_stride, o.m, 0) && ok;
            }
        } else if constexpr (diag) {
            // Border row coupling_row + i holds one entry, in column i.
            for (int i = 0; i < S; ++i) {
                for (int j = 0; j < S; ++j)
                    TDLS_ARROW_LUPP_DYN_W(o.m + j) =
                        j == i ? TDLS_ARROW_LUPP_DYN_W(o.acoef + i) : T(0);
                ok = multipliers_pass(z, o, work, work_stride, o.m, i) && ok;
            }
        } else {
            for (int r = 0; r < B; ++r) {
                for (int j = 0; j < S; ++j) {
                    if constexpr (streamed)
                        TDLS_ARROW_LUPP_DYN_W(o.Mb + r * S + j) =
                            TDLS_ARROW_LUPP_DYN_W(streamed_arg(so...)->Abk + r * S + j);
                    else
                        TDLS_ARROW_LUPP_DYN_W(o.Mb + r * S + j) =
                            TDLS_ARROW_LUPP_DYN_A(storage_map::abk(shape, k, r, j));
                }
                ok = multipliers_pass(z, o, work, work_stride, o.Mb + r * S, 0) && ok;
            }
        }
        return ok;
    }
 
    TDLS_HOST_DEVICE static constexpr int
    physical_factor(const Sizes& z, const Offsets& o, T* TDLS_RESTRICT work, const int work_stride,
                    int* TDLS_RESTRICT piv, const int piv_stride) noexcept {
        return PhysicalBlockSolver::physical_factorize_checked(
            z.S, work + TDLS_STRIDED(o.Dl, work_stride), work_stride, piv, piv_stride,
            condensation_pivot_threshold);
    }
 
    TDLS_HOST_DEVICE static constexpr void
    physical_solve(const Sizes& z, const Offsets& o, T* TDLS_RESTRICT work, const int work_stride,
                   const int* TDLS_RESTRICT piv, const int piv_stride, const int w,
                   const bool backward = !multiplier_update) noexcept {
        const int S              = z.S;
        const T* TDLS_RESTRICT F = work + TDLS_STRIDED(o.Dl, work_stride);
        T* TDLS_RESTRICT W       = work + TDLS_STRIDED(o.W, work_stride);
        const std::size_t cs     = TDLS_STRIDED(S, work_stride);
        PhysicalBlockSolver::physical_forward(S, w, F, work_stride, piv, piv_stride, W, work_stride,
                                              cs);
        if (col || backward)
            PhysicalBlockSolver::physical_backward(S, w, F, work_stride, piv, piv_stride, W,
                                                   work_stride, cs);
    }
 
    TDLS_HOST_DEVICE static constexpr bool multipliers_pass(const Sizes& z, const Offsets& o,
                                                            T* TDLS_RESTRICT work,
                                                            const int work_stride, const int m,
                                                            const int first) noexcept {
        const int S = z.S;
        const T u   = condensation_pivot_threshold;
        bool ok     = true;
        for (int l = first; l < S; ++l) {
            TDLS_ARROW_LUPP_DYN_W(m + l) *= TDLS_ARROW_LUPP_DYN_W(o.U + l * S + l);
            ok = ok && u * detail::abs(TDLS_ARROW_LUPP_DYN_W(m + l)) <= T(1);
            for (int j = l + 1; j < S; ++j)
                TDLS_ARROW_LUPP_DYN_W(m + j) -=
                    TDLS_ARROW_LUPP_DYN_W(m + l) * TDLS_ARROW_LUPP_DYN_W(o.U + l * S + j);
        }
        return ok;
    }
 
    template<int NC>
    TDLS_HOST_DEVICE static constexpr void
    schur_sub([[maybe_unused]] const ArrowheadShape& shape, const Sizes& z,
              [[maybe_unused]] const Offsets& o, [[maybe_unused]] T* TDLS_RESTRICT A,
              [[maybe_unused]] const int A_stride, [[maybe_unused]] T* TDLS_RESTRICT work,
              [[maybe_unused]] const int work_stride, [[maybe_unused]] T* TDLS_RESTRICT y,
              [[maybe_unused]] const int rhs_stride, [[maybe_unused]] const std::size_t xcol_stride,
              const int r, const int c, const T v) noexcept {
        [[maybe_unused]] constexpr int ycols = NC;
        if constexpr (in_place) {
            if (c < z.B)
                TDLS_ARROW_LUPP_DYN_A(storage_map::abb(shape, z.R0 + r, c)) -= v;
            else
                TDLS_ARROW_LUPP_DYN_YW(c - z.B, z.R0 + r) -= v;
        } else {
            if (c < z.B)
                TDLS_ARROW_LUPP_DYN_W(o.Sacc + r * z.B + c) -= v;
            else
                TDLS_ARROW_LUPP_DYN_W(o.racc + (NC == 1 ? 0 : (c - z.B) * z.R) + r) -= v;
        }
    }
 
    template<int NC, bool canonical, bool recompute = false, typename... Streamed>
    TDLS_HOST_DEVICE static constexpr void
    condense_pass(const ArrowheadShape& shape, const Sizes& z, const Offsets& o, T* TDLS_RESTRICT A,
                  const int A_stride, const int* TDLS_RESTRICT piv, const int piv_stride,
                  T* TDLS_RESTRICT work, const int work_stride, T* TDLS_RESTRICT y,
                  const int rhs_stride, [[maybe_unused]] const std::size_t xcol_stride, const int k,
                  const int c0, const int w, const bool solved, const Streamed*... so) noexcept {
        constexpr bool streamed              = sizeof...(Streamed) > 0;
        [[maybe_unused]] constexpr int ycols = NC;
        const int B = z.B, S = z.S, R = z.R;
        const T* const Dl = work + TDLS_STRIDED(o.Dl, work_stride);
        T* const W        = work + TDLS_STRIDED(o.W, work_stride);
        if (!solved) {
            if constexpr (streamed) {
                // Columns c0 .. c0 + w of the streamed block move to the
                // front of W: the columns before c0 are solved and stored
                // already, and c0 >= w when c0 > 0, so the two ranges do
                // not overlap.
                if (c0 > 0)
                    for (int i = 0; i < S; ++i)
                        for (int j = 0; j < w; ++j)
                            TDLS_ARROW_LUPP_DYN_W(o.W + j * S + i) =
                                TDLS_ARROW_LUPP_DYN_W(o.W + (c0 + j) * S + i);
            } else {
                for (int i = 0; i < S; ++i)
                    for (int j = 0; j < w; ++j)
                        TDLS_ARROW_LUPP_DYN_W(o.W + j * S + i) =
                            c0 + j < B
                                ? TDLS_ARROW_LUPP_DYN_A(storage_map::akb(shape, k, i, c0 + j))
                                : TDLS_ARROW_LUPP_DYN_YW(c0 + j - B, B + k * S + i);
            }
            if (z.reg) {
                physical_solve(z, o, work, work_stride, piv, piv_stride, w);
            } else if constexpr (multiplier_update) {
                BlockSolver::forward_inplace_multirhs(S, w, Dl, work_stride, piv, piv_stride, W,
                                                      work_stride, TDLS_STRIDED(S, work_stride));
            } else if constexpr (Config.gather_rhs) {
                T* const X = work + TDLS_STRIDED(o.X, work_stride);
                BlockSolver::substitute_multirhs(S, w, Dl, work_stride, piv, piv_stride, W, X,
                                                 work_stride, TDLS_STRIDED(S, work_stride));
                for (int e = 0; e < w * S; ++e)
                    TDLS_ARROW_LUPP_DYN_W(o.W + e) = TDLS_ARROW_LUPP_DYN_W(o.X + e);
            } else {
                BlockSolver::substitute_inplace_multirhs(S, w, Dl, work_stride, piv, piv_stride, W,
                                                         work_stride, TDLS_STRIDED(S, work_stride));
            }
        }
        if constexpr (!canonical && !recompute) {
            for (int i = 0; i < S; ++i)
                for (int j = 0; j < w; ++j) {
                    if (c0 + j < B)
                        TDLS_ARROW_LUPP_DYN_A(storage_map::akb(shape, k, i, c0 + j)) =
                            TDLS_ARROW_LUPP_DYN_W(o.W + j * S + i);
                    else
                        TDLS_ARROW_LUPP_DYN_YW(c0 + j - B, B + k * S + i) =
                            TDLS_ARROW_LUPP_DYN_W(o.W + j * S + i);
                }
        }
        if constexpr (multiplier_update) {
            for (int r = 0; r < R; ++r)
                for (int j = 0; j < w; ++j) {
                    T sum = T(0);
                    for (int l = 0; l < S; ++l)
                        sum += TDLS_ARROW_LUPP_DYN_W(o.Mb + r * S + l) *
                               TDLS_ARROW_LUPP_DYN_W(o.W + j * S + l);
                    schur_sub<NC>(shape, z, o, A, A_stride, work, work_stride, y, rhs_stride,
                                  xcol_stride, r, c0 + j, sum);
                }
        } else if constexpr (diag) {
            for (int i = 0; i < S; ++i)
                for (int j = 0; j < w; ++j)
                    schur_sub<NC>(shape, z, o, A, A_stride, work, work_stride, y, rhs_stride,
                                  xcol_stride, i, c0 + j,
                                  TDLS_ARROW_LUPP_DYN_W(o.acoef + i) *
                                      TDLS_ARROW_LUPP_DYN_W(o.W + j * S + i));
        } else {
            for (int r = 0; r < B; ++r) {
                for (int i = 0; i < S; ++i) {
                    if constexpr (streamed)
                        TDLS_ARROW_LUPP_DYN_W(o.a + i) =
                            TDLS_ARROW_LUPP_DYN_W(streamed_arg(so...)->Abk + r * S + i);
                    else
                        TDLS_ARROW_LUPP_DYN_W(o.a + i) =
                            TDLS_ARROW_LUPP_DYN_A(storage_map::abk(shape, k, r, i));
                }
                for (int j = 0; j < w; ++j) {
                    T sum = T(0);
                    for (int i = 0; i < S; ++i)
                        sum +=
                            TDLS_ARROW_LUPP_DYN_W(o.a + i) * TDLS_ARROW_LUPP_DYN_W(o.W + j * S + i);
                    schur_sub<NC>(shape, z, o, A, A_stride, work, work_stride, y, rhs_stride,
                                  xcol_stride, r, c0 + j, sum);
                }
            }
        }
    }
 
    template<int NC>
    TDLS_HOST_DEVICE static constexpr void
    recover(const ArrowheadShape& shape, const Sizes& z, const Offsets& o, const T* TDLS_RESTRICT A,
            const int A_stride, T* TDLS_RESTRICT work, const int work_stride, T* TDLS_RESTRICT y,
            const int rhs_stride, [[maybe_unused]] const std::size_t xcol_stride,
            [[maybe_unused]] const int nrhs) noexcept {
        const int B = z.B, S = z.S, NB = z.NB, R = z.R;
        const int nc                         = NC > 0 ? NC : nrhs;
        [[maybe_unused]] constexpr int ycols = NC;
        // w and s: column m at m * R and m * S.
        for (int e = 0; e < R * nc; ++e)
            TDLS_ARROW_LUPP_DYN_W(o.w + e) = T(0);
        for (int k = NB - 1; k >= 0; --k) {
            for (int i = 0; i < S; ++i) {
                for (int m = 0; m < nc; ++m)
                    TDLS_ARROW_LUPP_DYN_W(o.s + m * S + i) =
                        TDLS_ARROW_LUPP_DYN_YW(m, B + k * S + i);
                for (int c = 0; c < B; ++c) {
                    const T a = TDLS_ARROW_LUPP_DYN_A(storage_map::akb(shape, k, i, c));
                    for (int m = 0; m < nc; ++m)
                        TDLS_ARROW_LUPP_DYN_W(o.s + m * S + i) -=
                            a * TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + c);
                }
            }
            if (TDLS_ARROW_LUPP_DYN_W(o.refused + k) != T(0)) {
                // x_k = c_k - G_k xb - H_k w, then w += A_bk x_k.
                for (int i = 0; i < S; ++i)
                    for (int a = 0; a < R; ++a) {
                        const T h = a < S ? TDLS_ARROW_LUPP_DYN_A(storage_map::d(shape, k, i, a))
                                          : TDLS_ARROW_LUPP_DYN_W(
                                                o.extra + (k * S + i) * z.extra_cols + a - S);
                        for (int m = 0; m < nc; ++m)
                            TDLS_ARROW_LUPP_DYN_W(o.s + m * S + i) -=
                                h * TDLS_ARROW_LUPP_DYN_W(o.w + m * R + a);
                    }
                if constexpr (diag) {
                    for (int i = 0; i < S; ++i) {
                        const T coef = TDLS_ARROW_LUPP_DYN_A(storage_map::abk_diag(shape, k, i));
                        for (int m = 0; m < nc; ++m)
                            TDLS_ARROW_LUPP_DYN_W(o.w + m * R + i) +=
                                coef * TDLS_ARROW_LUPP_DYN_W(o.s + m * S + i);
                    }
                } else {
                    for (int r = 0; r < B; ++r)
                        for (int i = 0; i < S; ++i) {
                            const T coef = TDLS_ARROW_LUPP_DYN_A(storage_map::abk(shape, k, r, i));
                            for (int m = 0; m < nc; ++m)
                                TDLS_ARROW_LUPP_DYN_W(o.w + m * R + r) +=
                                    coef * TDLS_ARROW_LUPP_DYN_W(o.s + m * S + i);
                        }
                }
            } else if constexpr (multiplier_update) {
                // x_k = U^-1 s, U in the upper triangle of D_k.
                for (int l = S - 1; l >= 0; --l) {
                    for (int j = l + 1; j < S; ++j) {
                        const T ulj = TDLS_ARROW_LUPP_DYN_A(storage_map::d(shape, k, l, j));
                        for (int m = 0; m < nc; ++m)
                            TDLS_ARROW_LUPP_DYN_W(o.s + m * S + l) -=
                                ulj * TDLS_ARROW_LUPP_DYN_W(o.s + m * S + j);
                    }
                    const T ull = TDLS_ARROW_LUPP_DYN_A(storage_map::d(shape, k, l, l));
                    for (int m = 0; m < nc; ++m)
                        TDLS_ARROW_LUPP_DYN_W(o.s + m * S + l) *= ull;
                }
            }
            for (int i = 0; i < S; ++i)
                for (int m = 0; m < nc; ++m)
                    TDLS_ARROW_LUPP_DYN_YW(m, B + k * S + i) =
                        TDLS_ARROW_LUPP_DYN_W(o.s + m * S + i);
        }
    }
 
    template<int NC, bool canonical>
    TDLS_HOST_DEVICE TDLS_NO_INLINE static constexpr bool
    second_pass(const ArrowheadShape& shape, const Sizes& z, const Offsets& o, T* TDLS_RESTRICT A,
                const int A_stride, T* TDLS_RESTRICT work, const int work_stride,
                T* TDLS_RESTRICT y, const int rhs_stride,
                [[maybe_unused]] const std::size_t xcol_stride,
                [[maybe_unused]] const int nrhs) noexcept {
        const int B = z.B, S = z.S, NB = z.NB, R = z.R, R0 = z.R0;
        const int front_rows = z.front_rows, front_cols = z.front_cols;
        const int nc                         = NC > 0 ? NC : nrhs;
        const int rhs0                       = S + B + R;
        [[maybe_unused]] constexpr int ycols = NC;
        const T u                            = condensation_pivot_threshold;
        // The singularity floor of the block solver: a smaller pivot, zero
        // or subnormal, is singular there, and its reciprocal could overflow.
        constexpr T singular_floor = BlockSolver::singular_floor;
        for (int r = 0; r < B; ++r) {
            for (int c = 0; c < B; ++c)
                TDLS_ARROW_LUPP_DYN_F(S + r, S + c) = TDLS_ARROW_LUPP_DYN_W(o.Sbb + r * B + c);
            for (int a = 0; a < R; ++a)
                TDLS_ARROW_LUPP_DYN_F(S + r, S + B + a) = r == R0 + a ? T(1) : T(0);
            for (int m = 0; m < nc; ++m)
                TDLS_ARROW_LUPP_DYN_F(S + r, rhs0 + m) = TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + r);
        }
        for (int k = 0; k < NB; ++k) {
            if (TDLS_ARROW_LUPP_DYN_W(o.refused + k) == T(0)) continue;
            for (int i = 0; i < S; ++i) {
                for (int j = 0; j < S; ++j)
                    TDLS_ARROW_LUPP_DYN_F(i, j) =
                        TDLS_ARROW_LUPP_DYN_A(storage_map::d(shape, k, i, j));
                for (int c = 0; c < B; ++c)
                    TDLS_ARROW_LUPP_DYN_F(i, S + c) =
                        TDLS_ARROW_LUPP_DYN_A(storage_map::akb(shape, k, i, c));
                for (int a = 0; a < R; ++a)
                    TDLS_ARROW_LUPP_DYN_F(i, S + B + a) = T(0);
                for (int m = 0; m < nc; ++m) {
                    if constexpr (canonical)
                        TDLS_ARROW_LUPP_DYN_F(i, rhs0 + m) = T(0);
                    else
                        TDLS_ARROW_LUPP_DYN_F(i, rhs0 + m) =
                            TDLS_ARROW_LUPP_DYN_YW(m, B + k * S + i);
                }
            }
            for (int r = 0; r < B; ++r) {
                for (int j = 0; j < S; ++j) {
                    T v = T(0);
                    if constexpr (diag) {
                        v = TDLS_ARROW_LUPP_DYN_F(S + r, S + B + j) *
                            TDLS_ARROW_LUPP_DYN_A(storage_map::abk_diag(shape, k, j));
                    } else {
                        for (int a = 0; a < B; ++a)
                            v += TDLS_ARROW_LUPP_DYN_F(S + r, S + B + a) *
                                 TDLS_ARROW_LUPP_DYN_A(storage_map::abk(shape, k, a, j));
                    }
                    TDLS_ARROW_LUPP_DYN_F(S + r, j) = v;
                }
            }
            for (int e = 0; e < front_rows; ++e)
                TDLS_ARROW_LUPP_DYN_W(o.block_row + e) = e < S ? T(1) : T(0);
            for (int j = 0; j < S; ++j) {
                int pmax = j, pblk = -1;
                T vmax = T(0), vblk = T(0);
                for (int e = j; e < front_rows; ++e) {
                    const T v = detail::abs(TDLS_ARROW_LUPP_DYN_F(e, j));
                    if (v > vmax) {
                        vmax = v;
                        pmax = e;
                    }
                    if (TDLS_ARROW_LUPP_DYN_W(o.block_row + e) != T(0) && v > vblk) {
                        vblk = v;
                        pblk = e;
                    }
                }
                if (!(vmax >= singular_floor)) return false;
                const int p = pblk >= 0 && vblk >= u * vmax && vblk >= singular_floor ? pblk : pmax;
                if (p != j) {
                    for (int c = j; c < front_cols; ++c) {
                        const T t                   = TDLS_ARROW_LUPP_DYN_F(p, c);
                        TDLS_ARROW_LUPP_DYN_F(p, c) = TDLS_ARROW_LUPP_DYN_F(j, c);
                        TDLS_ARROW_LUPP_DYN_F(j, c) = t;
                    }
                    const T t                              = TDLS_ARROW_LUPP_DYN_W(o.block_row + p);
                    TDLS_ARROW_LUPP_DYN_W(o.block_row + p) = TDLS_ARROW_LUPP_DYN_W(o.block_row + j);
                    TDLS_ARROW_LUPP_DYN_W(o.block_row + j) = t;
                }
                const T inv                 = T(1) / TDLS_ARROW_LUPP_DYN_F(j, j);
                TDLS_ARROW_LUPP_DYN_F(j, j) = inv;
                for (int e = j + 1; e < front_rows; ++e) {
                    const T l = TDLS_ARROW_LUPP_DYN_F(e, j) * inv;
                    if (l == T(0)) continue;
                    for (int c = j + 1; c < front_cols; ++c)
                        TDLS_ARROW_LUPP_DYN_F(e, c) -= l * TDLS_ARROW_LUPP_DYN_F(j, c);
                }
            }
            if constexpr (!canonical) {
                for (int c = S; c < front_cols; ++c) {
                    for (int i = S - 1; i >= 0; --i) {
                        T t = TDLS_ARROW_LUPP_DYN_F(i, c);
                        for (int m = i + 1; m < S; ++m)
                            t -= TDLS_ARROW_LUPP_DYN_F(i, m) * TDLS_ARROW_LUPP_DYN_F(m, c);
                        TDLS_ARROW_LUPP_DYN_F(i, c) = t * TDLS_ARROW_LUPP_DYN_F(i, i);
                    }
                }
                for (int i = 0; i < S; ++i) {
                    for (int c = 0; c < B; ++c)
                        TDLS_ARROW_LUPP_DYN_A(storage_map::akb(shape, k, i, c)) =
                            TDLS_ARROW_LUPP_DYN_F(i, S + c);
                    for (int a = 0; a < R; ++a) {
                        if (a < S)
                            TDLS_ARROW_LUPP_DYN_A(storage_map::d(shape, k, i, a)) =
                                TDLS_ARROW_LUPP_DYN_F(i, S + B + a);
                        else
                            TDLS_ARROW_LUPP_DYN_W(o.extra + (k * S + i) * z.extra_cols + a - S) =
                                TDLS_ARROW_LUPP_DYN_F(i, S + B + a);
                    }
                    for (int m = 0; m < nc; ++m)
                        TDLS_ARROW_LUPP_DYN_YW(m, B + k * S + i) =
                            TDLS_ARROW_LUPP_DYN_F(i, rhs0 + m);
                }
            }
        }
        for (int r = 0; r < B; ++r) {
            for (int c = 0; c < B; ++c)
                TDLS_ARROW_LUPP_DYN_W(o.Sbb + r * B + c) = TDLS_ARROW_LUPP_DYN_F(S + r, S + c);
            for (int m = 0; m < nc; ++m)
                TDLS_ARROW_LUPP_DYN_W(o.rb + m * B + r) = TDLS_ARROW_LUPP_DYN_F(S + r, rhs0 + m);
        }
        return true;
    }
};
 
 
 
#undef TDLS_ARROW_LUPP_DYN_A
#undef TDLS_ARROW_LUPP_DYN_Y
#undef TDLS_ARROW_LUPP_DYN_YW
#undef TDLS_ARROW_LUPP_DYN_W
#undef TDLS_ARROW_LUPP_DYN_F
#undef TDLS_ARROW_LUPP_DYN_LU
#undef TDLS_ARROW_LUPP_DYN_RC
 
 
 
} // namespace tdls
 
 
 
#endif // TDLS_SOLVERS_ARROWHEAD_LUPP_SOLVER_DYNAMIC_HPP