gemm.hpp

Go to the documentation of this file.

 
#ifndef NUMERICAL_GEMM_HPP
#define NUMERICAL_GEMM_HPP
 
#include <algorithm>
#include <array>
#include <bit>
#include <cmath>
#include <atomic>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <memory>
#include <type_traits>
#include <utility>
#include <vector>
 
#if defined(__GNUC__) && defined(__aarch64__)
#  include <arm_neon.h>
#  define NUMERICAL_NEON_COMPILED // the arm64 fused tile is compiled in (tested with #ifdef only)
#endif
 
#include "parallel.hpp"
#include "amx.hpp"
#include "buffer.hpp"
#include "simd.hpp"
#include "x86.hpp"
 
// A function kept out of its callers: inlined, it spent GCC's growth budget for the whole translation unit, and the
// products, the Cholesky factorisation and the LU solve lost their own inlining (one arm64 thread, GCC 14: a product of
// order 8 0.29 → 0.33 µs, a Cholesky factorisation of order 17 0.77 → 0.87 µs, a solve of order 512 0.47 → 0.54 ms).
#if defined(_MSC_VER) && !defined(__clang__)
#  define NUMERICAL_OUT_OF_LINE __declspec(noinline)
#else
#  define NUMERICAL_OUT_OF_LINE [[gnu::noinline]] // NOLINT(cppcoreguidelines-macro-usage): an attribute, which no function can stand for
#endif
 
#if defined(__clang__)
// A block's depth loop kept scalar: vectorised, Clang sums its terms in order but splits each multiply-add (which
// -ffp-contract=on lets it fuse, not makes it) into a product and a sum, the last terms fused: one element rounded
// two ways, where the packed kernels fuse every term.
#  define NUMERICAL_SCALAR_DEPTH _Pragma("clang loop vectorize(disable) interleave(disable)")
#else
#  define NUMERICAL_SCALAR_DEPTH
#endif
 
#if defined(__GNUC__)
// A block's loops unrolled whole, so that each of its sums stays in a register: GCC at -O2 does not unroll them,
// and kept the sums in memory (a product of order 17 on one arm64 core: 2.0 µs, against 0.31 under Clang).
#  define NUMERICAL_UNROLL_BLOCK _Pragma("GCC unroll 32")
#else
#  define NUMERICAL_UNROLL_BLOCK
#endif
 
namespace numerical {
 
  enum class gemm_path : std::uint8_t
  {
    portable = 0, 
    neon     = 1, 
    avx2     = 2, 
    avx512   = 3, 
    amx      = 4, 
  };
} // namespace numerical
 
namespace numerical::detail {
 
#ifdef NUMERICAL_NEON_COMPILED
  inline constexpr bool neon_compiled {true};
#else
  inline constexpr bool neon_compiled {false};
#endif
 
  [[nodiscard]] inline bool gemm_path_usable(gemm_path path) noexcept
  {
    switch (path) {
#ifdef NUMERICAL_X86_COMPILED
      case gemm_path::avx2:   return cpu_features().avx2_fma;
      case gemm_path::avx512: return cpu_features().avx512f;
#endif
      case gemm_path::neon:   return neon_compiled;
      case gemm_path::amx:    return amx_usable();
      default:                return path == gemm_path::portable;
    }
  }
 
  [[nodiscard]] inline gemm_path best_gemm_path(gemm_path most) noexcept
  {
    for (auto p {static_cast<int>(most)}; p > 0; --p) {
      if (gemm_path_usable(static_cast<gemm_path>(p))) {
        return static_cast<gemm_path>(p);
      }
    }
    return gemm_path::portable;
  }
 
  [[nodiscard]] inline std::atomic<gemm_path>& gemm_path_cap() noexcept
  {
    static std::atomic<gemm_path> cap {gemm_path::amx};
    return cap;
  }
 
  [[nodiscard]] inline std::atomic<std::uint8_t>& gemm_path_found() noexcept
  {
    static std::atomic<std::uint8_t> found {0xFF};
    return found;
  }
 
  [[nodiscard]] inline gemm_path gemm_path_current() noexcept
  {
    const unsigned cap   {static_cast<unsigned>(gemm_path_cap().load(std::memory_order_relaxed))};
    const unsigned found {gemm_path_found().load(std::memory_order_relaxed)};
    if (found >> 3U == cap) {
      return static_cast<gemm_path>(found & 7U);
    }
    const gemm_path path {best_gemm_path(static_cast<gemm_path>(cap))};
    gemm_path_found().store(static_cast<std::uint8_t>((cap << 3U) | static_cast<unsigned>(path)), std::memory_order_relaxed);
    return path;
  }
 
  inline void subtract_product(double*       y,
                               const double& f,
                               const double* x,
                               std::size_t   count)
  {
    const double factor  {f}; // read once: it may be an element of y's array
#ifdef NUMERICAL_X86_COMPILED
    const gemm_path path {gemm_path_current()};
    if (path == gemm_path::avx512) {
      avx512_scaled_add(y, x, -factor, count);
      return;
    }
    if (path == gemm_path::avx2) {
      avx2_scaled_add(y, x, -factor, count);
      return;
    }
#endif
    subtract_product<double>(y, factor, x, count);
  }
 
  inline void subtract_product(double*       y,
                               const double* x,
                               const double& f,
                               std::size_t   count)
  {
    subtract_product(y, f, x, count);
  }
 
  inline void subtract_products(double*       y,
                                const double* f,
                                const double* x,
                                std::size_t   ldx,
                                std::size_t   terms,
                                std::size_t   count)
  {
#ifdef NUMERICAL_X86_COMPILED
    const gemm_path path {gemm_path_current()};
    if (path == gemm_path::avx512) {
      avx512_subtract_products(y, f, x, ldx, terms, count);
      return;
    }
    if (path == gemm_path::avx2) {
      avx2_subtract_products(y, f, x, ldx, terms, count);
      return;
    }
#endif
    subtract_products<double>(y, f, x, ldx, terms, count);
  }
 
  inline void subtract_rows(double*       y,
                            std::size_t   ldy,
                            const double* f,
                            std::size_t   ldf,
                            const double* x,
                            std::size_t   ldx,
                            std::size_t   terms,
                            std::size_t   count)
  {
#ifdef NUMERICAL_X86_COMPILED
    const gemm_path path {gemm_path_current()};
    if (path == gemm_path::avx512) {
      avx512_subtract_rows(y, ldy, f, ldf, x, ldx, terms, count);
      return;
    }
    if (path == gemm_path::avx2) {
      avx2_subtract_rows(y, ldy, f, ldf, x, ldx, terms, count);
      return;
    }
#endif
    subtract_rows<double>(y, ldy, f, ldf, x, ldx, terms, count);
  }
 
  inline void subtract_products(std::complex<double>*       y,
                                const std::complex<double>* f,
                                const std::complex<double>* x,
                                std::size_t                 ldx,
                                std::size_t                 terms,
                                std::size_t                 count)
  {
#ifdef NUMERICAL_X86_COMPILED
    // A complex number's parts are an array of two (the standard's layout guarantee): addressed as doubles.
    const gemm_path path {gemm_path_current()};
    if (path == gemm_path::avx512 || path == gemm_path::avx2) {
      auto* const       into   {reinterpret_cast<double*>(y)};       // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee above
      const auto* const scales {reinterpret_cast<const double*>(f)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee above
      const auto* const rows   {reinterpret_cast<const double*>(x)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee above
      if (path == gemm_path::avx512) {
        avx512_subtract_complex_products(into, scales, rows, ldx, terms, count);
      }
      else {
        avx2_subtract_complex_products(into, scales, rows, ldx, terms, count);
      }
      return;
    }
#endif
    subtract_products<std::complex<double>>(y, f, x, ldx, terms, count);
  }
 
  inline void subtract_product(std::complex<double>*       y,
                               const std::complex<double>& f,
                               const std::complex<double>* x,
                               std::size_t                 count)
  {
    const std::complex<double> factor {f}; // read once: it may be an element of y's array
#ifdef NUMERICAL_X86_COMPILED
    const gemm_path path              {gemm_path_current()};
    if (path == gemm_path::avx512 || path == gemm_path::avx2) {
      subtract_products(y, &factor, x, 0, 1, count);
      return;
    }
#endif
    subtract_product<std::complex<double>>(y, factor, x, count);
  }
 
  inline void subtract_product(std::complex<double>*       y,
                               const std::complex<double>* x,
                               const std::complex<double>& f,
                               std::size_t                 count)
  {
    const std::complex<double> factor {f}; // read once: it may be an element of y's array
#ifdef NUMERICAL_X86_COMPILED
    const gemm_path path              {gemm_path_current()};
    if (path == gemm_path::avx512 || path == gemm_path::avx2) {
      subtract_products(y, &factor, x, 0, 1, count);
      return;
    }
#endif
    subtract_product<std::complex<double>>(y, x, factor, count);
  }
 
  inline void add_product(double*       y,
                          const double& f,
                          const double* x,
                          std::size_t   count)
  {
    const double factor  {f}; // read once: it may be an element of y's array
#ifdef NUMERICAL_X86_COMPILED
    const gemm_path path {gemm_path_current()};
    if (path == gemm_path::avx512) {
      avx512_scaled_add(y, x, factor, count);
      return;
    }
    if (path == gemm_path::avx2) {
      avx2_scaled_add(y, x, factor, count);
      return;
    }
#endif
    add_product<double>(y, factor, x, count);
  }
 
  inline void add_product(std::complex<float>*       y,
                          const std::complex<float>& f,
                          const std::complex<float>* x,
                          std::size_t                count)
  {
    const std::complex<float> factor {f}; // read once: it may be an element of y's array
#ifdef NUMERICAL_X86_COMPILED
    const gemm_path path             {gemm_path_current()};
    if (path == gemm_path::avx512 || path == gemm_path::avx2) {
      for (std::size_t i {avx2_scaled_add_c32(y, x, factor, count)}; i < count; ++i) {
        complex_scaled_add(y[i], factor, x[i]);
      }
      return;
    }
#endif
    add_product<std::complex<float>>(y, factor, x, count);
  }
 
  inline void add_product(std::complex<double>*       y,
                          const std::complex<double>& f,
                          const std::complex<double>* x,
                          std::size_t                 count)
  {
    const std::complex<double> factor {f}; // read once: it may be an element of y's array
#ifdef NUMERICAL_X86_COMPILED
    const gemm_path path              {gemm_path_current()};
    if (path == gemm_path::avx512 || path == gemm_path::avx2) {
      for (std::size_t i {avx2_scaled_add_c64(y, x, factor, count)}; i < count; ++i) {
        complex_scaled_add(y[i], factor, x[i]);
      }
      return;
    }
#endif
    add_product<std::complex<double>>(y, factor, x, count);
  }
 
  [[nodiscard]] inline double dot_product(const double* x,
                                          const double* y,
                                          std::size_t   count)
  {
#ifdef NUMERICAL_X86_COMPILED
    const gemm_path path {gemm_path_current()};
    if (path == gemm_path::avx512) {
      return avx512_dot(x, y, count);
    }
    if (path == gemm_path::avx2) {
      return avx2_dot(x, y, count);
    }
#endif
    return dot_product<double>(x, y, count);
  }
 
  template <split_complex C>
  [[nodiscard]] C complex_dot_product(const C*    x,
                                      const C*    y,
                                      std::size_t count)
  {
#ifdef NUMERICAL_X86_COMPILED
    if (gemm_path_current() != gemm_path::portable) {
      return avx2_complex_dot(x, y, count);
    }
#endif
    return dot_product<C>(x, y, count); // the native vectors' (the overloads below come here)
  }
 
  [[nodiscard]] inline std::complex<double> dot_product(const std::complex<double>* x,
                                                        const std::complex<double>* y,
                                                        std::size_t                 count)
  {
    return complex_dot_product(x, y, count);
  }
 
  [[nodiscard]] inline std::complex<float> dot_product(const std::complex<float>* x,
                                                       const std::complex<float>* y,
                                                       std::size_t                count)
  {
    return complex_dot_product(x, y, count);
  }
 
  template <split_complex C>
  void complex_row_dots(const C*    m,
                        std::size_t ld,
                        std::size_t rows,
                        const C*    x,
                        std::size_t count,
                        C*          out)
  {
    for (std::size_t r {0}; r < rows; ++r) {
      out[r] = complex_dot_product(m + (r * ld), x, count);
    }
  }
 
  inline void row_dots(const std::complex<double>* m,
                       std::size_t                 ld,
                       std::size_t                 rows,
                       const std::complex<double>* x,
                       std::size_t                 count,
                       std::complex<double>*       out)
  {
    complex_row_dots(m, ld, rows, x, count, out);
  }
 
  inline void row_dots(const std::complex<float>* m,
                       std::size_t                ld,
                       std::size_t                rows,
                       const std::complex<float>* x,
                       std::size_t                count,
                       std::complex<float>*       out)
  {
    complex_row_dots(m, ld, rows, x, count, out);
  }
 
  inline void row_dots(const double* m,
                       std::size_t   ld,
                       std::size_t   rows,
                       const double* x,
                       std::size_t   count,
                       double*       out)
  {
#ifdef NUMERICAL_X86_COMPILED
    if (gemm_path_current() != gemm_path::portable) {
      for (std::size_t r {0}; r < rows; ++r) {
        out[r] = dot_product(m + (r * ld), x, count);
      }
      return;
    }
#endif
    row_dots<double>(m, ld, rows, x, count, out);
  }
 
  inline constexpr std::size_t gemm_mr {4};
 
  template <typename T>
  inline constexpr std::size_t gemm_nr {sizeof(T) >= 16 ? std::size_t {6} : std::size_t {96} / sizeof(T)};
 
  inline constexpr std::size_t gemm_kc {256};
 
  inline constexpr std::size_t gemm_mc {64};
 
  inline constexpr std::size_t gemm_slivers_per_task {4};
 
  inline constexpr std::size_t tasks_per_thread {4};
 
  inline constexpr std::size_t gemm_kernel_speedup {16};
 
  template <typename T>
  inline constexpr std::size_t gemm_nc {gemm_nr<T> * 128};
 
#ifdef NUMERICAL_X86_COMPILED
  inline constexpr std::size_t gemm_packed_threshold {std::size_t {1} << 11U};
#else
  inline constexpr std::size_t gemm_packed_threshold {std::size_t {1} << 13U};
#endif
 
  inline constexpr std::size_t gemm_min_depth {8};
 
  template <typename T>
  using gemm_tile = std::array<T, gemm_mr * gemm_nr<T>>;
 
  template <typename T>
  void gemm_kernel(std::size_t   kc,
                   const T*      ap,
                   const T*      bp,
                   std::size_t   ldb,
                   gemm_tile<T>& tile)
  {
    constexpr std::size_t nr {gemm_nr<T>};
    for (std::size_t p {0}; p < kc; ++p) {
      for (std::size_t i {0}; i < gemm_mr; ++i) {
        const T ai {ap[(p * gemm_mr) + i]};
        for (std::size_t j {0}; j < nr; ++j) {
          tile[(i * nr) + j] += ai * bp[(p * ldb) + j];
        }
      }
    }
  }
 
#if defined(__GNUC__)
  template <typename T>
  requires std::is_same_v<T, float> || std::is_same_v<T, double>
  void gemm_kernel(std::size_t   kc,
                   const T*      ap,
                   const T*      bp,
                   std::size_t   ldb,
                   gemm_tile<T>& tile)
  {
    using lane                 = typename native_vector<T>::type;
    constexpr std::size_t             nr    {gemm_nr<T>};
    constexpr std::size_t             lanes {nr * sizeof(T) / native_vector_bytes};
    // The tile moves in and out one vector at a time: a single copy of the whole array would leave it in
    // memory under GCC, which then keeps no accumulator in a register (7.4 GFLOP/s instead of 11.8 on one
    // x86-64 core).
    std::array<lane, gemm_mr * lanes> acc {};
#  pragma GCC unroll 32
    for (std::size_t q {0}; q < acc.size(); ++q) {
      lane v;
      std::memcpy(&v, tile.data() + (q * (native_vector_bytes / sizeof(T))), sizeof(v));
      acc[q] = v;
    }
    // Unrolled, so that every accumulator has a fixed index and stays in a register: GCC at -O2 does not
    // unroll on its own, and then reads and writes the tile through memory.
    for (std::size_t p {0}; p < kc; ++p) {
      std::array<lane, lanes> b {};
      std::memcpy(b.data(), bp + (p * ldb), sizeof(b));
#  pragma GCC unroll 8
      for (std::size_t i {0}; i < gemm_mr; ++i) {
        const T ai {ap[(p * gemm_mr) + i]};
#  pragma GCC unroll 8
        for (std::size_t j {0}; j < lanes; ++j) {
          acc[(i * lanes) + j] += ai * b[j];
        }
      }
    }
#  pragma GCC unroll 32
    for (std::size_t q {0}; q < acc.size(); ++q) {
      const lane v {acc[q]};
      std::memcpy(tile.data() + (q * (native_vector_bytes / sizeof(T))), &v, sizeof(v));
    }
  }
 
#endif
 
#ifdef NUMERICAL_NEON_COMPILED
  inline void neon_kernel(std::size_t        kc,
                          const double*      ap,
                          const double*      bp,
                          std::size_t        ldb,
                          gemm_tile<double>& tile) noexcept
  {
    static_assert(gemm_mr == 4 && gemm_nr<double> == 12, "two rows of A a vector, twelve columns six vectors");
    constexpr std::size_t              lanes {6};
    std::array<float64x2_t, 4 * lanes> acc   {};
    NUMERICAL_UNROLL_BLOCK
    for (std::size_t q {0}; q < acc.size(); ++q) {
      acc[q] = vld1q_f64(tile.data() + (2 * q));
    }
    for (std::size_t p {0}; p < kc; ++p) {
      const float64x2_t top    {vld1q_f64(ap + (p * gemm_mr))};
      const float64x2_t bottom {vld1q_f64(ap + (p * gemm_mr) + 2)};
      NUMERICAL_UNROLL_BLOCK
      for (std::size_t j {0}; j < lanes; ++j) {
        const float64x2_t b {vld1q_f64(bp + (p * ldb) + (2 * j))};
        acc[j]               = vfmaq_laneq_f64(acc[j], b, top, 0);
        acc[lanes + j]       = vfmaq_laneq_f64(acc[lanes + j], b, top, 1);
        acc[(2 * lanes) + j] = vfmaq_laneq_f64(acc[(2 * lanes) + j], b, bottom, 0);
        acc[(3 * lanes) + j] = vfmaq_laneq_f64(acc[(3 * lanes) + j], b, bottom, 1);
      }
    }
    NUMERICAL_UNROLL_BLOCK
    for (std::size_t q {0}; q < acc.size(); ++q) {
      vst1q_f64(tile.data() + (2 * q), acc[q]);
    }
  }
 
  inline void neon_kernel(std::size_t       kc,
                          const float*      ap,
                          const float*      bp,
                          std::size_t       ldb,
                          gemm_tile<float>& tile) noexcept
  {
    static_assert(gemm_mr == 4 && gemm_nr<float> == 24, "four rows of A a vector, 24 columns six vectors");
    constexpr std::size_t              lanes {6};
    std::array<float32x4_t, 4 * lanes> acc   {};
    NUMERICAL_UNROLL_BLOCK
    for (std::size_t q {0}; q < acc.size(); ++q) {
      acc[q] = vld1q_f32(tile.data() + (4 * q));
    }
    for (std::size_t p {0}; p < kc; ++p) {
      const float32x4_t column {vld1q_f32(ap + (p * gemm_mr))};
      NUMERICAL_UNROLL_BLOCK
      for (std::size_t j {0}; j < lanes; ++j) {
        const float32x4_t b {vld1q_f32(bp + (p * ldb) + (4 * j))};
        acc[j]               = vfmaq_laneq_f32(acc[j], b, column, 0);
        acc[lanes + j]       = vfmaq_laneq_f32(acc[lanes + j], b, column, 1);
        acc[(2 * lanes) + j] = vfmaq_laneq_f32(acc[(2 * lanes) + j], b, column, 2);
        acc[(3 * lanes) + j] = vfmaq_laneq_f32(acc[(3 * lanes) + j], b, column, 3);
      }
    }
    NUMERICAL_UNROLL_BLOCK
    for (std::size_t q {0}; q < acc.size(); ++q) {
      vst1q_f32(tile.data() + (4 * q), acc[q]);
    }
  }
 
#endif
 
  template <std::size_t Mr,
            typename T>
  void gemm_pack_a(const T*    a,
                   std::size_t lda,
                   std::size_t rows,
                   std::size_t kc,
                   T*          out)
  {
    for (std::size_t ir {0}; ir < rows; ir += Mr) {
      const std::size_t mr     {std::min(Mr, rows - ir)};
      T*                sliver {out + (ir * kc)};
      if (mr == Mr) { // a whole sliver, its loop free of the padding's test (M1 Pro, order 64 for the AMX: 1.21 → 0.96 µs)
        for (std::size_t p {0}; p < kc; ++p) {
          for (std::size_t i {0}; i < Mr; ++i) {
            sliver[(p * Mr) + i] = a[((ir + i) * lda) + p];
          }
        }
        continue;
      }
      for (std::size_t p {0}; p < kc; ++p) {
        for (std::size_t i {0}; i < Mr; ++i) {
          sliver[(p * Mr) + i] = i < mr ? a[((ir + i) * lda) + p] : T {};
        }
      }
    }
  }
 
  template <std::size_t Nr,
            typename T>
  void gemm_pack_b(const T*    b,
                   std::size_t ldb,
                   std::size_t kc,
                   std::size_t cols,
                   T*          out)
  {
    for (std::size_t jr {0}; jr < cols; jr += Nr) {
      const std::size_t width  {std::min(Nr, cols - jr)};
      T*                sliver {out + (jr * kc)};
      for (std::size_t p {0}; p < kc; ++p) {
        for (std::size_t j {0}; j < Nr; ++j) {
          sliver[(p * Nr) + j] = j < width ? b[(p * ldb) + jr + j] : T {};
        }
      }
    }
  }
 
  template <std::size_t Nr,
            typename T>
  void gemm_pack_bt(const T*    bt,
                    std::size_t ldbt,
                    std::size_t kc,
                    std::size_t cols,
                    T*          out)
  {
    for (std::size_t jr {0}; jr < cols; jr += Nr) {
      const std::size_t width  {std::min(Nr, cols - jr)};
      T*                sliver {out + (jr * kc)};
      for (std::size_t j {0}; j < Nr; ++j) {
        const T* const row {bt + ((jr + j) * ldbt)};
        for (std::size_t p {0}; p < kc; ++p) {
          sliver[(p * Nr) + j] = j < width ? row[p] : T {};
        }
      }
    }
  }
 
  template <typename T>
  [[nodiscard]] const T* rows_in_place(const T*     x,
                                       std::size_t  ld,
                                       std::size_t& ld_out) noexcept
  {
    ld_out = ld;
    return x;
  }
 
  template <typename T>
  [[nodiscard]] const T* rows_from_transpose(const T*     x,
                                             std::size_t  ld,
                                             std::size_t  rows,
                                             std::size_t  cols,
                                             buffer<T>&   scratch,
                                             std::size_t& ld_out)
  {
    scratch.resize(rows * cols);
    for (std::size_t j {0}; j < cols; ++j) {
      for (std::size_t i {0}; i < rows; ++i) {
        scratch[(i * cols) + j] = x[(j * ld) + i];
      }
    }
    ld_out = cols;
    return scratch.data();
  }
 
  struct b_rows
  {
    template <std::size_t Nr,
              typename T>
    static void pack(const T*    b,
                     std::size_t ldb,
                     std::size_t p0,
                     std::size_t j0,
                     std::size_t kc,
                     std::size_t cols,
                     T*          out)
    {
      gemm_pack_b<Nr>(b + (p0 * ldb) + j0, ldb, kc, cols, out);
    }
 
    template <typename Kernel>
    [[nodiscard]] static constexpr bool in_place(std::size_t ldb) noexcept
    {
      return Kernel::reads_b_in_place(ldb);
    }
 
    template <typename T>
    [[nodiscard]] static const T* rows_of(const T*     x,
                                          std::size_t  ld,
                                          std::size_t  /*rows*/,
                                          std::size_t  /*cols*/,
                                          buffer<T>&   /*scratch*/,
                                          std::size_t& ld_out)
    {
      return rows_in_place(x, ld, ld_out);
    }
 
    template <typename T>
    [[nodiscard]] static const T& at(const T*    b,
                                     std::size_t ldb,
                                     std::size_t p,
                                     std::size_t j) noexcept
    {
      return b[(p * ldb) + j];
    }
  };
 
  struct b_transposed
  {
    template <std::size_t Nr,
              typename T>
    static void pack(const T*    bt,
                     std::size_t ldbt,
                     std::size_t p0,
                     std::size_t j0,
                     std::size_t kc,
                     std::size_t cols,
                     T*          out)
    {
      gemm_pack_bt<Nr>(bt + (j0 * ldbt) + p0, ldbt, kc, cols, out);
    }
 
    template <typename Kernel>
    [[nodiscard]] static constexpr bool in_place(std::size_t /*ldbt*/) noexcept
    {
      return false;
    }
 
    template <typename T>
    [[nodiscard]] static const T* rows_of(const T*     x,
                                          std::size_t  ld,
                                          std::size_t  rows,
                                          std::size_t  cols,
                                          buffer<T>&   scratch,
                                          std::size_t& ld_out)
    {
      return rows_from_transpose(x, ld, rows, cols, scratch, ld_out);
    }
 
    template <typename T>
    [[nodiscard]] static const T& at(const T*    bt,
                                     std::size_t ldbt,
                                     std::size_t p,
                                     std::size_t j) noexcept
    {
      return bt[(j * ldbt) + p];
    }
  };
 
  template <std::size_t Mr,
            typename T>
  void gemm_pack_at(const T*    at,
                    std::size_t ldat,
                    std::size_t rows,
                    std::size_t kc,
                    T*          out)
  {
    for (std::size_t ir {0}; ir < rows; ir += Mr) {
      const std::size_t mr     {std::min(Mr, rows - ir)};
      T*                sliver {out + (ir * kc)};
      for (std::size_t p {0}; p < kc; ++p) {
        const T* const column {at + (p * ldat) + ir};
        for (std::size_t i {0}; i < Mr; ++i) {
          sliver[(p * Mr) + i] = i < mr ? column[i] : T {};
        }
      }
    }
  }
 
  struct a_rows
  {
    template <typename Kernel,
              typename T>
    static void pack(const T*    a,
                     std::size_t lda,
                     std::size_t i0,
                     std::size_t p0,
                     std::size_t rows,
                     std::size_t kc,
                     T*          out)
    {
      Kernel::pack_a(a + (i0 * lda) + p0, lda, rows, kc, out);
    }
 
    template <typename T>
    [[nodiscard]] static const T* rows_of(const T*     a,
                                          std::size_t  lda,
                                          std::size_t  /*rows*/,
                                          std::size_t  /*cols*/,
                                          buffer<T>&   /*scratch*/,
                                          std::size_t& ld_out)
    {
      return rows_in_place(a, lda, ld_out);
    }
 
    template <typename T>
    [[nodiscard]] static const T& at(const T*    a,
                                     std::size_t lda,
                                     std::size_t i,
                                     std::size_t p) noexcept
    {
      return a[(i * lda) + p];
    }
  };
 
  struct a_transposed
  {
    template <typename Kernel,
              typename T>
    static void pack(const T*    at,
                     std::size_t ldat,
                     std::size_t i0,
                     std::size_t p0,
                     std::size_t rows,
                     std::size_t kc,
                     T*          out)
    {
      gemm_pack_at<Kernel::mr>(at + (p0 * ldat) + i0, ldat, rows, kc, out);
    }
 
    template <typename T>
    [[nodiscard]] static const T* rows_of(const T*     at,
                                          std::size_t  ldat,
                                          std::size_t  rows,
                                          std::size_t  cols,
                                          buffer<T>&   scratch,
                                          std::size_t& ld_out)
    {
      return rows_from_transpose(at, ldat, rows, cols, scratch, ld_out);
    }
 
    template <typename T>
    [[nodiscard]] static const T& at(const T*    at_,
                                     std::size_t ldat,
                                     std::size_t i,
                                     std::size_t p) noexcept
    {
      return at_[(p * ldat) + i];
    }
  };
 
  template <typename T>
  struct vector_tile_kernel
  {
    static constexpr std::size_t mr           {gemm_mr};             
    static constexpr std::size_t nr           {gemm_nr<T>};          
    static constexpr std::size_t depth        {gemm_kc};             
    static constexpr std::size_t block_rows   {gemm_mc};             
    static constexpr std::size_t b_in_place_rows {mr};
    static constexpr std::size_t speedup         {gemm_kernel_speedup}; 
    static constexpr std::size_t serial_below    {1};                   
 
    struct session
    {};
 
    static constexpr bool reads_b_in_place(std::size_t ldb) noexcept
    {
      return ldb <= nr;
    }
 
    static void pack_a(const T*    a,
                       std::size_t lda,
                       std::size_t rows,
                       std::size_t kc,
                       T*          out)
    {
      gemm_pack_a<mr>(a, lda, rows, kc, out);
    }
 
    static void tile(std::size_t kc,
                     const T*    ap,
                     const T*    bp,
                     std::size_t ldb,
                     T*          c,
                     std::size_t ldc,
                     std::size_t rows,
                     std::size_t cols)
    {
      gemm_tile<T> accumulators;
      accumulators.fill(T {});
      for (std::size_t i {0}; i < rows; ++i) {
        std::copy_n(c + (i * ldc), cols, accumulators.data() + (i * nr));
      }
      gemm_kernel(kc, ap, bp, ldb, accumulators);
      for (std::size_t i {0}; i < rows; ++i) {
        std::copy_n(accumulators.data() + (i * nr), cols, c + (i * ldc));
      }
    }
  };
 
#ifdef NUMERICAL_NEON_COMPILED
  template <typename T>
  struct neon_tile_kernel : vector_tile_kernel<T>
  {
    NUMERICAL_ALWAYS_INLINE static void tile(std::size_t kc,
                                             const T*    ap,
                                             const T*    bp,
                                             std::size_t ldb,
                                             T*          c,
                                             std::size_t ldc,
                                             std::size_t rows,
                                             std::size_t cols) noexcept
    {
      constexpr std::size_t nr           {vector_tile_kernel<T>::nr};
      gemm_tile<T>          accumulators {};
      for (std::size_t i {0}; i < rows; ++i) {
        std::copy_n(c + (i * ldc), cols, accumulators.data() + (i * nr));
      }
      neon_kernel(kc, ap, bp, ldb, accumulators);
      for (std::size_t i {0}; i < rows; ++i) {
        std::copy_n(accumulators.data() + (i * nr), cols, c + (i * ldc));
      }
    }
  };
#endif
 
#ifdef NUMERICAL_AMX_COMPILED
  struct amx_tile_kernel
  {
    static constexpr std::size_t mr {amx_mr}; 
    static constexpr std::size_t nr {amx_nr}; 
    static constexpr std::size_t depth      {512};
    static constexpr std::size_t block_rows {gemm_mc};             
    static constexpr std::size_t b_in_place_rows {block_rows};
    static constexpr std::size_t speedup         {gemm_kernel_speedup}; 
    static constexpr std::size_t serial_below {std::size_t {1} << 24};
    using session = amx_session; 
 
    static constexpr bool reads_b_in_place(std::size_t ldb) noexcept
    {
      return ldb <= amx_in_place_stride;
    }
 
    static void pack_a(const double* a,
                       std::size_t   lda,
                       std::size_t   rows,
                       std::size_t   kc,
                       double*       out) noexcept
    {
      amx_pack_a(a, lda, rows, kc, out);
    }
 
    static void tile(std::size_t   kc,
                     const double* ap,
                     const double* bp,
                     std::size_t   ldb,
                     double*       c,
                     std::size_t   ldc,
                     std::size_t   rows,
                     std::size_t   cols) noexcept
    {
      amx_tile_f64(kc, ap, bp, ldb, c, ldc, rows, cols);
    }
  };
 
  struct amx_tile_kernel_f32
  {
    static constexpr std::size_t mr              {amx_f32_side};
    static constexpr std::size_t nr              {amx_f32_side};
    static constexpr std::size_t depth           {512};
    static constexpr std::size_t block_rows      {gemm_mc};
    static constexpr std::size_t b_in_place_rows {block_rows};
    static constexpr std::size_t speedup         {gemm_kernel_speedup};
    static constexpr std::size_t serial_below    {std::size_t {1} << 24};
    using session = amx_session;
 
    static constexpr bool reads_b_in_place(std::size_t ldb) noexcept
    {
      return ldb <= 2 * amx_in_place_stride;
    }
 
    static void pack_a(const float* a,
                       std::size_t  lda,
                       std::size_t  rows,
                       std::size_t  kc,
                       float*       out) noexcept
    {
      amx_pack_a_f32(a, lda, rows, kc, out);
    }
 
    static void tile(std::size_t  kc,
                     const float* ap,
                     const float* bp,
                     std::size_t  ldb,
                     float*       c,
                     std::size_t  ldc,
                     std::size_t  rows,
                     std::size_t  cols) noexcept
    {
      amx_tile_f32(kc, ap, bp, ldb, c, ldc, rows, cols);
    }
  };
#endif
 
#ifdef NUMERICAL_X86_COMPILED
  inline constexpr std::size_t narrow_least_depth {8};
#else
  inline constexpr std::size_t narrow_least_depth {2};
#endif
 
  inline constexpr std::size_t narrow_columns {8};
 
  template <typename T>
  [[nodiscard]] inline T rounded_apart_add(T sum,
                                           T a,
                                           T b) noexcept
  {
#if defined(__GNUC__) && !defined(__clang__)
    const T product {a * b};
    return sum + product;
#else
    return sum + (a * b);
#endif
  }
 
  template <typename R>
  [[nodiscard]] inline std::complex<R> rounded_apart_add(std::complex<R> sum,
                                                         std::complex<R> a,
                                                         std::complex<R> b) noexcept
  {
    return sum + operator_times(a, b);
  }
 
  inline constexpr std::size_t narrow_group {4};
 
  enum class term_rounding : unsigned char
  {
    compiler,  
    fused,     
    x86_fused, 
  };
 
  template <term_rounding Mode>
  using rounding_is = std::integral_constant<term_rounding, Mode>;
 
  template <typename T>
  [[nodiscard]] inline T rounded_term(rounding_is<term_rounding::compiler> /*mode*/,
                                      T sum,
                                      T a,
                                      T b) noexcept
  {
    return rounded_apart_add(sum, a, b);
  }
 
  template <typename T>
  [[nodiscard]] inline T rounded_term(rounding_is<term_rounding::fused> /*mode*/,
                                      T sum,
                                      T a,
                                      T b) noexcept
  {
    return std::fma(a, b, sum);
  }
 
#ifdef NUMERICAL_X86_COMPILED
  template <std::size_t R, std::size_t N, typename T>
  NUMERICAL_OUT_OF_LINE NUMERICAL_TARGET("avx2,fma") void narrow_block_x86(std::size_t k,
                                                                           const T*    a,
                                                                           std::size_t lda,
                                                                           const T*    b,
                                                                           std::size_t ldb,
                                                                           T*          c,
                                                                           std::size_t ldc) noexcept
  {
    std::array<T, R * N> sums {};
    NUMERICAL_UNROLL_BLOCK
    for (std::size_t e {0}; e < R * N; ++e) {
      sums[e] = c[((e / N) * ldc) + (e % N)];
    }
    for (std::size_t p {0}; p < k; ++p) {
      NUMERICAL_UNROLL_BLOCK
      for (std::size_t r {0}; r < R; ++r) {
        const T ap {a[(r * lda) + p]};
        NUMERICAL_UNROLL_BLOCK
        for (std::size_t j {0}; j < N; ++j) {
          sums[(r * N) + j] = std::fma(ap, b[(p * ldb) + j], sums[(r * N) + j]);
        }
      }
    }
    NUMERICAL_UNROLL_BLOCK
    for (std::size_t e {0}; e < R * N; ++e) {
      c[((e / N) * ldc) + (e % N)] = sums[e];
    }
  }
 
  template <typename T>
  NUMERICAL_TARGET("avx2,fma") void plain_row_x86(std::size_t n,
                                                  std::size_t k,
                                                  const T*    a,
                                                  const T*    b,
                                                  std::size_t ldb,
                                                  T*          c) noexcept
  {
    for (std::size_t p {0}; p < k; ++p) {
      const T ap {a[p]};
      for (std::size_t j {0}; j < n; ++j) {
        c[j] = std::fma(ap, b[(p * ldb) + j], c[j]);
      }
    }
  }
 
#endif
 
  template <std::size_t R, std::size_t N, term_rounding Mode, typename T>
  NUMERICAL_OUT_OF_LINE void narrow_block(std::size_t k,
                                          const T*    a,
                                          std::size_t lda,
                                          const T*    b,
                                          std::size_t ldb,
                                          T*          c,
                                          std::size_t ldc) noexcept
  {
#ifdef NUMERICAL_X86_COMPILED
    if constexpr (Mode == term_rounding::x86_fused) {
      narrow_block_x86<R, N>(k, a, lda, b, ldb, c, ldc);
    }
    else
#endif
    {
      std::array<T, R * N> sums {};
      NUMERICAL_UNROLL_BLOCK
      for (std::size_t e {0}; e < R * N; ++e) {
        sums[e] = c[((e / N) * ldc) + (e % N)];
      }
      NUMERICAL_SCALAR_DEPTH
      for (std::size_t p {0}; p < k; ++p) {
        NUMERICAL_UNROLL_BLOCK
        for (std::size_t r {0}; r < R; ++r) {
          const T ap {a[(r * lda) + p]};
          NUMERICAL_UNROLL_BLOCK
          for (std::size_t j {0}; j < N; ++j) {
            sums[(r * N) + j] = rounded_term(rounding_is<Mode> {}, sums[(r * N) + j], ap, b[(p * ldb) + j]);
          }
        }
      }
      NUMERICAL_UNROLL_BLOCK
      for (std::size_t e {0}; e < R * N; ++e) {
        c[((e / N) * ldc) + (e % N)] = sums[e];
      }
    }
  }
 
  template <std::size_t R, term_rounding Mode, typename T>
  void narrow_rows(std::size_t n,
                   std::size_t k,
                   const T*    a,
                   std::size_t lda,
                   const T*    b,
                   std::size_t ldb,
                   T*          c,
                   std::size_t ldc) noexcept
  {
    static_assert(narrow_columns == 8, "the switch below names each narrower last block");
    std::size_t j {0};
    for (; j + narrow_columns <= n; j += narrow_columns) {
      narrow_block<R, narrow_columns, Mode>(k, a, lda, b + j, ldb, c + j, ldc);
    }
    switch (n - j) {
      case 1: narrow_block<R, 1, Mode>(k, a, lda, b + j, ldb, c + j, ldc); break;
      case 2: narrow_block<R, 2, Mode>(k, a, lda, b + j, ldb, c + j, ldc); break;
      case 3: narrow_block<R, 3, Mode>(k, a, lda, b + j, ldb, c + j, ldc); break;
      case 4: narrow_block<R, 4, Mode>(k, a, lda, b + j, ldb, c + j, ldc); break;
      case 5: narrow_block<R, 5, Mode>(k, a, lda, b + j, ldb, c + j, ldc); break;
      case 6: narrow_block<R, 6, Mode>(k, a, lda, b + j, ldb, c + j, ldc); break;
      case 7: narrow_block<R, 7, Mode>(k, a, lda, b + j, ldb, c + j, ldc); break;
      default: break; // a whole number of blocks
    }
  }
 
  template <term_rounding Mode, typename T>
  void gemm_rows_rounded(std::size_t m,
                         std::size_t n,
                         std::size_t k,
                         const T*    a,
                         std::size_t lda,
                         const T*    b,
                         std::size_t ldb,
                         T*          c,
                         std::size_t ldc)
  {
    if (k >= (n <= narrow_columns ? std::size_t {2} : narrow_least_depth)) {
      parallel_for((m + narrow_group - 1) / narrow_group, narrow_group * k * n / vector_speedup, [&](std::size_t g) {
                     const std::size_t first {g * narrow_group};
                     if (first + narrow_group <= m) {
                       narrow_rows<narrow_group, Mode>(n, k, a + (first * lda), lda, b, ldb, c + (first * ldc), ldc);
                       return;
                     }
                     for (std::size_t i {first}; i < m; ++i) {
                       narrow_rows<1, Mode>(n, k, a + (i * lda), lda, b, ldb, c + (i * ldc), ldc);
                     }
                   });
      return;
    }
    parallel_for(m, k * n / vector_speedup, [&](std::size_t i) {
#ifdef NUMERICAL_X86_COMPILED
                   if constexpr (Mode == term_rounding::x86_fused) {
                     plain_row_x86(n, k, a + (i * lda), b, ldb, c + (i * ldc));
                   }
                   else
#endif
                   {
                     for (std::size_t p {0}; p < k; ++p) {
                       const T aip {a[(i * lda) + p]};
                       for (std::size_t j {0}; j < n; ++j) {
                         c[(i * ldc) + j] = rounded_term(rounding_is<Mode> {}, c[(i * ldc) + j], aip, b[(p * ldb) + j]);
                       }
                     }
                   }
                 });
  }
 
  [[nodiscard]] inline bool rows_fuse() noexcept
  {
#ifdef NUMERICAL_NEON_COMPILED
    return gemm_path_cap().load(std::memory_order_relaxed) != gemm_path::portable;
#else
    return false;
#endif
  }
 
  template <typename T>
  void gemm_rows(std::size_t m,
                 std::size_t n,
                 std::size_t k,
                 const T*    a,
                 std::size_t lda,
                 const T*    b,
                 std::size_t ldb,
                 T*          c,
                 std::size_t ldc)
  {
#ifdef NUMERICAL_X86_COMPILED
    if constexpr (std::is_same_v<T, double> || std::is_same_v<T, float>) {
      const gemm_path path {gemm_path_current()};
      if (path == gemm_path::avx512 || path == gemm_path::avx2) {
        gemm_rows_rounded<term_rounding::x86_fused>(m, n, k, a, lda, b, ldb, c, ldc);
        return;
      }
    }
#endif
    if constexpr (std::is_floating_point_v<T>) {
      if (rows_fuse()) {
        gemm_rows_rounded<term_rounding::fused>(m, n, k, a, lda, b, ldb, c, ldc);
        return;
      }
    }
    gemm_rows_rounded<term_rounding::compiler>(m, n, k, a, lda, b, ldb, c, ldc);
  }
 
  inline constexpr std::size_t packed_line {64};
 
  template <typename T>
  [[nodiscard]] T* line_aligned(T* at) noexcept
  {
    const std::size_t miss {static_cast<std::size_t>(std::bit_cast<std::uintptr_t>(at) % packed_line)};
    return at + (((packed_line - miss) % packed_line) / sizeof(T));
  }
 
  template <typename T>
  class packing_storage
  {
  public:
 
    static constexpr std::size_t kept_limit {(std::size_t {32} << 20U) / sizeof(T)};
 
    packing_storage(std::size_t a_count,
                    std::size_t b_count)
      : storage_ {grown(std::exchange(kept(), buffer<T> {}), a_count + padding + b_count + padding)},
        a_ {aligned(storage_.data())},
        b_ {aligned(storage_.data() + a_count + padding)}
    {}
 
    packing_storage(const packing_storage&)            = delete;
    packing_storage& operator=(const packing_storage&) = delete;
    packing_storage(packing_storage&&)                 = delete;
    packing_storage& operator=(packing_storage&&)      = delete;
 
    ~packing_storage()
    {
      kept() = storage_.size() <= kept_limit ? std::move(storage_) : buffer<T> {};
    }
 
    [[nodiscard]] T* a() const noexcept
    {
      return a_;
    }
 
    [[nodiscard]] T* b() const noexcept
    {
      return b_;
    }
 
    [[nodiscard]] static buffer<T>& kept() noexcept
    {
      thread_local buffer<T> storage;
      return storage;
    }
 
  private:
 
    static constexpr std::size_t padding {(packed_line / sizeof(T)) + 1}; 
 
    [[nodiscard]] static buffer<T> grown(buffer<T>      storage,
                                         std::size_t    count)
    {
      storage.resize(std::max(storage.size(), count), T {});
      return storage;
    }
 
    [[nodiscard]] static T* aligned(T* at) noexcept
    {
      return line_aligned(at);
    }
 
    buffer<T>      storage_;
    T*             a_;
    T*             b_;
  };
 
  template <typename Kernel>
  [[nodiscard]] constexpr bool gemm_shares(std::size_t m,
                                           std::size_t n,
                                           std::size_t k) noexcept
  {
    return m * n * k >= Kernel::serial_below;
  }
 
  template <typename Kernel>
  inline constexpr bool prefetches_next_row_of_c {requires(const double* c, std::size_t ldc) {
                                                    Kernel::prefetch_rows(c, ldc);
                                                  }};
 
  template <typename Kernel,
            typename T>
  void prefetch_next_row(const T* /*c*/,
                         std::size_t /*ldc*/) noexcept
  {
  }
 
  template <typename Kernel,
            typename T>
  requires prefetches_next_row_of_c<Kernel>
  void prefetch_next_row(const T*    c,
                         std::size_t ldc) noexcept
  {
    Kernel::prefetch_rows(c, ldc);
  }
 
  template <typename Kernel>
  [[nodiscard]] constexpr std::size_t task_rows_cap(std::size_t /*kc*/) noexcept
  {
    return Kernel::block_rows;
  }
 
  template <typename Kernel>
  requires requires(std::size_t kc)
  {
    Kernel::task_rows_for(kc);
  }
 
  [[nodiscard]] constexpr std::size_t task_rows_cap(std::size_t kc) noexcept
  {
    return Kernel::task_rows_for(kc);
  }
 
  template <typename Kernel,
            typename T,
            typename B = b_rows,
            typename A = a_rows>
  void gemm_packed(std::size_t m,
                   std::size_t n,
                   std::size_t k,
                   const T*    a,
                   std::size_t lda,
                   const T*    b,
                   std::size_t ldb,
                   T*          c,
                   std::size_t ldc,
                   const T*    a_packed = nullptr)
  {
    constexpr std::size_t    mr     {Kernel::mr};
    constexpr std::size_t    nr     {Kernel::nr};
    constexpr std::size_t    nc_max {nr * 128};
    const std::size_t        a_rows {((m + mr) - 1) / mr * mr};
    const std::size_t        b_cols {((std::min(nc_max, n) + nr) - 1) / nr * nr};
    const packing_storage<T> packed {a_packed != nullptr ? 0 : a_rows * std::min(Kernel::depth, k), b_cols * std::min(Kernel::depth, k)};
    const std::size_t        a_blocks {a_packed != nullptr ? 0 : ((m + Kernel::block_rows) - 1) / Kernel::block_rows};
    const T* const           a_tiles  {a_packed != nullptr ? a_packed : packed.a()};
    const std::size_t        shared   {gemm_shares<Kernel>(m, n, k) ? 1U : 0U}; // a weight of zero starts no thread
    // B's whole slivers unpacked: rows close enough, or few enough rows of A whatever B's stride
    // (Kernel::b_in_place_rows), where the kernel reads B in place at all.
    const bool               in_place {B::template in_place<Kernel>(ldb) || (m <= Kernel::b_in_place_rows && B::template in_place<Kernel>(1))};
    for (std::size_t jc {0}; jc < n; jc += nc_max) {
      const std::size_t nc      {std::min(nc_max, n - jc)};
      const std::size_t slivers {((nc + nr) - 1) / nr};
      for (std::size_t pc {0}; pc < k; pc += Kernel::depth) {
        const std::size_t kc         {std::min(Kernel::depth, k - pc)};
        const std::size_t groups     {((slivers + gemm_slivers_per_task) - 1) / gemm_slivers_per_task};
        // One pass packs both: B's groups of slivers first, then A's row blocks.
        // The packing weighed as the tile's work is (\p Kernel's speedup): a small product then packs in parallel too
        // (eight AVX2 cores: order 128 40 → 30 µs, 192 92 → 62).
        parallel_for(groups + a_blocks, shared * Kernel::block_rows * kc * gemm_kernel_speedup / Kernel::speedup, [&](std::size_t piece) {
                       if (piece < groups) {
                         const std::size_t jr    {piece * gemm_slivers_per_task * nr};
                         const std::size_t width {std::min(gemm_slivers_per_task * nr, nc - jr)};
                         const std::size_t whole {in_place ? width - (width % nr) : 0}; // read in place
                         if (whole < width) { // a partial sliver is packed even so: its outside zero
                           B::template pack<nr>(b, ldb, pc, jc + jr + whole, kc, width - whole, packed.b() + ((jr + whole) * kc));
                         }
                       }
                       else {
                         const std::size_t ic {(piece - groups) * Kernel::block_rows};
                         A::template pack<Kernel>(a, lda, ic, pc, std::min(Kernel::block_rows, m - ic), kc, packed.a() + (ic * kc));
                       }
                     });
        // A task covers task_rows rows: a whole block, or fewer whole tiles when blocks would leave a thread fewer
        // than tasks_per_thread tasks (order 256 on twenty x86-64 threads: twelve tasks took 116 µs, 44 took 104).
        const std::size_t tiles      {((m + mr) - 1) / mr};
        const std::size_t row_tasks  {(((tasks_per_thread * max_threads()) + groups) - 1) / groups};
        const std::size_t task_rows  {std::clamp(((tiles + row_tasks) - 1) / row_tasks, std::size_t {1}, task_rows_cap<Kernel>(kc) / mr) * mr};
        const std::size_t task_count {((m + task_rows) - 1) / task_rows};
        parallel_for(task_count * groups, shared * task_rows * gemm_slivers_per_task * nr * kc / Kernel::speedup, [&](std::size_t task) {
                       const std::size_t                   ic     {(task / groups) * task_rows};
                       const std::size_t                   first  {(task % groups) * gemm_slivers_per_task};
                       const std::size_t                   last   {std::min(first + gemm_slivers_per_task, slivers)};
                       const std::size_t                   ic_end {std::min(ic + task_rows, m)};
                       [[maybe_unused]] const typename Kernel::session on {};
                       // A sliver (the operand each step broadcasts) stays in the L1 cache across the task's B slivers,
                       // which stream from the L2, as BLIS orders its loops: B sliver outer, its 32 KiB filled the L1
                       // (a product of order 2048 on one AVX-512 core: 72 GFLOP/s, 77 in this order; AVX2 and AMX gain too).
                       for (std::size_t ir {ic}; ir < ic_end; ir += mr) {
                         for (std::size_t s {first}; s < last; ++s) {
                           const std::size_t jr {s * nr};
                           const bool direct {in_place && jr + nr <= nc};
                           if (ir + mr < ic_end) { // the next row's tile of C, while this one computes
                             prefetch_next_row<Kernel>(c + ((ir + mr) * ldc) + jc + jr, ldc);
                           }
                           Kernel::tile(kc, a_tiles + (ir * kc), direct ? b + (pc * ldb) + jc + jr : packed.b() + (jr * kc),
                                        direct ? ldb : nr, c + (ir * ldc) + jc + jr, ldc, std::min(mr, m - ir), std::min(nr, nc - jr));
                         }
                       }
                     });
      }
    }
  }
 
  template <typename Kernel,
            typename T,
            typename B = b_rows,
            typename A = a_rows>
  void gemm_packed_triangle(std::size_t n,
                            std::size_t k,
                            const T*    a,
                            std::size_t lda,
                            const T*    b,
                            std::size_t ldb,
                            T*          c,
                            std::size_t ldc,
                            bool        lower)
  {
    constexpr std::size_t    mr       {Kernel::mr};
    constexpr std::size_t    nr       {Kernel::nr};
    constexpr std::size_t    nc_max   {nr * 128};
    const std::size_t        a_rows   {((n + mr) - 1) / mr * mr};
    const std::size_t        b_cols   {((std::min(nc_max, n) + nr) - 1) / nr * nr};
    const packing_storage<T> packed   {a_rows * std::min(Kernel::depth, k), b_cols * std::min(Kernel::depth, k)};
    const std::size_t        a_blocks {((n + Kernel::block_rows) - 1) / Kernel::block_rows};
    const std::size_t        shared   {gemm_shares<Kernel>(n, n / 2, k) ? 1U : 0U};
    const bool               in_place {B::template in_place<Kernel>(ldb)};
    // Whether rows [r0, r1) and columns [j0, j1) of C meet the triangle, and whether they lie wholly in it.
    const auto meets  = [lower](std::size_t r0, std::size_t r1, std::size_t j0, std::size_t j1) { return lower ? j0 < r1 : r0 < j1; };
    const auto inside = [lower](std::size_t r0, std::size_t r1, std::size_t j0, std::size_t j1) { return lower ? j1 <= r0 + 1 : r1 <= j0 + 1; };
    for (std::size_t jc {0}; jc < n; jc += nc_max) {
      const std::size_t nc      {std::min(nc_max, n - jc)};
      const std::size_t slivers {((nc + nr) - 1) / nr};
      for (std::size_t pc {0}; pc < k; pc += Kernel::depth) {
        const std::size_t kc     {std::min(Kernel::depth, k - pc)};
        const std::size_t groups {((slivers + gemm_slivers_per_task) - 1) / gemm_slivers_per_task};
        parallel_for(groups + a_blocks, shared * Kernel::block_rows * kc * gemm_kernel_speedup / Kernel::speedup, [&](std::size_t piece) {
                       if (piece < groups) {
                         const std::size_t jr    {piece * gemm_slivers_per_task * nr};
                         const std::size_t width {std::min(gemm_slivers_per_task * nr, nc - jr)};
                         const std::size_t whole {in_place ? width - (width % nr) : 0}; // read in place
                         if (whole < width) { // a partial sliver is packed even so: its outside zero
                           B::template pack<nr>(b, ldb, pc, jc + jr + whole, kc, width - whole, packed.b() + ((jr + whole) * kc));
                         }
                         return;
                       }
                       const std::size_t ic   {(piece - groups) * Kernel::block_rows};
                       const std::size_t rows {std::min(Kernel::block_rows, n - ic)};
                       if (meets(ic, ic + rows, jc, jc + nc)) {
                         A::template pack<Kernel>(a, lda, ic, pc, rows, kc, packed.a() + (ic * kc));
                       }
                     });
        const std::size_t tiles      {((n + mr) - 1) / mr};
        const std::size_t row_tasks  {(((tasks_per_thread * max_threads()) + groups) - 1) / groups};
        const std::size_t task_rows  {std::clamp(((tiles + row_tasks) - 1) / row_tasks, std::size_t {1}, task_rows_cap<Kernel>(kc) / mr) * mr};
        const std::size_t task_count {((n + task_rows) - 1) / task_rows};
        parallel_for(task_count * groups, shared * task_rows * gemm_slivers_per_task * nr * kc / (2 * Kernel::speedup), [&](std::size_t task) {
                       const std::size_t                   ic     {(task / groups) * task_rows};
                       const std::size_t                   first  {(task % groups) * gemm_slivers_per_task};
                       const std::size_t                   last   {std::min(first + gemm_slivers_per_task, slivers)};
                       const std::size_t                   ic_end {std::min(ic + task_rows, n)};
                       [[maybe_unused]] const typename Kernel::session on {};
                       for (std::size_t ir {ic}; ir < ic_end; ir += mr) {
                         const std::size_t rows {std::min(mr, n - ir)};
                         for (std::size_t s {first}; s < last; ++s) {
                           const std::size_t jr   {s * nr};
                           const std::size_t j0   {jc + jr};
                           const std::size_t cols {std::min(nr, nc - jr)};
                           if (!meets(ir, ir + rows, j0, j0 + cols)) {
                             continue;
                           }
                           const bool     direct {in_place && jr + nr <= nc};
                           const T* const bp     {direct ? b + (pc * ldb) + j0 : packed.b() + (jr * kc)};
                           const std::size_t ldbp {direct ? ldb : nr};
                           T* const       block  {c + (ir * ldc) + j0};
                           if (inside(ir, ir + rows, j0, j0 + cols)) {
                             Kernel::tile(kc, packed.a() + (ir * kc), bp, ldbp, block, ldc, rows, cols);
                             continue;
                           }
                           std::array<T, mr * nr> copy {};
                           for (std::size_t i {0}; i < rows; ++i) {
                             std::copy_n(block + (i * ldc), cols, copy.data() + (i * nr));
                           }
                           Kernel::tile(kc, packed.a() + (ir * kc), bp, ldbp, copy.data(), nr, rows, cols);
                           for (std::size_t i {0}; i < rows; ++i) {
                             for (std::size_t j {0}; j < cols; ++j) {
                               if (lower ? j0 + j <= ir + i : j0 + j >= ir + i) {
                                 block[(i * ldc) + j] = copy[(i * nr) + j];
                               }
                             }
                           }
                         }
                       }
                     });
      }
    }
  }
 
#ifdef NUMERICAL_X86_COMPILED
  struct avx512_tile_kernel
  {
    static constexpr std::size_t mr              {avx512_mr}; 
    static constexpr std::size_t nr              {avx512_nr}; 
    static constexpr std::size_t depth           {gemm_kc};   
    static constexpr std::size_t block_rows      {96};        
    static constexpr std::size_t b_in_place_rows {mr};        
    static constexpr std::size_t speedup      {gemm_kernel_speedup / 4};
    static constexpr std::size_t serial_below {1}; 
 
    struct session
    {};
 
    static constexpr bool reads_b_in_place(std::size_t ldb) noexcept
    {
      return ldb <= nr;
    }
 
    static void pack_a(const double* a,
                       std::size_t   lda,
                       std::size_t   rows,
                       std::size_t   kc,
                       double*       out)
    {
      gemm_pack_a<mr>(a, lda, rows, kc, out);
    }
 
    static void tile(std::size_t   kc,
                     const double* ap,
                     const double* bp,
                     std::size_t   ldb,
                     double*       c,
                     std::size_t   ldc,
                     std::size_t   rows,
                     std::size_t   cols) noexcept
    {
      avx512_tile_f64(kc, ap, bp, ldb, c, ldc, rows, cols);
    }
  };
 
  struct avx2_tile_kernel
  {
    static constexpr std::size_t mr              {avx2_mr}; 
    static constexpr std::size_t nr              {avx2_nr}; 
    static constexpr std::size_t depth           {gemm_kc}; 
    static constexpr std::size_t block_rows      {72};      
    static constexpr std::size_t b_in_place_rows {mr};      
 
    [[nodiscard]] static constexpr std::size_t task_rows_for(std::size_t kc) noexcept
    {
      return kc <= 128 ? 36 : block_rows;
    }
 
    static void prefetch_rows(const double* c,
                              std::size_t   ldc) noexcept
    {
      for (std::size_t r {0}; r < mr; ++r) {
        _mm_prefetch(reinterpret_cast<const char*>(c + (r * ldc)), _MM_HINT_T0); // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the intrinsic takes a byte address
      }
    }
 
    static constexpr std::size_t speedup      {gemm_kernel_speedup / 4};
    static constexpr std::size_t serial_below {1}; 
 
    struct session
    {};
 
    static constexpr bool reads_b_in_place(std::size_t ldb) noexcept
    {
      return ldb <= nr;
    }
 
    static void pack_a(const double* a,
                       std::size_t   lda,
                       std::size_t   rows,
                       std::size_t   kc,
                       double*       out)
    {
      gemm_pack_a<mr>(a, lda, rows, kc, out);
    }
 
    static void tile(std::size_t   kc,
                     const double* ap,
                     const double* bp,
                     std::size_t   ldb,
                     double*       c,
                     std::size_t   ldc,
                     std::size_t   rows,
                     std::size_t   cols) noexcept
    {
      avx2_tile_f64(kc, ap, bp, ldb, c, ldc, rows, cols);
    }
  };
  struct avx512_tile_kernel_f32
  {
    static constexpr std::size_t mr              {avx512_mr};
    static constexpr std::size_t nr              {avx512_nr_f32};
    static constexpr std::size_t depth           {gemm_kc};
    static constexpr std::size_t block_rows      {96};
    static constexpr std::size_t b_in_place_rows {mr};
    static constexpr std::size_t speedup         {gemm_kernel_speedup / 4};
    static constexpr std::size_t serial_below    {1};
 
    struct session
    {};
 
    static constexpr bool reads_b_in_place(std::size_t ldb) noexcept
    {
      return ldb <= nr;
    }
 
    static void pack_a(const float* a,
                       std::size_t  lda,
                       std::size_t  rows,
                       std::size_t  kc,
                       float*       out)
    {
      gemm_pack_a<mr>(a, lda, rows, kc, out);
    }
 
    static void tile(std::size_t  kc,
                     const float* ap,
                     const float* bp,
                     std::size_t  ldb,
                     float*       c,
                     std::size_t  ldc,
                     std::size_t  rows,
                     std::size_t  cols) noexcept
    {
      avx512_tile_f32(kc, ap, bp, ldb, c, ldc, rows, cols);
    }
  };
 
  struct avx2_tile_kernel_f32
  {
    static constexpr std::size_t mr              {avx2_mr};
    static constexpr std::size_t nr              {avx2_nr_f32};
    static constexpr std::size_t depth           {gemm_kc};
    static constexpr std::size_t block_rows      {72};
    static constexpr std::size_t b_in_place_rows {mr};
    static constexpr std::size_t speedup         {gemm_kernel_speedup / 4};
    static constexpr std::size_t serial_below    {1};
 
    [[nodiscard]] static constexpr std::size_t task_rows_for(std::size_t kc) noexcept
    {
      return kc <= 128 ? 36 : block_rows;
    }
 
    struct session
    {};
 
    static constexpr bool reads_b_in_place(std::size_t ldb) noexcept
    {
      return ldb <= nr;
    }
 
    static void pack_a(const float* a,
                       std::size_t  lda,
                       std::size_t  rows,
                       std::size_t  kc,
                       float*       out)
    {
      gemm_pack_a<mr>(a, lda, rows, kc, out);
    }
 
    static void tile(std::size_t  kc,
                     const float* ap,
                     const float* bp,
                     std::size_t  ldb,
                     float*       c,
                     std::size_t  ldc,
                     std::size_t  rows,
                     std::size_t  cols) noexcept
    {
      avx2_tile_f32(kc, ap, bp, ldb, c, ldc, rows, cols);
    }
  };
 
  struct avx512_complex_tile_kernel
  {
    static constexpr std::size_t mr              {avx512_complex_mr};        
    static constexpr std::size_t nr              {avx512_complex_nr};        
    static constexpr std::size_t depth           {gemm_kc / 2};              
    static constexpr std::size_t block_rows      {96};                       
    static constexpr std::size_t b_in_place_rows {mr};                       
    static constexpr std::size_t speedup         {gemm_kernel_speedup / 16}; 
    static constexpr std::size_t serial_below    {1};                        
 
    struct session
    {};
 
    static constexpr bool reads_b_in_place(std::size_t ldb) noexcept
    {
      return ldb <= nr;
    }
 
    static void pack_a(const std::complex<double>* a,
                       std::size_t                 lda,
                       std::size_t                 rows,
                       std::size_t                 kc,
                       std::complex<double>*       out)
    {
      gemm_pack_a<mr>(a, lda, rows, kc, out);
    }
 
    static void tile(std::size_t                 kc,
                     const std::complex<double>* ap,
                     const std::complex<double>* bp,
                     std::size_t                 ldb,
                     std::complex<double>*       c,
                     std::size_t                 ldc,
                     std::size_t                 rows,
                     std::size_t                 cols) noexcept
    {
      avx512_tile_c64(kc, ap, bp, ldb, c, ldc, rows, cols);
    }
  };
 
  struct avx2_complex_tile_kernel
  {
    static constexpr std::size_t mr              {avx2_complex_mr};          
    static constexpr std::size_t nr              {avx2_complex_nr};          
    static constexpr std::size_t depth           {gemm_kc / 2};              
    static constexpr std::size_t block_rows      {72};                       
    static constexpr std::size_t b_in_place_rows {mr};                       
    static constexpr std::size_t speedup         {gemm_kernel_speedup / 16}; 
    static constexpr std::size_t serial_below    {1};                        
 
    struct session
    {};
 
    static constexpr bool reads_b_in_place(std::size_t ldb) noexcept
    {
      return ldb <= nr;
    }
 
    static void pack_a(const std::complex<double>* a,
                       std::size_t                 lda,
                       std::size_t                 rows,
                       std::size_t                 kc,
                       std::complex<double>*       out)
    {
      gemm_pack_a<mr>(a, lda, rows, kc, out);
    }
 
    static void tile(std::size_t                 kc,
                     const std::complex<double>* ap,
                     const std::complex<double>* bp,
                     std::size_t                 ldb,
                     std::complex<double>*       c,
                     std::size_t                 ldc,
                     std::size_t                 rows,
                     std::size_t                 cols) noexcept
    {
      avx2_tile_c64(kc, ap, bp, ldb, c, ldc, rows, cols);
    }
  };
  struct avx512_complex_tile_kernel_f32
  {
    static constexpr std::size_t mr              {avx512_complex_mr};
    static constexpr std::size_t nr              {avx512_complex_nr_f32};
    static constexpr std::size_t depth           {gemm_kc}; // an element is two `float`, the size of a `double`
    static constexpr std::size_t block_rows      {96};
    static constexpr std::size_t b_in_place_rows {mr};
    static constexpr std::size_t speedup         {gemm_kernel_speedup / 16};
    static constexpr std::size_t serial_below    {1};
 
    struct session
    {};
 
    static constexpr bool reads_b_in_place(std::size_t ldb) noexcept
    {
      return ldb <= nr;
    }
 
    static void pack_a(const std::complex<float>* a,
                       std::size_t                lda,
                       std::size_t                rows,
                       std::size_t                kc,
                       std::complex<float>*       out)
    {
      gemm_pack_a<mr>(a, lda, rows, kc, out);
    }
 
    static void tile(std::size_t                kc,
                     const std::complex<float>* ap,
                     const std::complex<float>* bp,
                     std::size_t                ldb,
                     std::complex<float>*       c,
                     std::size_t                ldc,
                     std::size_t                rows,
                     std::size_t                cols) noexcept
    {
      avx512_tile_c32(kc, ap, bp, ldb, c, ldc, rows, cols);
    }
  };
 
  struct avx2_complex_tile_kernel_f32
  {
    static constexpr std::size_t mr              {avx2_complex_mr};
    static constexpr std::size_t nr              {avx2_complex_nr_f32};
    static constexpr std::size_t depth           {gemm_kc}; // an element is two `float`, the size of a `double`
    static constexpr std::size_t block_rows      {72};
    static constexpr std::size_t b_in_place_rows {mr};
    static constexpr std::size_t speedup         {gemm_kernel_speedup / 16};
    static constexpr std::size_t serial_below    {1};
 
    struct session
    {};
 
    static constexpr bool reads_b_in_place(std::size_t ldb) noexcept
    {
      return ldb <= nr;
    }
 
    static void pack_a(const std::complex<float>* a,
                       std::size_t                lda,
                       std::size_t                rows,
                       std::size_t                kc,
                       std::complex<float>*       out)
    {
      gemm_pack_a<mr>(a, lda, rows, kc, out);
    }
 
    static void tile(std::size_t                kc,
                     const std::complex<float>* ap,
                     const std::complex<float>* bp,
                     std::size_t                ldb,
                     std::complex<float>*       c,
                     std::size_t                ldc,
                     std::size_t                rows,
                     std::size_t                cols) noexcept
    {
      avx2_tile_c32(kc, ap, bp, ldb, c, ldc, rows, cols);
    }
  };
 
#endif
 
  template <typename T>
  [[nodiscard]] bool gemm_elsewhere(std::size_t m,
                                    std::size_t n,
                                    std::size_t k,
                                    const T*    a,
                                    std::size_t lda,
                                    const T*    b,
                                    std::size_t ldb,
                                    T*          c,
                                    std::size_t ldc)
  {
    (void)m, (void)n, (void)k, (void)a, (void)lda, (void)b, (void)ldb, (void)c, (void)ldc;
    return false;
  }
 
  template <typename B,
            typename A = a_rows>
  [[nodiscard]] bool gemm_double_elsewhere(std::size_t   m,
                                           std::size_t   n,
                                           std::size_t   k,
                                           const double* a,
                                           std::size_t   lda,
                                           const double* b,
                                           std::size_t   ldb,
                                           double*       c,
                                           std::size_t   ldc)
  {
    const gemm_path path {gemm_path_current()};
#ifdef NUMERICAL_AMX_COMPILED
    if (path == gemm_path::amx && amx_suits(m, n, k)) { // else the NEON tile
      gemm_packed<amx_tile_kernel, double, B, A>(m, n, k, a, lda, b, ldb, c, ldc);
      return true;
    }
#endif
#ifdef NUMERICAL_NEON_COMPILED
    if (path == gemm_path::neon || path == gemm_path::amx) {
      gemm_packed<neon_tile_kernel<double>, double, B, A>(m, n, k, a, lda, b, ldb, c, ldc);
      return true;
    }
#endif
#ifdef NUMERICAL_X86_COMPILED
    if (path == gemm_path::avx512) {
      gemm_packed<avx512_tile_kernel, double, B, A>(m, n, k, a, lda, b, ldb, c, ldc);
      return true;
    }
    if (path == gemm_path::avx2) {
      gemm_packed<avx2_tile_kernel, double, B, A>(m, n, k, a, lda, b, ldb, c, ldc);
      return true;
    }
#endif
    (void)path;
    (void)m, (void)n, (void)k, (void)a, (void)lda, (void)b, (void)ldb, (void)c, (void)ldc;
    return false;
  }
 
  template <typename A,
            typename B>
  void gemm_triangle(std::size_t   n,
                     std::size_t   k,
                     const double* a,
                     std::size_t   lda,
                     const double* b,
                     std::size_t   ldb,
                     double*       c,
                     std::size_t   ldc,
                     bool          lower)
  {
    const gemm_path path {gemm_path_current()};
#ifdef NUMERICAL_AMX_COMPILED
    if (path == gemm_path::amx && amx_suits(n, n, k)) {
      gemm_packed_triangle<amx_tile_kernel, double, B, A>(n, k, a, lda, b, ldb, c, ldc, lower);
      return;
    }
#endif
#ifdef NUMERICAL_NEON_COMPILED
    if (path == gemm_path::neon || path == gemm_path::amx) {
      gemm_packed_triangle<neon_tile_kernel<double>, double, B, A>(n, k, a, lda, b, ldb, c, ldc, lower);
      return;
    }
#endif
#ifdef NUMERICAL_X86_COMPILED
    if (path == gemm_path::avx512) {
      gemm_packed_triangle<avx512_tile_kernel, double, B, A>(n, k, a, lda, b, ldb, c, ldc, lower);
      return;
    }
    if (path == gemm_path::avx2) {
      gemm_packed_triangle<avx2_tile_kernel, double, B, A>(n, k, a, lda, b, ldb, c, ldc, lower);
      return;
    }
#endif
    (void)path;
    gemm_packed_triangle<vector_tile_kernel<double>, double, B, A>(n, k, a, lda, b, ldb, c, ldc, lower);
  }
 
  [[nodiscard]] inline bool gemm_elsewhere(std::size_t   m,
                                           std::size_t   n,
                                           std::size_t   k,
                                           const double* a,
                                           std::size_t   lda,
                                           const double* b,
                                           std::size_t   ldb,
                                           double*       c,
                                           std::size_t   ldc)
  {
    return gemm_double_elsewhere<b_rows>(m, n, k, a, lda, b, ldb, c, ldc);
  }
 
  template <typename B,
            typename A = a_rows>
  [[nodiscard]] bool gemm_float_elsewhere(std::size_t  m,
                                          std::size_t  n,
                                          std::size_t  k,
                                          const float* a,
                                          std::size_t  lda,
                                          const float* b,
                                          std::size_t  ldb,
                                          float*       c,
                                          std::size_t  ldc)
  {
#ifdef NUMERICAL_NEON_COMPILED
    const gemm_path path {gemm_path_current()};
#  ifdef NUMERICAL_AMX_COMPILED
    if (path == gemm_path::amx && amx_suits_f32(m, n, k)) {
      gemm_packed<amx_tile_kernel_f32, float, B, A>(m, n, k, a, lda, b, ldb, c, ldc);
      return true;
    }
#  endif
    if (path == gemm_path::neon || path == gemm_path::amx) {
      gemm_packed<neon_tile_kernel<float>, float, B, A>(m, n, k, a, lda, b, ldb, c, ldc);
      return true;
    }
#endif
#ifdef NUMERICAL_X86_COMPILED
    const gemm_path x86_path {gemm_path_current()};
    if (x86_path == gemm_path::avx512) {
      gemm_packed<avx512_tile_kernel_f32, float, B, A>(m, n, k, a, lda, b, ldb, c, ldc);
      return true;
    }
    if (x86_path == gemm_path::avx2) {
      gemm_packed<avx2_tile_kernel_f32, float, B, A>(m, n, k, a, lda, b, ldb, c, ldc);
      return true;
    }
#endif
    (void)m, (void)n, (void)k, (void)a, (void)lda, (void)b, (void)ldb, (void)c, (void)ldc;
    return false;
  }
 
  [[nodiscard]] inline bool gemm_elsewhere(std::size_t  m,
                                           std::size_t  n,
                                           std::size_t  k,
                                           const float* a,
                                           std::size_t  lda,
                                           const float* b,
                                           std::size_t  ldb,
                                           float*       c,
                                           std::size_t  ldc)
  {
    return gemm_float_elsewhere<b_rows>(m, n, k, a, lda, b, ldb, c, ldc);
  }
 
  template <typename T>
  void gemm(std::size_t m,
            std::size_t n,
            std::size_t k,
            const T*    a,
            std::size_t lda,
            const T*    b,
            std::size_t ldb,
            T*          c,
            std::size_t ldc);
 
  inline constexpr std::size_t split_columns {512};
 
  inline constexpr std::size_t split_depth {256};
 
  inline constexpr std::size_t split_kept {std::size_t {1} << 18};
 
  inline constexpr std::size_t split_parts_rows {64};
 
  [[nodiscard]] inline bool split_by_parts(std::size_t m) noexcept
  {
    [[maybe_unused]] const gemm_path path {gemm_path_current()};
#ifdef NUMERICAL_AMX_COMPILED
    if (path == gemm_path::amx) {
      return m <= split_parts_rows;
    }
#endif
    return m <= gemm_mr;
  }
 
  template <typename R>
  struct split_operands
  {
    buffer<R> a {}; 
    buffer<R> b {}; 
    buffer<R> c {}; 
  };
 
  template <typename R>
  [[nodiscard]] split_operands<R>& split_storage()
  {
    thread_local split_operands<R> kept;
    return kept;
  }
 
  template <split_complex Z,
            typename B = b_rows,
            typename A = a_rows>
  void split_parts_block(std::size_t                 m,
                         std::size_t                 nb,
                         std::size_t                 kb,
                         std::size_t                 pc,
                         std::size_t                 jc,
                         const Z*                    a,
                         std::size_t                 lda,
                         const Z*                    b,
                         std::size_t                 ldb,
                         typename Z::value_type*     ar,
                         typename Z::value_type*     br,
                         typename Z::value_type*     cr)
  {
    using R = typename Z::value_type;
    parallel_for(m, kb, [&](std::size_t i) {
                   R* const row {ar + (i * 3 * kb)};
                   for (std::size_t p {0}; p < kb; ++p) {
                     const Z value {A::at(a, lda, i, pc + p)};
                     row[p]            = value.real();
                     row[kb + p]       = -value.imag();
                     row[(2 * kb) + p] = value.imag();
                   }
                 });
    parallel_for(kb, nb, [&](std::size_t p) {
                   R* const row {br + (p * 2 * nb)};
                   for (std::size_t j {0}; j < nb; ++j) {
                     const Z value {B::at(b, ldb, pc + p, jc + j)};
                     row[j]      = value.real();
                     row[nb + j] = value.imag();
                   }
                 });
    gemm(m, nb, kb, ar, 3 * kb, br, 2 * nb, cr, 2 * nb);                 // Re C += Re A·Re B
    gemm(m, nb, kb, ar + kb, 3 * kb, br + nb, 2 * nb, cr, 2 * nb);       // Re C += −Im A·Im B
    gemm(m, nb, kb, ar, 3 * kb, br + nb, 2 * nb, cr + nb, 2 * nb);       // Im C += Re A·Im B
    gemm(m, nb, kb, ar + (2 * kb), 3 * kb, br, 2 * nb, cr + nb, 2 * nb); // Im C += Im A·Re B
  }
 
  template <split_complex Z,
            typename B = b_rows,
            typename A = a_rows>
  void split_expanded_block(std::size_t                 m,
                            std::size_t                 nb,
                            std::size_t                 kb,
                            std::size_t                 pc,
                            std::size_t                 jc,
                            const Z*                    a,
                            std::size_t                 lda,
                            const Z*                    b,
                            std::size_t                 ldb,
                            typename Z::value_type*     ar,
                            typename Z::value_type*     br,
                            typename Z::value_type*     cr)
  {
    using R = typename Z::value_type;
    parallel_for(m, kb, [&](std::size_t i) {
                   for (std::size_t p {0}; p < kb; ++p) {
                     const Z value {A::at(a, lda, i, pc + p)};
                     ar[(i * 2 * kb) + p]      = value.real();
                     ar[(i * 2 * kb) + kb + p] = -value.imag();
                   }
                 });
    parallel_for(kb, 2 * nb, [&](std::size_t p) {
                   R* const top {br + (p * 2 * nb)};
                   R* const bottom {br + ((kb + p) * 2 * nb)};
                   for (std::size_t j {0}; j < nb; ++j) {
                     const Z value {B::at(b, ldb, pc + p, jc + j)};
                     top[j]         = value.real();
                     top[nb + j]    = value.imag();
                     bottom[j]      = value.imag();
                     bottom[nb + j] = -value.real();
                   }
                 });
    gemm(m, 2 * nb, 2 * kb, ar, 2 * kb, br, 2 * nb, cr, 2 * nb);
  }
 
  template <split_complex Z,
            typename B = b_rows,
            typename A = a_rows>
  void gemm_split(std::size_t m,
                  std::size_t n,
                  std::size_t k,
                  const Z*    a,
                  std::size_t lda,
                  const Z*    b,
                  std::size_t ldb,
                  Z*          c,
                  std::size_t ldc)
  {
    using R = typename Z::value_type;
    const bool                  parts   {split_by_parts(m)};
    const std::size_t           widest  {std::min(n, split_columns)};
    const std::size_t           deepest {std::min(k, split_depth)};
    split_operands<R>&          kept    {split_storage<R>()};
    kept.a.resize(std::max(kept.a.size(), m * (parts ? 3 : 2) * deepest)); // grown only, written before they are read
    kept.b.resize(std::max(kept.b.size(), (parts ? 2 : 4) * deepest * widest));
    kept.c.resize(std::max(kept.c.size(), m * 2 * widest));
    R* const                    ar      {kept.a.data()};
    R* const                    br      {kept.b.data()};
    R* const                    cr      {kept.c.data()};
    for (std::size_t jc {0}; jc < n; jc += split_columns) {
      const std::size_t nb {std::min(split_columns, n - jc)};
      parallel_for(m, nb, [&](std::size_t i) {
                     for (std::size_t j {0}; j < nb; ++j) {
                       cr[(i * 2 * nb) + j]      = c[(i * ldc) + jc + j].real();
                       cr[(i * 2 * nb) + nb + j] = c[(i * ldc) + jc + j].imag();
                     }
                   });
      for (std::size_t pc {0}; pc < k; pc += split_depth) {
        const std::size_t kb {std::min(split_depth, k - pc)};
        if (parts) {
          split_parts_block<Z, B, A>(m, nb, kb, pc, jc, a, lda, b, ldb, ar, br, cr);
        }
        else {
          split_expanded_block<Z, B, A>(m, nb, kb, pc, jc, a, lda, b, ldb, ar, br, cr);
        }
      }
      parallel_for(m, nb, [&](std::size_t i) {
                     for (std::size_t j {0}; j < nb; ++j) {
                       c[(i * ldc) + jc + j] = Z {cr[(i * 2 * nb) + j], cr[(i * 2 * nb) + nb + j]};
                     }
                   });
    }
    for (buffer<R>* operand : {&kept.a, &kept.b, &kept.c}) {
      if (operand->size() > split_kept) {
        *operand = buffer<R> {};
      }
    }
  }
 
  inline constexpr std::size_t plain_loop_least_width {64};
 
  [[nodiscard]] inline bool plain_loop_pays(std::size_t m,
                                            std::size_t n,
                                            std::size_t k) noexcept
  {
    [[maybe_unused]] const gemm_path path {gemm_path_current()};
#ifdef NUMERICAL_NEON_COMPILED
    if (path == gemm_path::neon) {
      return k <= 16 || (m <= 64 && k <= 64);
    }
#endif
#ifdef NUMERICAL_X86_COMPILED
    if (path == gemm_path::avx2 || path == gemm_path::avx512) {
      return n >= plain_loop_least_width && ((m <= 16 && k <= 64) || (m <= 32 && k <= 32));
    }
#endif
    (void)m, (void)n, (void)k;
    return false;
  }
 
  template <typename T>
  [[nodiscard]] bool plain_loop_suits(std::size_t /*m*/,
                                      std::size_t /*n*/,
                                      std::size_t /*k*/) noexcept
  {
    return false;
  }
 
  template <>
  [[nodiscard]] inline bool plain_loop_suits<double>(std::size_t m,
                                                     std::size_t n,
                                                     std::size_t k) noexcept
  {
    return plain_loop_pays(m, n, k);
  }
 
  template <typename T>
  void gemm_large(std::size_t m,
                  std::size_t n,
                  std::size_t k,
                  const T*    a,
                  std::size_t lda,
                  const T*    b,
                  std::size_t ldb,
                  T*          c,
                  std::size_t ldc)
  {
    if (plain_loop_suits<T>(m, n, k)) { // the plain loop rounds as the tiles do: the same bits either way
      gemm_rows(m, n, k, a, lda, b, ldb, c, ldc);
      return;
    }
    if (!gemm_elsewhere(m, n, k, a, lda, b, ldb, c, ldc)) {
      gemm_packed<vector_tile_kernel<T>>(m, n, k, a, lda, b, ldb, c, ldc);
    }
  }
 
#ifdef NUMERICAL_X86_COMPILED
  template <split_complex Z,
            typename Double,
            typename Float>
  using x86_complex_kernel = std::conditional_t<std::is_same_v<Z, std::complex<double>>, Double, Float>;
#endif
 
  template <typename B,
            typename A,
            split_complex T>
  void gemm_large_complex(std::size_t m,
                          std::size_t n,
                          std::size_t k,
                          const T*    a,
                          std::size_t lda,
                          const T*    b,
                          std::size_t ldb,
                          T*          c,
                          std::size_t ldc)
  {
    const gemm_path path {gemm_path_current()};
    if (path == gemm_path::portable) {
      gemm_packed<vector_tile_kernel<T>, T, B, A>(m, n, k, a, lda, b, ldb, c, ldc);
      return;
    }
#ifdef NUMERICAL_X86_COMPILED
    if (path == gemm_path::avx512) { // a complex tile of its own on the x86-64 units
      gemm_packed<x86_complex_kernel<T, avx512_complex_tile_kernel, avx512_complex_tile_kernel_f32>, T, B, A>(m, n, k, a, lda, b, ldb, c, ldc);
      return;
    }
    if (path == gemm_path::avx2) {
      gemm_packed<x86_complex_kernel<T, avx2_complex_tile_kernel, avx2_complex_tile_kernel_f32>, T, B, A>(m, n, k, a, lda, b, ldb, c, ldc);
      return;
    }
#endif
    gemm_split<T, B, A>(m, n, k, a, lda, b, ldb, c, ldc);
  }
 
  template <split_complex T>
  void gemm_large(std::size_t m,
                  std::size_t n,
                  std::size_t k,
                  const T*    a,
                  std::size_t lda,
                  const T*    b,
                  std::size_t ldb,
                  T*          c,
                  std::size_t ldc)
  {
    gemm_large_complex<b_rows, a_rows>(m, n, k, a, lda, b, ldb, c, ldc);
  }
 
  template <typename T>
  void gemm(std::size_t m,
            std::size_t n,
            std::size_t k,
            const T*    a,
            std::size_t lda,
            const T*    b,
            std::size_t ldb,
            T*          c,
            std::size_t ldc)
  {
    if constexpr (split_complex<T>) {
      if (gemm_path_current() != gemm_path::portable) {
        if (m != 0 && n != 0 && k != 0) {
          gemm_large(m, n, k, a, lda, b, ldb, c, ldc);
        }
        return;
      }
    }
    if (k < gemm_min_depth || m * n < ((gemm_packed_threshold + k) - 1) / k) {
      gemm_rows(m, n, k, a, lda, b, ldb, c, ldc);
      return;
    }
    gemm_large(m, n, k, a, lda, b, ldb, c, ldc);
  }
 
  template <typename B,
            typename A>
  void gemm_large_layouts(std::size_t   m,
                          std::size_t   n,
                          std::size_t   k,
                          const double* a,
                          std::size_t   lda,
                          const double* b,
                          std::size_t   ldb,
                          double*       c,
                          std::size_t   ldc)
  {
    if (!gemm_double_elsewhere<B, A>(m, n, k, a, lda, b, ldb, c, ldc)) {
      gemm_packed<vector_tile_kernel<double>, double, B, A>(m, n, k, a, lda, b, ldb, c, ldc);
    }
  }
 
  template <typename B,
            typename A>
  void gemm_large_layouts(std::size_t  m,
                          std::size_t  n,
                          std::size_t  k,
                          const float* a,
                          std::size_t  lda,
                          const float* b,
                          std::size_t  ldb,
                          float*       c,
                          std::size_t  ldc)
  {
    if (!gemm_float_elsewhere<B, A>(m, n, k, a, lda, b, ldb, c, ldc)) {
      gemm_packed<vector_tile_kernel<float>, float, B, A>(m, n, k, a, lda, b, ldb, c, ldc);
    }
  }
 
  template <typename B,
            typename A,
            split_complex T>
  void gemm_large_layouts(std::size_t m,
                          std::size_t n,
                          std::size_t k,
                          const T*    a,
                          std::size_t lda,
                          const T*    b,
                          std::size_t ldb,
                          T*          c,
                          std::size_t ldc)
  {
    gemm_large_complex<B, A>(m, n, k, a, lda, b, ldb, c, ldc);
  }
 
  template <typename A,
            typename B,
            typename T>
  void gemm_layouts(std::size_t m,
                    std::size_t n,
                    std::size_t k,
                    const T*    a,
                    std::size_t lda,
                    const T*    b,
                    std::size_t ldb,
                    T*          c,
                    std::size_t ldc)
  {
    if constexpr (split_complex<T>) {
      if (gemm_path_current() != gemm_path::portable) {
        if (m != 0 && n != 0 && k != 0) {
          gemm_large_complex<B, A>(m, n, k, a, lda, b, ldb, c, ldc);
        }
        return;
      }
    }
    if (k < gemm_min_depth || m * n < ((gemm_packed_threshold + k) - 1) / k) {
      buffer<T>      a_rows_copy;
      buffer<T>      b_rows_copy;
      std::size_t    a_ld  {};
      std::size_t    b_ld  {};
      const T* const left  {A::rows_of(a, lda, m, k, a_rows_copy, a_ld)};
      const T* const right {B::rows_of(b, ldb, k, n, b_rows_copy, b_ld)};
      gemm_rows(m, n, k, left, a_ld, right, b_ld, c, ldc);
      return;
    }
    gemm_large_layouts<B, A>(m, n, k, a, lda, b, ldb, c, ldc);
  }
 
  template <typename T>
  void gemm_bt(std::size_t m,
               std::size_t n,
               std::size_t k,
               const T*    a,
               std::size_t lda,
               const T*    bt,
               std::size_t ldbt,
               T*          c,
               std::size_t ldc)
  {
    gemm_layouts<a_rows, b_transposed>(m, n, k, a, lda, bt, ldbt, c, ldc);
  }
 
  template <typename T>
  void gemm_at(std::size_t m,
               std::size_t n,
               std::size_t k,
               const T*    at,
               std::size_t ldat,
               const T*    b,
               std::size_t ldb,
               T*          c,
               std::size_t ldc)
  {
    gemm_layouts<a_transposed, b_rows>(m, n, k, at, ldat, b, ldb, c, ldc);
  }
 
  template <typename T>
  void gemm_atbt(std::size_t m,
                 std::size_t n,
                 std::size_t k,
                 const T*    at,
                 std::size_t ldat,
                 const T*    bt,
                 std::size_t ldbt,
                 T*          c,
                 std::size_t ldc)
  {
    gemm_layouts<a_transposed, b_transposed>(m, n, k, at, ldat, bt, ldbt, c, ldc);
  }
 
  enum class left_kernel : unsigned char
  {
    plain,  
    vector, 
    neon,   
    amx,    
    avx512, 
    avx2    
  };
 
  template <typename T>
  struct packed_left
  {
    left_kernel kernel;  
    std::size_t m;       
    std::size_t k;       
    const T*    a;       
    std::size_t lda;     
    const T*    packed;  
    std::size_t granule; 
  };
 
  template <typename T>
  [[nodiscard]] packed_left<T> read_in_place(std::size_t m,
                                             std::size_t k,
                                             const T*    a,
                                             std::size_t lda) noexcept
  {
    return packed_left<T> {.kernel = left_kernel::plain, .m = m, .k = k, .a = a, .lda = lda, .packed = nullptr, .granule = 1};
  }
 
  template <typename T>
  [[nodiscard]] packed_left<T> left_rows(const packed_left<T>& left,
                                         std::size_t           r0,
                                         std::size_t           rows) noexcept
  {
    packed_left<T> part {left};
    part.m      = rows;
    part.a      = left.a + (r0 * left.lda);
    part.packed = left.packed == nullptr ? nullptr : left.packed + (r0 * left.k);
    return part;
  }
 
  template <typename T>
  [[nodiscard]] left_kernel left_kernel_for(std::size_t /*m*/,
                                            std::size_t /*n*/,
                                            std::size_t /*k*/) noexcept
  {
    return left_kernel::plain;
  }
 
  template <>
  [[nodiscard]] inline left_kernel left_kernel_for<double>(std::size_t m,
                                                           std::size_t n,
                                                           std::size_t k) noexcept
  {
    if (k < gemm_min_depth || m * n < ((gemm_packed_threshold + k) - 1) / k || k > gemm_kc) {
      return left_kernel::plain;
    }
    const gemm_path path {gemm_path_current()};
#ifdef NUMERICAL_AMX_COMPILED
    if (path == gemm_path::amx && amx_suits(m, n, k)) {
      return left_kernel::amx;
    }
#endif
#ifdef NUMERICAL_NEON_COMPILED
    if (path == gemm_path::neon || path == gemm_path::amx) {
      return left_kernel::neon;
    }
#endif
#ifdef NUMERICAL_X86_COMPILED
    if (path == gemm_path::avx512) {
      return left_kernel::avx512;
    }
    if (path == gemm_path::avx2) {
      return left_kernel::avx2;
    }
#endif
    (void)path;
    return left_kernel::vector;
  }
 
  template <typename Kernel,
            typename T>
  [[nodiscard]] const T* pack_left_for(std::size_t m,
                                       std::size_t k,
                                       const T*    a,
                                       std::size_t lda,
                                       buffer<T>&  storage)
  {
    const std::size_t a_rows {((m + Kernel::mr) - 1) / Kernel::mr * Kernel::mr};
    storage.resize(std::max(storage.size(), (a_rows * k) + (packed_line / sizeof(T)) + 1));
    T* const          out    {line_aligned(storage.data())};
    const std::size_t blocks {((m + Kernel::block_rows) - 1) / Kernel::block_rows};
    parallel_for(blocks, Kernel::block_rows * k * gemm_kernel_speedup / Kernel::speedup, [&](std::size_t block) {
                   const std::size_t ic {block * Kernel::block_rows};
                   Kernel::pack_a(a + (ic * lda), lda, std::min(Kernel::block_rows, m - ic), k, out + (ic * k));
                 });
    return out;
  }
 
  template <typename T>
  [[nodiscard]] packed_left<T> pack_left(std::size_t m,
                                         std::size_t /*n*/,
                                         std::size_t k,
                                         const T*    a,
                                         std::size_t lda,
                                         buffer<T>&  /*storage*/)
  {
    return read_in_place(m, k, a, lda);
  }
 
  template <>
  [[nodiscard]] inline packed_left<double> pack_left<double>(std::size_t       m,
                                                             std::size_t       n,
                                                             std::size_t       k,
                                                             const double*     a,
                                                             std::size_t       lda,
                                                             buffer<double>&   storage)
  {
    packed_left<double> left {.kernel = left_kernel_for<double>(m, n, k), .m = m, .k = k, .a = a, .lda = lda, .packed = nullptr, .granule = 1};
#ifdef NUMERICAL_AMX_COMPILED
    if (left.kernel == left_kernel::amx) {
      left.packed  = pack_left_for<amx_tile_kernel>(m, k, a, lda, storage);
      left.granule = amx_tile_kernel::mr;
    }
#endif
#ifdef NUMERICAL_NEON_COMPILED
    if (left.kernel == left_kernel::neon) {
      left.packed  = pack_left_for<neon_tile_kernel<double>>(m, k, a, lda, storage);
      left.granule = neon_tile_kernel<double>::mr;
    }
#endif
#ifdef NUMERICAL_X86_COMPILED
    if (left.kernel == left_kernel::avx512) {
      left.packed  = pack_left_for<avx512_tile_kernel>(m, k, a, lda, storage);
      left.granule = avx512_tile_kernel::mr;
    }
    if (left.kernel == left_kernel::avx2) {
      left.packed  = pack_left_for<avx2_tile_kernel>(m, k, a, lda, storage);
      left.granule = avx2_tile_kernel::mr;
    }
#endif
    if (left.kernel == left_kernel::vector) {
      left.packed  = pack_left_for<vector_tile_kernel<double>>(m, k, a, lda, storage);
      left.granule = vector_tile_kernel<double>::mr;
    }
    return left;
  }
 
  template <typename T>
  void gemm_with(const packed_left<T>& left,
                 std::size_t           n,
                 const T*              b,
                 std::size_t           ldb,
                 T*                    c,
                 std::size_t           ldc)
  {
    gemm(left.m, n, left.k, left.a, left.lda, b, ldb, c, ldc);
  }
 
  template <>
  inline void gemm_with<double>(const packed_left<double>& left,
                                std::size_t                n,
                                const double*              b,
                                std::size_t                ldb,
                                double*                    c,
                                std::size_t                ldc)
  {
    if (left.kernel == left_kernel::plain || left_kernel_for<double>(left.m, n, left.k) != left.kernel) {
      gemm(left.m, n, left.k, left.a, left.lda, b, ldb, c, ldc);
      return;
    }
#ifdef NUMERICAL_AMX_COMPILED
    if (left.kernel == left_kernel::amx) {
      gemm_packed<amx_tile_kernel>(left.m, n, left.k, left.a, left.lda, b, ldb, c, ldc, left.packed);
      return;
    }
#endif
#ifdef NUMERICAL_NEON_COMPILED
    if (left.kernel == left_kernel::neon) {
      gemm_packed<neon_tile_kernel<double>>(left.m, n, left.k, left.a, left.lda, b, ldb, c, ldc, left.packed);
      return;
    }
#endif
#ifdef NUMERICAL_X86_COMPILED
    if (left.kernel == left_kernel::avx512) {
      gemm_packed<avx512_tile_kernel>(left.m, n, left.k, left.a, left.lda, b, ldb, c, ldc, left.packed);
      return;
    }
    if (left.kernel == left_kernel::avx2) {
      gemm_packed<avx2_tile_kernel>(left.m, n, left.k, left.a, left.lda, b, ldb, c, ldc, left.packed);
      return;
    }
#endif
    gemm_packed<vector_tile_kernel<double>>(left.m, n, left.k, left.a, left.lda, b, ldb, c, ldc, left.packed);
  }
 
  inline constexpr std::size_t dot_chunk {std::size_t {1} << 14U};
 
  template <typename T>
  [[nodiscard]] T chunked_dot(const T*    a,
                              const T*    b,
                              std::size_t count)
  {
    const std::size_t chunks {(count + dot_chunk - 1) / dot_chunk};
    if (chunks <= 1) {
      return dot_product(a, b, count);
    }
    buffer<T> partial(chunks);
    parallel_for(chunks, dot_chunk, [&](std::size_t c) {
                   const std::size_t first {c * dot_chunk};
                   partial[c] = dot_product(a + first, b + first, std::min(dot_chunk, count - first));
                 });
    T total {};
    for (const T& sum : partial) {
      total += sum;
    }
    return total;
  }
 
  template <typename T>
  inline constexpr std::size_t multiply_add_weight {1};
 
  template <>
  inline constexpr std::size_t multiply_add_weight<std::complex<double>> {4};
 
  template <>
  inline constexpr std::size_t multiply_add_weight<std::complex<float>> {4};
 
  inline constexpr std::size_t row_group {64};
 
#ifdef NUMERICAL_AMX_COMPILED
  inline constexpr std::size_t amx_dot_least {std::size_t {1} << 15U};
 
  inline constexpr std::size_t amx_dot_most {std::size_t {1} << 18U};
#endif
 
  template <typename T>
  void matrix_times_vector(const T*    a,
                           std::size_t rows,
                           std::size_t cols,
                           const T*    x,
                           T*          y)
  {
#ifdef NUMERICAL_AMX_COMPILED
    if constexpr (std::is_same_v<T, double>) {
      const std::size_t elements {rows * cols};
      if (cols <= dot_chunk && elements >= amx_dot_least && elements < amx_dot_most && gemm_path_current() == gemm_path::amx) {
        const amx_session on {}; // on the caller alone
        amx_dot_rows_f64(a, cols, rows, x, cols, y);
        return;
      }
    }
#endif
    if (cols > dot_chunk) {   // rows longer than a chunk: chunked_dot's partial sums
      parallel_for(rows, cols, [&](std::size_t r) { y[r] = chunked_dot(a + (r * cols), x, cols); });
      return;
    }
    // A complex multiply-add weighs four real ones. Groups of fewer rows where the matrix has too few for each thread
    // to take some (189 rows of 64 were three groups: six AVX2 cores' L3 read by one, a complex matrix of 189 × 206
    // times a vector 24 µs); each row its own dot product whatever the grouping, so the same bits.
    constexpr std::size_t weight {multiply_add_weight<T>};
    const std::size_t     group  {std::clamp<std::size_t>(rows / (4 * max_threads()), 8, row_group)};
    parallel_for((rows + group - 1) / group, group * cols * weight, [&](std::size_t g) {
                   const std::size_t first {g * group};
                   row_dots(a + (first * cols), cols, std::min(group, rows - first), x, cols, y + first);
                 });
  }
} // namespace numerical::detail
 
namespace numerical {
 
  inline gemm_path limit_gemm_path(gemm_path most) noexcept
  {
    detail::gemm_path_cap().store(most, std::memory_order_relaxed);
    return detail::best_gemm_path(most); // probes only the paths at or below the cap
  }
 
  [[nodiscard]] inline gemm_path gemm_path_in_use() noexcept
  {
    return detail::gemm_path_current();
  }
} // namespace numerical
 
#endif // NUMERICAL_GEMM_HPP