// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-FileCopyrightText: The Eigen Authors
// SPDX-License-Identifier: MPL-2.0

// SME GEMM kernel tests.
// Requires a compiler targeting SME2, e.g. -march=armv9.2-a+sme2, which selects
// the backend on its own; -DEIGEN_ARM64_USE_SME then makes a toolchain that
// cannot provide it fail the build rather than fall back to NEON.
// double and complex<double> additionally need FEAT_SME_F64F64 (+sme-f64f64, or
// a -mcpu that implies it); without it EIGEN_VECTORIZE_SME_F64F64 is undefined
// and they keep the generic kernel, so their subtests pack their cases through
// that path instead.

#ifdef EIGEN_TEST_PART_6
#undef EIGEN_SME_NEON_SMALL_PRODUCT
#define EIGEN_SME_NEON_SMALL_PRODUCT 65536
#undef EIGEN_DEFAULT_DENSE_INDEX_TYPE
#define EIGEN_DEFAULT_DENSE_INDEX_TYPE int
#endif

#ifdef EIGEN_TEST_PART_7
#ifndef EIGEN_SME_NO_NEON_SMALL_BLOCKS
#define EIGEN_SME_NO_NEON_SMALL_BLOCKS
#endif
#endif

#ifdef EIGEN_TEST_PART_8
#ifndef EIGEN_SME_FORCE_NEON_SMALL_BLOCKS
#define EIGEN_SME_FORCE_NEON_SMALL_BLOCKS
#endif
#endif

#include "product.h"

#include <array>
#include <cfenv>

// Without the right -march flags, __ARM_FEATURE_SME is undefined and
// EIGEN_VECTORIZE_SME never fires - the test would silently compile
// against the NEON GEBP kernel and pass, making this a useless no-op.
// Fail the build instead.
#if !defined(EIGEN_VECTORIZE_SME)
#error \
    "product_sme requires the SME backend.  Build with -march=armv9.2-a+sme2 " \
    "(see -DEIGEN_TEST_SME=ON in test/CMakeLists.txt for the typical CMake " \
    "invocation)."
#endif

// The backend is chosen from the compiler's target features; EIGEN_ARM64_USE_SME,
// which this test's CMake option passes, only turns a toolchain that cannot
// provide it into a build error. Assert that here, so a regression making
// selection depend on the opt-in again is caught rather than masked by it.
#if !defined(EIGEN_ARM64_SME_SELECTED)
#error "the SME backend was forced by EIGEN_ARM64_USE_SME rather than selected from the compiler's target features."
#endif

template <typename Scalar>
using SmeColMajorMat = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
template <typename Scalar>
using SmeRowMajorMat = Matrix<Scalar, Dynamic, Dynamic, RowMajor>;
template <typename Scalar>
using SmeVector = Matrix<Scalar, Dynamic, 1>;
template <typename Scalar>
using SmeColMajorStridedMat = Map<SmeColMajorMat<Scalar>, 0, Stride<Dynamic, Dynamic>>;
template <typename Scalar>
using SmeRowMajorStridedMat = Map<SmeRowMajorMat<Scalar>, 0, Stride<Dynamic, Dynamic>>;

// The logical micro-kernel block for Scalar (kSmeMr/kSmeNr for float and so
// on).  Sizes below are expressed in terms of it so each scalar sweeps its own
// block boundaries; sme_tile() is the ZA tile side at the SVL=512 design point,
// which is where the intra-block splits fall.  Complex scalars pair two tiles
// per accumulator, so their block is not square.
template <typename Scalar>
static constexpr int sme_mr() {
  return internal::sme_block<Scalar>::mr;
}
template <typename Scalar>
static constexpr int sme_nr() {
  return internal::sme_block<Scalar>::nr;
}
template <typename Scalar>
static constexpr int sme_tile() {
  return internal::sme_block<Scalar>::mr / internal::sme_block<Scalar>::kGridRows;
}

// svwhilelt_b* is overloaded only on the four fixed-width types, so forwarding
// the caller's own type would leave Index ambiguous wherever it is `long` while
// int64_t is `long long`, and deducing one type from both arguments would reject
// the literal-plus-Index bounds the packers pass.  Name every spelling Index
// takes, mix the two argument types, and check the predicate so the traits'
// widening cannot change which lanes are active.
// EIGEN_DONT_INLINE here and below: GCC (14.1 through trunk) inlines a __arm_locally_streaming function into a
// non-streaming caller without its mode switch, running the body as non-streaming SVE at that vector length.
template <typename Scalar, typename Begin, typename End = Begin>
EIGEN_DONT_INLINE __arm_locally_streaming static bool sme_whilelt_covers_first_lane_only() {
  using Traits = internal::sme_packet_traits<Scalar>;
  const svbool_t pg = Traits::whilelt(Begin(0), End(1));
  return svptest_first(Traits::ptrue(), pg) && svcntp_b8(Traits::ptrue(), pg) == 1;
}

template <typename Scalar>
static void test_whilelt_operand_types() {
  VERIFY((sme_whilelt_covers_first_lane_only<Scalar, int>()));
  VERIFY((sme_whilelt_covers_first_lane_only<Scalar, long>()));
  VERIFY((sme_whilelt_covers_first_lane_only<Scalar, long long>()));
  VERIFY((sme_whilelt_covers_first_lane_only<Scalar, Index>()));
  VERIFY((sme_whilelt_covers_first_lane_only<Scalar, int, Index>()));
}

// Column j of out (n rows, one streaming vector) receives the j-th operation's result.
template <typename Scalar>
EIGEN_DONT_INLINE __arm_locally_streaming static Scalar sme_packet_ops(Scalar* out, const Scalar* a, const Scalar* b,
                                                                       const Scalar* c, Index active) {
  using Traits = internal::sme_packet_traits<Scalar>;
  using Vec = typename Traits::type;
  const svbool_t all = Traits::ptrue();
  const svbool_t pg = Traits::whilelt(0, active);
  const Index n = Traits::size();
  const Vec va = internal::ploadu(all, a), vb = internal::ploadu(all, b), vc = internal::ploadu(all, c);
  internal::pstoreu(all, out, internal::pmadd(all, va, vb, vc));
  internal::pstoreu(all, out + n, internal::pnmadd(all, va, vb, vc));
  internal::pstoreu(all, out + 2 * n, internal::pmadd_m(pg, va, vb, vc));
  internal::pstoreu(all, out + 3 * n, internal::padd(all, va, vb));
  internal::pstoreu(all, out + 4 * n, internal::pmul(all, va, vb));
  internal::pstoreu(all, out + 5 * n, internal::pnegate(all, va));
  internal::pstoreu(pg, out + 6 * n, internal::pset1<Vec>(Scalar(7)));
  return internal::predux(pg, va);
}

// Operand order and predication of the packet layer, on small integers so every result is exact.
// a_i = i + 2, b_i = -(i + 3), c_i = 2i + 5 are pairwise distinct and never +-1, so swapping two
// operands of pmadd or pnmadd changes every lane: ab + c - (ac + b) = (b - c)(a - 1), and so on.
template <typename Scalar>
static void test_packet_ops() {
  using Vec = Matrix<Scalar, Dynamic, 1>;
  using Mat = Matrix<Scalar, Dynamic, Dynamic>;
  const Index n = internal::sme_packet_traits<Scalar>::size();
  const Vec a = Vec::LinSpaced(n, Scalar(2), Scalar(n + 1));
  const Vec b = Vec::LinSpaced(n, Scalar(-3), Scalar(-(n + 2)));
  const Vec c = Vec::LinSpaced(n, Scalar(5), Scalar(2 * n + 3));
  for (Index active = 0; active <= n; ++active) {
    Mat out = Mat::Constant(n, 7, Scalar(-1));
    const Scalar sum = sme_packet_ops(out.data(), a.data(), b.data(), c.data(), active);
    Mat ref(n, 7);
    ref.col(0) = a.cwiseProduct(b) + c;
    ref.col(1) = c - a.cwiseProduct(b);
    ref.col(2) = c;
    ref.col(2).head(active) = ref.col(0).head(active);
    ref.col(3) = a + b;
    ref.col(4) = a.cwiseProduct(b);
    ref.col(5) = -a;
    ref.col(6).setConstant(Scalar(-1));
    ref.col(6).head(active).setConstant(Scalar(7));
    VERIFY_IS_EQUAL(out, ref);
    VERIFY_IS_EQUAL(sum, a.head(active).sum());
  }
}

// Write one element into a packed panel of width w, in the layout the SME
// kernel reads: a real scalar lands at dst[k*w + r], a complex one splits into
// the depth step's real and imaginary halves.
template <typename Scalar>
static void set_packed(Scalar* panel, Index w, Index k, Index r, const Scalar& v) {
  panel[k * w + r] = v;
}
template <typename RealScalar>
static void set_packed(std::complex<RealScalar>* panel, Index w, Index k, Index r, const std::complex<RealScalar>& v) {
  RealScalar* p = reinterpret_cast<RealScalar*>(panel + k * w);
  p[r] = numext::real(v);
  p[w + r] = numext::imag(v);
}

template <typename InputMat, typename ResultMat, typename ResultMap>
static void verify_strided_result(int n, ResultMat& storage, const Stride<Dynamic, Dynamic>& stride) {
  InputMat A = InputMat::Random(n, n);
  InputMat B = InputMat::Random(n, n);
  ResultMap C(storage.data(), n, n, stride);
  C = ResultMat::Random(n, n);
  ResultMat c_before = C.eval();

  C.noalias() += A * B;

  ResultMat ref = c_before + (A.lazyProduct(B)).eval();
  ResultMat got = C;
  VERIFY_IS_APPROX(got, ref);
}

template <typename Scalar, typename InputMat>
static void test_general_strided_result(int n) {
  // General-stride C path: InputMat selects the source packers, while both C
  // strides are non-unit so sme_store_za_tile uses scalar scatter.
  SmeColMajorMat<Scalar> storage = SmeColMajorMat<Scalar>::Zero(2 * n, n);
  verify_strided_result<InputMat, SmeColMajorMat<Scalar>, SmeColMajorStridedMat<Scalar>>(
      n, storage, Stride<Dynamic, Dynamic>(/*outer=*/2 * n, /*inner=*/2));

  // Padding rows skipped by the strided Map should not be touched.
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < n; ++j) {
      VERIFY(storage(2 * i + 1, j) == Scalar(0));
    }
  }
}

template <typename Scalar>
static void test_rowmajor_strided_result(int n) {
  // RowMajor C path: inner stride is one, with padded columns after the Map.
  SmeRowMajorMat<Scalar> storage = SmeRowMajorMat<Scalar>::Zero(n, 2 * n);
  verify_strided_result<SmeRowMajorMat<Scalar>, SmeRowMajorMat<Scalar>, SmeRowMajorStridedMat<Scalar>>(
      n, storage, Stride<Dynamic, Dynamic>(/*outer=*/2 * n, /*inner=*/1));

  // Padding columns skipped by the strided Map should not be touched.
  for (int i = 0; i < n; ++i) {
    for (int j = n; j < 2 * n; ++j) {
      VERIFY(storage(i, j) == Scalar(0));
    }
  }
}

// A non-trivial alpha: complex where the scalar is, so the store's rotation of
// the accumulator is exercised rather than a plain rescale.
template <typename Scalar>
struct nontrivial_alpha_impl {
  static Scalar run() { return Scalar(1.375); }
};
template <typename RealScalar>
struct nontrivial_alpha_impl<std::complex<RealScalar>> {
  static std::complex<RealScalar> run() { return std::complex<RealScalar>(RealScalar(1.375), RealScalar(-0.625)); }
};

// Conjugated operands and a non-real alpha.
//
// blas_traits folds .adjoint()/.conjugate() and a scalar factor into the
// gebp_kernel's ConjugateLhs/ConjugateRhs and alpha rather than materializing a
// temporary, so these expressions are the only thing that reaches the complex
// kernel's FMOPA/FMOPS sign choices and its alpha-scaled store. The references
// materialize the conjugation first and multiply coefficient-wise, so they do
// not share a code path with what they check.
template <typename Scalar, typename LhsMat, typename RhsMat, typename ResMat>
static void verify_conjugated_products(int n) {
  const LhsMat A = LhsMat::Random(n, n);
  const RhsMat B = RhsMat::Random(n, n);
  const SmeColMajorMat<Scalar> Ap = A, Bp = B;
  const SmeColMajorMat<Scalar> Aa = Ap.adjoint().eval(), Ac = Ap.conjugate().eval();
  const SmeColMajorMat<Scalar> Ba = Bp.adjoint().eval(), Bc = Bp.conjugate().eval();
  const Scalar alpha = nontrivial_alpha_impl<Scalar>::run();

  ResMat C(n, n);
  C.setZero();
  C.noalias() += A.adjoint() * B;
  VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), SmeColMajorMat<Scalar>(Aa.lazyProduct(Bp)));
  C.setZero();
  C.noalias() += A * B.adjoint();
  VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), SmeColMajorMat<Scalar>(Ap.lazyProduct(Ba)));
  C.setZero();
  C.noalias() += A.adjoint() * B.adjoint();
  VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), SmeColMajorMat<Scalar>(Aa.lazyProduct(Ba)));
  C.setZero();
  C.noalias() += A.conjugate() * B;
  VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), SmeColMajorMat<Scalar>(Ac.lazyProduct(Bp)));
  C.setZero();
  C.noalias() += A * B.conjugate();
  VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), SmeColMajorMat<Scalar>(Ap.lazyProduct(Bc)));
  C.setZero();
  C.noalias() += A.conjugate() * B.conjugate();
  VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), SmeColMajorMat<Scalar>(Ac.lazyProduct(Bc)));
  C.setZero();
  C.noalias() += alpha * (A.adjoint() * B);
  VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), SmeColMajorMat<Scalar>(alpha * Aa.lazyProduct(Bp)));

  // A deep rectangular product reaches GEMM even for half-tile output widths.
  const LhsMat deepA = LhsMat::Random(n, 65);
  const RhsMat deepB = RhsMat::Random(65, n);
  const ResMat initial = ResMat::Random(n, n);
  const SmeColMajorMat<Scalar> reference = deepA.lazyProduct(deepB);
  C = initial;
  C.noalias() += deepA * deepB;
  VERIFY(verifyProductUpdate(C, initial + reference, Scalar(1), deepA, deepB, initial));
  C = initial;
  C.noalias() += alpha * (deepA * deepB);
  VERIFY(verifyProductUpdate(C, initial + alpha * reference, alpha, deepA, deepB, initial));
  // A real alpha takes the complex store's fused path: -1 from -=, and a
  // non-unit real scale.
  const Scalar real_alpha(typename NumTraits<Scalar>::Real(1.375));
  C = initial;
  C.noalias() -= deepA * deepB;
  VERIFY(verifyProductUpdate(C, initial - reference, Scalar(-1), deepA, deepB, initial));
  C = initial;
  C.noalias() += real_alpha * (deepA * deepB);
  VERIFY(verifyProductUpdate(C, initial + real_alpha * reference, real_alpha, deepA, deepB, initial));

#ifdef EIGEN_VECTORIZE_SME_F64F64
  // RowMajor GEMM normally transposes the computation. Call the packed kernel
  // too, so RowMajor C actually reaches its horizontal ZA slice store.
  if (n <= sme_tile<Scalar>()) {
    std::vector<Scalar> packedA(n * 65), packedB(n * 65);
    for (Index k = 0; k < 65; ++k) {
      for (Index i = 0; i < n; ++i) {
        set_packed(packedA.data(), Index(n), k, i, deepA(i, k));
        set_packed(packedB.data(), Index(n), k, i, deepB(k, i));
      }
    }
    for (Scalar scale : {Scalar(1), alpha, Scalar(-1), real_alpha}) {
      C = initial;
      internal::sme_gebp_dispatch<Scalar, false, false>::run(
          C.data(), Index(ResMat::IsRowMajor ? n : 1), Index(ResMat::IsRowMajor ? 1 : n), packedA.data(),
          packedB.data(), Index(n), Index(65), Index(n), scale, Index(65), Index(65), Index(0), Index(0));
      VERIFY(verifyProductUpdate(C, initial + scale * reference, scale, deepA, deepB, initial));
    }
  }
#endif
}

// Exercise the kc split path just above the SME blocking heuristic's depth cap
// (sme_max_kc in GeneralBlockPanelKernel.h, scaled by the scalar width).
template <typename Scalar>
static void test_deep_k_split() {
  constexpr int rows = 64;
  const int depth = int(2 * (Index(EIGEN_SME_MAX_KC) * Index(sizeof(float)) / Index(sizeof(Scalar)))) + 2;
  constexpr int cols = 64;
  SmeColMajorMat<Scalar> A = SmeColMajorMat<Scalar>::Random(rows, depth);
  SmeColMajorMat<Scalar> B = SmeColMajorMat<Scalar>::Random(depth, cols);
  SmeColMajorMat<Scalar> C = SmeColMajorMat<Scalar>::Random(rows, cols);
  SmeColMajorMat<Scalar> c_before = C;

  C.noalias() += A * B;

  VERIFY_IS_APPROX(C, c_before + (A.lazyProduct(B)).eval());
}

// LHS read in place: strides on either side of the 4 KB-multiple rule (packed versus in place), an odd stride
// and base offset, a sub-block at a row offset, and a depth tail, each checked against lazyProduct.
template <typename Scalar>
static void test_direct_lhs() {
  const Index kPage = 4096 / Index(sizeof(Scalar));
  const Index rows = 97, depth = 101, cols = 70;
  const Index strides[] = {rows + 32 + 1, kPage, kPage + 16, 2 * kPage + 16};
  for (Index lda : strides) {
    for (Index offset = 0; offset <= 1; ++offset) {
      SmeVector<Scalar> storage = SmeVector<Scalar>::Zero(lda * depth + 1);
      Map<SmeColMajorMat<Scalar>, 0, OuterStride<>> A(storage.data() + offset, rows + 32, depth, OuterStride<>(lda));
      A.setRandom();
      const SmeColMajorMat<Scalar> B = SmeColMajorMat<Scalar>::Random(depth, cols);
      const SmeColMajorMat<Scalar> C0 = SmeColMajorMat<Scalar>::Random(rows + 32, cols);
      SmeColMajorMat<Scalar> C = C0;
      C.noalias() += A * B;
      VERIFY_IS_APPROX(C, C0 + A.lazyProduct(B));
      SmeColMajorMat<Scalar> D = C0.middleRows(3, rows);
      D.noalias() += A.middleRows(3, rows) * B;
      VERIFY_IS_APPROX(D, C0.middleRows(3, rows) + A.middleRows(3, rows).lazyProduct(B));
    }
  }
#ifdef EIGEN_VECTORIZE_SME
  const int units = nbSmeUnits();
  VERIFY(units >= 0);
  setNbSmeUnits(1);
  VERIFY_IS_EQUAL(nbSmeUnits(), 1);
  setNbSmeUnits(units);
  VERIFY_IS_EQUAL(nbSmeUnits(), units);
#endif
  // An LHS packed once for all column blocks, where the first block reads it in place and a narrower last one, whose
  // panel spans more than EIGEN_SME_DIRECT_LHS_MAX_SPAN_BYTES, packs.
  {
    const Index big = Index(EIGEN_SME_DIRECT_LHS_MAX_SPAN_BYTES) / (2048 * Index(sizeof(Scalar))) + 2;
    const SmeColMajorMat<Scalar> Abig = SmeColMajorMat<Scalar>::Random(big, 2048);
    const auto A = Abig.topRows(256);
    const SmeColMajorMat<Scalar> B = SmeColMajorMat<Scalar>::Random(2048, 532);
    SmeColMajorMat<Scalar> C(256, 532);
    C.noalias() = A * B;
    VERIFY_IS_APPROX(C, A.lazyProduct(B));
  }
}

// Small K against large M and N (NEON depths and just past them), and a RHS of one panel read in place: up to
// half a panel wide always, up to a full panel within the span limit, with 4 KB-multiple strides.
template <typename Scalar>
static void test_small_k_and_single_panel_rhs() {
  const Index MR = sme_mr<Scalar>(), NR = sme_nr<Scalar>();
  for (Index k : {1, 2, 3, 5, 8, 17, 24, 25, 40}) {
    for (Index mn : {Index(130), 4 * MR + 3}) {
      const SmeColMajorMat<Scalar> A = SmeColMajorMat<Scalar>::Random(mn, k), B = SmeColMajorMat<Scalar>::Random(k, mn);
      SmeColMajorMat<Scalar> C = SmeColMajorMat<Scalar>::Random(mn, mn);
      const SmeColMajorMat<Scalar> C0 = C;
      C.noalias() += A * B;
      VERIFY_IS_APPROX(C, C0 + A.lazyProduct(B));
    }
  }
  const Index page = 4096 / Index(sizeof(Scalar));
  for (Index n : {Index(1), NR / 2 - 1, NR / 2, NR / 2 + 1, NR}) {
    for (Index lda : {5 * MR + 3, page, 2 * page + 1}) {
      for (Index depth : {Index(40), Index(301)}) {
        const Index rows = numext::mini(lda, 5 * MR + 3);
        SmeVector<Scalar> storage = SmeVector<Scalar>::Zero(lda * depth);
        Map<SmeColMajorMat<Scalar>, 0, OuterStride<>> A(storage.data(), rows, depth, OuterStride<>(lda));
        A.setRandom();
        const SmeColMajorMat<Scalar> B = SmeColMajorMat<Scalar>::Random(depth, n);
        SmeColMajorMat<Scalar> C = SmeColMajorMat<Scalar>::Random(rows, n);
        const SmeColMajorMat<Scalar> C0 = C;
        C.noalias() += A * B;
        VERIFY_IS_APPROX(C, C0 + A.lazyProduct(B));
      }
    }
  }
  // A single-panel RHS whose depth block spans more than EIGEN_SME_DIRECT_LHS_MAX_SPAN_BYTES packs the LHS.
  const Index kc = Index(EIGEN_SME_MAX_KC) * Index(sizeof(float)) / Index(sizeof(Scalar));
  const Index wide = page * Index(std::size_t(EIGEN_SME_DIRECT_LHS_MAX_SPAN_BYTES) / (std::size_t(kc) * 4096) + 1);
  SmeVector<Scalar> storage = SmeVector<Scalar>::Zero(wide * kc);
  Map<SmeColMajorMat<Scalar>, 0, OuterStride<>> A(storage.data(), 2 * MR + 5, kc, OuterStride<>(wide));
  A.setRandom();
  const SmeColMajorMat<Scalar> B = SmeColMajorMat<Scalar>::Random(kc, NR);
  SmeColMajorMat<Scalar> C = SmeColMajorMat<Scalar>::Zero(2 * MR + 5, NR);
  C.noalias() += A * B;
  VERIFY_IS_APPROX(C, A.lazyProduct(B));
}

// Disjoint parts per SME unit (2 and 3 units) against a serial product. product_threaded covers the thread pool;
// this reaches the OpenMP branch in an EIGEN_TEST_OPENMP build and runs serially otherwise.
template <typename Scalar>
static void test_disjoint_parts() {
  using Col = SmeColMajorMat<Scalar>;
  using Row = SmeRowMajorMat<Scalar>;
  const int saved = nbSmeUnits();
  const int shapes[][3] = {{301, 259, 131}, {200, 40, 1024}, {20, 700, 300}};
  for (int units : {2, 3}) {
    setNbSmeUnits(units);
    for (const auto& sh : shapes) {
      const Col A = Col::Random(sh[0], sh[2]), B = Col::Random(sh[2], sh[1]);
      const Col ref = A.lazyProduct(B);
      Col C(sh[0], sh[1]);
      C.noalias() = A * B;
      VERIFY_IS_APPROX(C, ref);
      Row R(sh[0], sh[1]);
      R.noalias() = A * B;
      VERIFY_IS_APPROX(Col(R), ref);
      Col big = Col::Random(sh[0] + 5, sh[1] + 3);
      Col expected = big;
      expected.block(2, 1, sh[0], sh[1]) += Scalar(1.5) * ref;
      big.block(2, 1, sh[0], sh[1]).noalias() += Scalar(1.5) * A * B;
      VERIFY_IS_APPROX(big, expected);
    }
  }
  setNbSmeUnits(saved);
}

// Streaming-mode entry and exit set every FP exception flag: products on each SME path, with exact results, must leave
// the caller's flags as they were.
template <typename Scalar, typename Product>
static void verify_fp_flags_kept(Product product) {
#if defined(EIGEN_HAS_OPENMP) || defined(EIGEN_GEMM_THREADPOOL)
  // A threaded build's parallelize_gemm sizes the thread count in floating point, which can raise FE_INEXACT.
  const int checked = FE_ALL_EXCEPT & ~FE_INEXACT;
#else
  const int checked = FE_ALL_EXCEPT;
#endif
  std::feclearexcept(FE_ALL_EXCEPT);
  product();
  VERIFY_IS_EQUAL(std::fetestexcept(checked), 0);
  std::feraiseexcept(FE_INVALID);
  product();
  VERIFY_IS_EQUAL(std::fetestexcept(checked), FE_INVALID);
  std::feclearexcept(FE_ALL_EXCEPT);
}

template <typename Scalar>
static void test_fp_flags() {
  using Col = SmeColMajorMat<Scalar>;
  using Row = SmeRowMajorMat<Scalar>;
  const Index MR = sme_mr<Scalar>(), NR = sme_nr<Scalar>();
  const Index m = 4 * MR + 3, n = 3 * NR + 5, k = 131;
  // Small integers keep every product exact, so no flag is legitimately raised.
  const Col A = Col::Random(m, k).unaryExpr([](const Scalar& x) { return Scalar(numext::round(numext::real(x) * 4)); });
  const Col B = Col::Random(k, n).unaryExpr([](const Scalar& x) { return Scalar(numext::round(numext::real(x) * 4)); });
  const Col S = Col::Identity(m, m) * Scalar(2) + Col::Ones(m, m), O = Col::Ones(m, n);
  // I plus a subdiagonal of ones: the solution of L x = 1 is 1, 0, 1, 0, ..., exact.
  Col L = Col::Identity(m, m);
  L.diagonal(-1).setOnes();
  const Col Tall = Col::Ones(1024, k), Thin = Col::Ones(k, NR / 2);
  Col C(m, n), D(1024, NR / 2);
  Row R(m, n);
  verify_fp_flags_kept<Scalar>([&] { C.noalias() = A * B; });
  verify_fp_flags_kept<Scalar>([&] { C.noalias() += A.conjugate() * B; });
  verify_fp_flags_kept<Scalar>([&] { R.noalias() = A * B; });
  verify_fp_flags_kept<Scalar>([&] { R.noalias() = Row(A) * Row(B); });
  verify_fp_flags_kept<Scalar>([&] { D.noalias() = Tall * Thin; });
  verify_fp_flags_kept<Scalar>([&] { C.noalias() = S.template selfadjointView<Lower>() * O; });
  verify_fp_flags_kept<Scalar>([&] { C.noalias() = S.template triangularView<Upper>() * O; });
  verify_fp_flags_kept<Scalar>([&] {
    Col X = O;
    L.template triangularView<Lower>().solveInPlace(X);
  });
}

// Tiny results (the NEON kernel and both sides of its routing), every storage order of A, B and C.
template <typename Scalar, int AO, int BO, int CO>
static void test_tiny_results_order() {
  using MA = Matrix<Scalar, Dynamic, Dynamic, AO>;
  using MB = Matrix<Scalar, Dynamic, Dynamic, BO>;
  using MC = Matrix<Scalar, Dynamic, Dynamic, CO>;
  const Index ps = Index(16 / sizeof(Scalar));
  for (Index m : {Index(1), Index(2), Index(3), ps, ps + 1, 2 * ps, 2 * ps + 1})
    for (Index n : {1, 2, 3, 5, 8, 9})
      for (Index k : {1, 5, 6, 7, 15, 16, 17, 64, 2049}) {
        const MA A = MA::Random(m, k);
        const MB B = MB::Random(k, n);
        const MC ref = A.lazyProduct(B);
        MC C = MC::Random(m, n);
        const MC C0 = C;
        C.noalias() = A * B;
        VERIFY(verifyProduct(C, ref, A, B));
        C = C0;
        C.noalias() += A * B;
        VERIFY(verifyProductUpdate(C, C0 + ref, Scalar(1), A, B, C0));
        C = C0;
        C.noalias() -= Scalar(2) * A * B;
        VERIFY(verifyProductUpdate(C, C0 - Scalar(2) * ref, Scalar(-2), A, B, C0));
        const MA At = A.transpose();
        C.noalias() = At.transpose() * B;
        VERIFY(verifyProduct(C, ref, A, B));
      }
}

template <typename Scalar>
static void test_tiny_results() {
  test_tiny_results_order<Scalar, ColMajor, ColMajor, ColMajor>();
  test_tiny_results_order<Scalar, RowMajor, ColMajor, ColMajor>();
  test_tiny_results_order<Scalar, ColMajor, RowMajor, ColMajor>();
  test_tiny_results_order<Scalar, RowMajor, RowMajor, ColMajor>();
  test_tiny_results_order<Scalar, ColMajor, ColMajor, RowMajor>();
  test_tiny_results_order<Scalar, RowMajor, ColMajor, RowMajor>();
  test_tiny_results_order<Scalar, ColMajor, RowMajor, RowMajor>();
  test_tiny_results_order<Scalar, RowMajor, RowMajor, RowMajor>();
  // Sub-blocks of larger operands, a result with an inner stride, and a Hankel view (outer stride 1 < rows): the
  // kernel reads only the rows and columns of the result.
  for (Index m : {Index(3), Index(16 / sizeof(Scalar)) + 1})
    for (Index n : {3, 7})
      for (Index k : {16, 19}) {
        const SmeColMajorMat<Scalar> Abig = SmeColMajorMat<Scalar>::Random(m + 5, k + 2);
        const SmeRowMajorMat<Scalar> Bbig = SmeRowMajorMat<Scalar>::Random(k + 3, n + 4);
        const auto A = Abig.block(2, 1, m, k);
        const auto B = Bbig.block(1, 2, k, n);
        SmeColMajorMat<Scalar> Cbig = SmeColMajorMat<Scalar>::Random(2 * m, n);
        const SmeColMajorMat<Scalar> Cbig0 = Cbig;
        Map<SmeColMajorMat<Scalar>, 0, Stride<Dynamic, 2>> C(Cbig.data(), m, n, Stride<Dynamic, 2>(2 * m, 2));
        C.noalias() += Scalar(3) * A * B;
        SmeColMajorMat<Scalar> ref = Cbig0;
        Map<SmeColMajorMat<Scalar>, 0, Stride<Dynamic, 2>>(ref.data(), m, n, Stride<Dynamic, 2>(2 * m, 2)) +=
            Scalar(3) * A.lazyProduct(B);
        VERIFY_IS_APPROX(Cbig, ref);
        std::vector<Scalar> h(std::size_t(k - 1 + m)), g(std::size_t(k - 1 + n));
        for (auto& x : h) x = internal::random<Scalar>();
        for (auto& x : g) x = internal::random<Scalar>();
        Map<const SmeColMajorMat<Scalar>, 0, OuterStride<>> H(h.data(), m, k, OuterStride<>(1));
        Map<const SmeRowMajorMat<Scalar>, 0, OuterStride<>> G(g.data(), k, n, OuterStride<>(1));
        SmeColMajorMat<Scalar> D(m, n);
        D.noalias() = H * G;
        VERIFY_IS_APPROX(D, H.lazyProduct(G));
      }
  // Operands packed tight at the end of their buffer, so a read past the last row or column leaves it.
  const Index ps = Index(16 / sizeof(Scalar));
  for (Index m : {Index(2), Index(3), ps + 1, 2 * ps})
    for (Index n : {2, 3, 5, 7})
      for (Index k : {16, 17, 18, 19}) {
        std::vector<Scalar> a(std::size_t(m * k)), b(std::size_t(k * n));
        Map<SmeColMajorMat<Scalar>> A(a.data(), m, k);
        Map<SmeRowMajorMat<Scalar>> B(b.data(), k, n);
        A.setRandom();
        B.setRandom();
        SmeColMajorMat<Scalar> C = SmeColMajorMat<Scalar>::Zero(m, n);
        C.noalias() += A * B;
        VERIFY_IS_APPROX(C, A.lazyProduct(B));
        Map<SmeRowMajorMat<Scalar>> Ar(a.data(), m, k);
        Map<SmeColMajorMat<Scalar>> Bc(b.data(), k, n);
        C.noalias() = Ar * Bc;
        VERIFY_IS_APPROX(C, Ar.lazyProduct(Bc));
      }
}

// Whole products on the NEON kernel (sme_small_gemm) and just past its bound, every storage order of A, B and C.
template <typename Scalar, int AO, int BO, int CO>
static void verify_small_products(Index m, Index n, Index k) {
  using MA = Matrix<Scalar, Dynamic, Dynamic, AO>;
  using MB = Matrix<Scalar, Dynamic, Dynamic, BO>;
  using MC = Matrix<Scalar, Dynamic, Dynamic, CO>;
  const MA A = MA::Random(m, k);
  const MB B = MB::Random(k, n);
  const MC ref = A.lazyProduct(B);
  const Scalar alpha = nontrivial_alpha_impl<Scalar>::run();
  MC C = MC::Random(m, n);
  const MC C0 = C;
  C.noalias() = A * B;
  VERIFY_IS_APPROX(C, ref);
  C = C0;
  C.noalias() += alpha * A * B;
  VERIFY_IS_APPROX(C, C0 + alpha * ref);
  const SmeColMajorMat<Scalar> Ac = A.conjugate(), Bc = B.conjugate();
  C.noalias() = A.conjugate() * B;
  VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), SmeColMajorMat<Scalar>(Ac.lazyProduct(B)));
  C.noalias() = A * B.conjugate();
  VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), SmeColMajorMat<Scalar>(A.lazyProduct(Bc)));
  const MA At = A.adjoint();
  C.noalias() = At.adjoint() * B;
  VERIFY_IS_APPROX(C, ref);
}

template <typename Scalar>
static void test_small_products() {
#if defined(EIGEN_SME_NO_NEON_SMALL_BLOCKS) || defined(EIGEN_SME_FORCE_NEON_SMALL_BLOCKS) || defined(EIGEN_USE_BLAS)
  constexpr bool kSmallProductEnabled = false;
#else
  constexpr bool kSmallProductEnabled = true;
#endif
  const Index size = internal::sme_neon_small_product<Scalar>::value;
  const Index thin = Index(32 / sizeof(Scalar));
  // Both sides of the RHS bound for a thin LHS and of the size bound past it, a tall real LHS against three columns,
  // then squares and thin shapes on either side.
  std::vector<std::array<Index, 3>> shapes;
  const Index rhs_depth =
      numext::mini(Index(65536 / sizeof(Scalar)) / 37, Index(NumTraits<Scalar>::IsComplex ? 4096 : 128));
  shapes.push_back({thin, 37, rhs_depth});
  shapes.push_back({thin, 37, rhs_depth + 1});
  // Deeper than the blocking's kc, which packs into a buffer of its own.
  shapes.push_back({thin, 4, 2048});
  for (Index n : {Index(16), Index(37)}) {
    const Index k = size / ((thin + 1) * n);
    shapes.push_back({thin + 1, n, k});
    shapes.push_back({thin + 1, n, k + 1});
  }
  shapes.push_back({3 * sme_mr<Scalar>(), 3, 16});
  for (Index s : {12, 17, 20, 25, 29, 31}) shapes.push_back({s, s, s});
  for (auto s : std::vector<std::array<Index, 3>>{
           {2, 64, 48}, {4, 32, 64}, {8, 64, 64}, {32, 4, 32}, {64, 2, 64}, {16, 16, 32}, {12, 12, 64}, {3, 33, 97}})
    shapes.push_back(s);
  int on_neon = 0, past = 0;
  for (const auto& s : shapes) {
    const bool neon = internal::sme_small_gemm_wins<Scalar>(s[0], s[1], s[2]);
    on_neon += neon;
    past += !neon;
    verify_small_products<Scalar, ColMajor, ColMajor, ColMajor>(s[0], s[1], s[2]);
    verify_small_products<Scalar, RowMajor, ColMajor, ColMajor>(s[0], s[1], s[2]);
    verify_small_products<Scalar, ColMajor, RowMajor, ColMajor>(s[0], s[1], s[2]);
    verify_small_products<Scalar, RowMajor, RowMajor, ColMajor>(s[0], s[1], s[2]);
    verify_small_products<Scalar, ColMajor, ColMajor, RowMajor>(s[0], s[1], s[2]);
    verify_small_products<Scalar, RowMajor, RowMajor, RowMajor>(s[0], s[1], s[2]);
  }
  VERIFY_IS_EQUAL(on_neon > 0, kSmallProductEnabled);
  VERIFY(past > 0);
  // Sub-blocks of larger operands and a result with an inner stride.
  const Index m = thin + 3, n = 9, k = 33;
  VERIFY_IS_EQUAL(internal::sme_small_gemm_wins<Scalar>(m, n, k), kSmallProductEnabled && m * n * k <= size);
  const SmeColMajorMat<Scalar> Abig = SmeColMajorMat<Scalar>::Random(m + 5, k + 2);
  const SmeRowMajorMat<Scalar> Bbig = SmeRowMajorMat<Scalar>::Random(k + 3, n + 4);
  const auto A = Abig.block(2, 1, m, k);
  const auto B = Bbig.block(1, 2, k, n);
  SmeColMajorMat<Scalar> Cbig = SmeColMajorMat<Scalar>::Random(2 * m, n);
  SmeColMajorMat<Scalar> ref = Cbig;
  Map<SmeColMajorMat<Scalar>, 0, Stride<Dynamic, 2>> C(Cbig.data(), m, n, Stride<Dynamic, 2>(2 * m, 2));
  C.noalias() -= Scalar(3) * A * B;
  Map<SmeColMajorMat<Scalar>, 0, Stride<Dynamic, 2>>(ref.data(), m, n, Stride<Dynamic, 2>(2 * m, 2)) -=
      Scalar(3) * A.lazyProduct(B);
  VERIFY_IS_APPROX(Cbig, ref);
}

template <typename TestIndex>
static void test_small_product_overflow() {
  // With a 65536 threshold, both the two-factor and three-factor products can wrap a 32-bit index to zero.
  VERIFY(!internal::sme_small_gemm_wins<float>(TestIndex(65536), TestIndex(65536), TestIndex(1)));
  VERIFY(!internal::sme_small_gemm_wins<float>(TestIndex(256), TestIndex(256), TestIndex(65536)));
  VERIFY(!internal::sme_small_gemm_wins<float>(TestIndex(16), TestIndex(16), TestIndex(257)));
#if !defined(EIGEN_SME_NO_NEON_SMALL_BLOCKS) && !defined(EIGEN_SME_FORCE_NEON_SMALL_BLOCKS) && !defined(EIGEN_USE_BLAS)
  VERIFY(internal::sme_small_gemm_wins<float>(TestIndex(16), TestIndex(16), TestIndex(256)));
#endif
}

// Every tail width of a deep block, which the ZA transposer packs: a ColMajor RHS tail and, through a RowMajor
// product, a transposed LHS tail.
template <typename Scalar>
static void test_deep_tail_panels() {
  const Index MR = sme_mr<Scalar>(), NR = sme_nr<Scalar>();
  for (Index depth : {4 * NR, 6 * NR + 3}) {
    for (Index tail = 1; tail < NR; ++tail) {
      for (Index n : {tail, NR + tail}) {
        const SmeColMajorMat<Scalar> A = SmeColMajorMat<Scalar>::Random(MR + 3, depth);
        const SmeColMajorMat<Scalar> B = SmeColMajorMat<Scalar>::Random(depth, n);
        SmeColMajorMat<Scalar> C = SmeColMajorMat<Scalar>::Random(MR + 3, n);
        const SmeColMajorMat<Scalar> C0 = C;
        C.noalias() += A * B;
        VERIFY_IS_APPROX(C, C0 + A.lazyProduct(B));
        const SmeRowMajorMat<Scalar> Ar = SmeRowMajorMat<Scalar>::Random(n, depth);
        const SmeRowMajorMat<Scalar> Br = SmeRowMajorMat<Scalar>::Random(depth, MR + 3);
        SmeRowMajorMat<Scalar> Cr = SmeRowMajorMat<Scalar>::Random(n, MR + 3);
        const SmeRowMajorMat<Scalar> Cr0 = Cr;
        Cr.noalias() += Ar * Br;
        VERIFY_IS_APPROX(Cr, Cr0 + Ar.lazyProduct(Br));
      }
    }
  }
}

// ---------------------------------------------------------------------------
// Raw packed-buffer tests.
//
// The product tests above validate the packers only transitively through a full
// product, where the gebp_kernel can mask a mispack, or where a bug only shows
// at a specific SVL/region boundary.  The tests below call the SME packers
// directly and compare the packed buffer exactly against a scalar reference
// (equality is exact).  This pins every region -- in particular a dropped
// row-group in the two-pass trailing transpose, and a SYMM packer that copies
// the unused triangle instead of mirroring -- at whatever SVL the run uses.
// ---------------------------------------------------------------------------

// A distinctive marker for buffer cells the packer must leave untouched, and
// for the unused triangle of a lower-triangular operand. Random values live in
// [-1, 1], so it never collides with a real packed value. The complex marker's
// imaginary part differs from its real one, so a packer that fills only the
// real half of a split depth step is caught.
template <typename Scalar>
struct pack_sentinel_impl {
  static Scalar run() { return Scalar(98765); }
};
template <typename RealScalar>
struct pack_sentinel_impl<std::complex<RealScalar>> {
  static std::complex<RealScalar> run() { return std::complex<RealScalar>(RealScalar(98765), RealScalar(-54321)); }
};
template <typename Scalar>
static Scalar pack_sentinel() {
  return pack_sentinel_impl<Scalar>::run();
}

// ---------------------------------------------------------------------------
// NEON small-block path: sme_gebp_neon called directly on hand-packed panels,
// over every row, column and depth tail of the tile ladder, both conjugation
// flags, a non-trivial alpha (and a real one for complex scalars), panel-mode
// stride/offset and a general-stride C.
// ---------------------------------------------------------------------------
template <typename Scalar, bool ConjLhs, bool ConjRhs>
static void verify_neon_small_block(Index rows, Index cols, Index depth, const Scalar& alpha, bool panel_mode,
                                    bool strided_c) {
  const Index MR = sme_mr<Scalar>();
  const Index NR = sme_nr<Scalar>();
  const Index strideA = panel_mode ? depth + 3 : depth;
  const Index strideB = panel_mode ? depth + 5 : depth;
  const Index offsetA = panel_mode ? 2 : 0;
  const Index offsetB = panel_mode ? 1 : 0;
  const SmeColMajorMat<Scalar> A = SmeColMajorMat<Scalar>::Random(rows, depth);
  const SmeColMajorMat<Scalar> B = SmeColMajorMat<Scalar>::Random(depth, cols);

  SmeVector<Scalar> packedA =
      SmeVector<Scalar>::Constant((rows + MR) * strideA + offsetA * MR, pack_sentinel<Scalar>());
  SmeVector<Scalar> packedB =
      SmeVector<Scalar>::Constant((cols + NR) * strideB + offsetB * NR, pack_sentinel<Scalar>());
  for (Index i = 0; i < rows; i += MR) {
    const Index w = numext::mini(MR, rows - i);
    for (Index k = 0; k < depth; ++k)
      for (Index r = 0; r < w; ++r) set_packed(packedA.data() + i * strideA + offsetA * w, w, k, r, A(i + r, k));
  }
  for (Index j = 0; j < cols; j += NR) {
    const Index w = numext::mini(NR, cols - j);
    for (Index k = 0; k < depth; ++k)
      for (Index c = 0; c < w; ++c) set_packed(packedB.data() + j * strideB + offsetB * w, w, k, c, B(k, j + c));
  }

  const Index rs = strided_c ? 2 : 1;
  const Index cs = strided_c ? 3 * rows : rows;
  SmeVector<Scalar> storage = SmeVector<Scalar>::Random(cs * cols);
  SmeColMajorStridedMat<Scalar> C(storage.data(), rows, cols, Stride<Dynamic, Dynamic>(cs, rs));
  const SmeColMajorMat<Scalar> c_before = C;
  const SmeVector<Scalar> storage_before = storage;

  internal::sme_gebp_neon<Scalar, ConjLhs, ConjRhs, Index>(storage.data(), rs, cs, packedA.data(), packedB.data(), rows,
                                                           depth, cols, alpha, strideA, strideB, offsetA, offsetB);

  const SmeColMajorMat<Scalar> Ac = ConjLhs ? A.conjugate().eval() : A;
  const SmeColMajorMat<Scalar> Bc = ConjRhs ? B.conjugate().eval() : B;
  const SmeColMajorMat<Scalar> ref = c_before + alpha * Ac.lazyProduct(Bc);
  VERIFY(verifyProductUpdate(C, ref, alpha, Ac, Bc, c_before));

  // Cells the strided map skips stay untouched.
  if (strided_c) {
    for (Index j = 0; j < cols; ++j)
      for (Index i = 0; i < rows; ++i)
        VERIFY_IS_EQUAL(storage(j * cs + i * rs + 1), storage_before(j * cs + i * rs + 1));
  }
}

template <typename Scalar>
static void test_neon_small_blocks() {
  const Index PS = Index(internal::packet_traits<typename NumTraits<Scalar>::Real>::size);
  const Index MR = sme_mr<Scalar>();
  const Index NR = sme_nr<Scalar>();
  const Index widths[] = {1, 2, 3, PS, PS + 1, 2 * PS + 1, 4 * PS, 4 * PS + 3, MR - 1, MR, MR + 1, 2 * MR + 5};
  const Index col_widths[] = {1, 2, 3, 4, 5, 9, NR - 1, NR, NR + 1, 2 * NR + 2};
  const Index depths[] = {1, 2, 3, 4, 5, 8, 9, 35};
  const Scalar alpha = nontrivial_alpha_impl<Scalar>::run();
  const Scalar real_alpha = Scalar(numext::real(alpha));
  for (Index rows : widths) {
    for (Index cols : col_widths) {
      for (Index depth : depths) {
        const bool panel_mode = (rows + cols + depth) % 2 == 0;
        const bool strided_c = (rows + depth) % 3 == 0;
        verify_neon_small_block<Scalar, false, false>(rows, cols, depth, alpha, panel_mode, strided_c);
        if (NumTraits<Scalar>::IsComplex) {
          // A real alpha has its own store; conjugation is folded into the accumulators before it.
          verify_neon_small_block<Scalar, false, false>(rows, cols, depth, real_alpha, panel_mode, strided_c);
          verify_neon_small_block<Scalar, true, false>(rows, cols, depth, alpha, panel_mode, strided_c);
          verify_neon_small_block<Scalar, false, true>(rows, cols, depth, alpha, panel_mode, strided_c);
          verify_neon_small_block<Scalar, true, true>(rows, cols, depth, alpha, panel_mode, strided_c);
        }
      }
    }
  }
}

// Lower-triangular n x n operand plus the dense selfadjoint reference the packer
// must emit. The unused triangle is filled with the sentinel so a packer that
// copies the dense matrix and never mirrors fails VERIFY_IS_EQUAL.
// product_selfadjoint_matrix stores the valid triangle where row >= col
// (after the Upper/RowMajor xor), so the packer must read stored(row,col)
// below the diagonal and conj(stored(col,row)) above it. Averaging with the
// adjoint cancels the diagonal's imaginary part exactly, which is what a
// selfadjoint view defines it to be.
template <typename Scalar, int StorageOrder>
static void make_lower_stored_selfadjoint(Index n, Matrix<Scalar, Dynamic, Dynamic, StorageOrder>& stored,
                                          Matrix<Scalar, Dynamic, Dynamic, StorageOrder>& full) {
  full = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>::Random(n, n);
  full = ((full + full.adjoint()) * Scalar(0.5)).eval();
  stored = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>::Constant(n, n, pack_sentinel<Scalar>());
  for (Index i = 0; i < n; ++i)
    for (Index j = 0; j <= i; ++j) stored(i, j) = full(i, j);
}

// LHS SYMM packer: a square selfadjoint diagonal block of size n, packed into
// uniform mr-wide depth-major panels.  Reference: full(i+r, k).
template <typename Scalar, int StorageOrder>
static void verify_symm_pack_lhs(Index n) {
  const Index MR = sme_mr<Scalar>();
  Matrix<Scalar, Dynamic, Dynamic, StorageOrder> stored, full;
  make_lower_stored_selfadjoint<Scalar, StorageOrder>(n, stored, full);

  SmeVector<Scalar> packed = SmeVector<Scalar>::Constant(n * n, pack_sentinel<Scalar>());
  SmeVector<Scalar> ref = SmeVector<Scalar>::Constant(n * n, pack_sentinel<Scalar>());
  for (Index i = 0; i < n; i += MR) {
    const Index w = numext::mini(MR, n - i);
    for (Index k = 0; k < n; ++k)
      for (Index r = 0; r < w; ++r) set_packed(ref.data() + i * n, w, k, r, Scalar(full(i + r, k)));
  }

  internal::symm_pack_lhs<Scalar, Index, sme_mr<Scalar>(), 1, StorageOrder> pack;
  pack(packed.data(), stored.data(), stored.outerStride(), /*cols(depth)=*/n, /*rows=*/n);
  VERIFY_IS_EQUAL(packed, ref);
}

// RHS SYMM packer: a depth block [k2, k2 + depth) x cols columns of an N x N
// selfadjoint matrix, packed into nr-wide depth-major panels.  Reference:
// full(k2 + k, j + c).  A k2 > 0 offset makes the transposed region non-empty,
// so partial-width panels reach the two-pass transpose.
template <typename Scalar, int StorageOrder>
static void verify_symm_pack_rhs(Index N, Index depth, Index cols, Index k2) {
  eigen_assert(k2 + depth <= N && cols <= N);
  const Index NR = sme_nr<Scalar>();
  Matrix<Scalar, Dynamic, Dynamic, StorageOrder> stored, full;
  make_lower_stored_selfadjoint<Scalar, StorageOrder>(N, stored, full);

  SmeVector<Scalar> packed = SmeVector<Scalar>::Constant(cols * depth, pack_sentinel<Scalar>());
  SmeVector<Scalar> ref = SmeVector<Scalar>::Constant(cols * depth, pack_sentinel<Scalar>());
  for (Index j = 0; j < cols; j += NR) {
    const Index w = numext::mini(NR, cols - j);
    for (Index k = 0; k < depth; ++k)
      for (Index c = 0; c < w; ++c) set_packed(ref.data() + j * depth, w, k, c, Scalar(full(k2 + k, j + c)));
  }

  internal::symm_pack_rhs<Scalar, Index, sme_nr<Scalar>(), StorageOrder> pack;
  pack(packed.data(), stored.data(), stored.outerStride(), /*rows(depth)=*/depth, /*cols=*/cols, k2);
  VERIFY_IS_EQUAL(packed, ref);
}

template <typename Scalar>
static void test_symm_pack() {
  // The last panel width sweeps a range of partial widths; at each SVL the
  // two-pass trailing transpose (the if->loop fix) fires when a partial width
  // leaves a trailing row-group remainder in (svl, 2*svl).  The spread below
  // hits that for svl in {2, 4, 8, 16, 32, 64} -- fp32 SVL 128..2048 and the
  // fp64 lane counts, which are half of those.
  const int sizes[] = {1, 5, 7, 17, 31, 32, 33, 37, 39, 45, 48, 49, 55, 57, 63, 64, 65, 79, 96, 97};
  for (int n : sizes) {
    verify_symm_pack_lhs<Scalar, ColMajor>(n);
    verify_symm_pack_lhs<Scalar, RowMajor>(n);
    // RHS, single depth block anchored at the diagonal (k2 == 0).
    verify_symm_pack_rhs<Scalar, ColMajor>(n, n, n, 0);
    verify_symm_pack_rhs<Scalar, RowMajor>(n, n, n, 0);
  }

  // RHS depth blocks offset from the diagonal (k2 > 0): the transposed region is
  // non-empty, so the RowMajor operand drives partial-width panels through the
  // two-pass transpose and the ColMajor operand through the partial copy.
  struct RhsCase {
    int N, depth, cols, k2;
  };
  const RhsCase rhs_cases[] = {
      {100, 32, 39, 16}, {100, 24, 39, 32}, {100, 40, 64, 8}, {100, 39, 39, 33}, {128, 57, 57, 40}, {128, 33, 45, 60},
  };
  for (const RhsCase& c : rhs_cases) {
    verify_symm_pack_rhs<Scalar, ColMajor>(c.N, c.depth, c.cols, c.k2);
    verify_symm_pack_rhs<Scalar, RowMajor>(c.N, c.depth, c.cols, c.k2);
  }

  // Full-panel direct regions [s_end, depth) of a ColMajor operand start at s_end = j + w - k2, off the depth chunk
  // grid, with real-view lengths 4*svl*a + 2*svl + t that pass through the four-chunk (64-bit), paired and single-tile
  // transposer loops in turn. The panel width w is nr and, where smaller, the 2*svl the full-panel path needs.
  const int svl = internal::sme_packet_traits<typename NumTraits<Scalar>::Real>::size();
  const int NR = sme_nr<Scalar>();
  const int scale = NumTraits<Scalar>::IsComplex ? 2 : 1;
  for (int real_range : {10 * svl + 3, 14 * svl + 1}) {
    const int range = (real_range + scale - 1) / scale;
    for (int k2 : {5, NR - 3}) {
      const int depth = NR - k2 + range, N = numext::maxi(k2 + depth, 2 * NR + 1);
      verify_symm_pack_rhs<Scalar, ColMajor>(N, depth, 2 * NR + 1, k2);
      verify_symm_pack_rhs<Scalar, RowMajor>(N, depth, 2 * NR + 1, k2);
    }
    if (svl >= 4 && 2 * svl < NR) {
      const int w = 2 * svl, k2 = w - 3, depth = 3 + range;
      verify_symm_pack_rhs<Scalar, ColMajor>(k2 + depth, depth, w, k2);
    }
  }
}

// ---------------------------------------------------------------------------
// Mapper-based packing fallback (sme_pack_{lhs,rhs}_fallback).
//
// Taken by by-value tensor sub-mappers and inner-strided blas mappers, which
// the raw pointer + stride packers cannot walk.  The product suite reaches these
// only through tensor contractions (nightly, SVL=512) and TriangularSolver
// (random sizes), so here we drive the fallback directly with both mapper
// families and compare the packed buffer against a scalar reference.
// ---------------------------------------------------------------------------

// Minimal stand-ins for by-value sub-mappers. ColMajor packets advance the
// first index, while RowMajor packets follow the normal storage-inner second
// index. operator() returns by value so both take the no-direct-access dispatch.
template <typename Scalar>
struct ByValueColMajorLhsMapper {
  const Scalar* data;
  Index stride;  // element(i, k) = data[i + k * stride], contiguous in i
  Scalar operator()(Index i, Index k) const { return data[i + k * stride]; }
  template <typename Packet>
  EIGEN_ALWAYS_INLINE Packet loadPacket(Index i, Index k) const {
    return internal::ploadu<Packet>(data + i + k * stride);
  }
};

template <typename Scalar>
struct ByValueRowMajorLhsMapper {
  const Scalar* data;
  Index stride;  // element(i, k) = data[i * stride + k], contiguous in k
  Scalar operator()(Index i, Index k) const { return data[i * stride + k]; }
  template <typename Packet>
  EIGEN_ALWAYS_INLINE Packet loadPacket(Index i, Index k) const {
    return internal::ploadu<Packet>(data + i * stride + k);
  }
};

template <typename Scalar>
struct ByValueColMajorRhsMapper {
  struct LinearMapper {
    const Scalar* p;  // &element(0, col); contiguous in depth
    Scalar operator()(Index k) const { return p[k]; }
    template <typename Packet>
    EIGEN_ALWAYS_INLINE Packet loadPacket(Index k) const {
      return internal::ploadu<Packet>(p + k);
    }
  };
  const Scalar* data;
  Index stride;  // element(k, col) = data[k + col * stride], contiguous in k
  Scalar operator()(Index k, Index col) const { return data[k + col * stride]; }
  LinearMapper getLinearMapper(Index k, Index col) const { return LinearMapper{data + k + col * stride}; }
};

template <typename Scalar>
struct ByValueRowMajorRhsMapper {
  struct LinearMapper {
    const Scalar* p;  // &element(row, col); packet offsets advance columns
    Scalar operator()(Index offset) const { return p[offset]; }
    template <typename Packet>
    EIGEN_ALWAYS_INLINE Packet loadPacket(Index offset) const {
      return internal::ploadu<Packet>(p + offset);
    }
  };
  const Scalar* data;
  Index stride;  // element(k, col) = data[k * stride + col], contiguous in col
  Scalar operator()(Index k, Index col) const { return data[k * stride + col]; }
  LinearMapper getLinearMapper(Index k, Index col) const { return LinearMapper{data + k * stride + col}; }
};

// Length of the packed LHS buffer for the given panel-mode layout, matching the
// dst_panel formula in sme_pack_lhs_fallback.
template <bool PanelMode>
static Index packed_len(Index outer, Index depth, Index unit, Index dst_stride, Index dst_offset) {
  // `outer` is rows (LHS) or cols (RHS); `unit` is the panel width mr or nr.
  if (!PanelMode) return outer * depth;
  Index end = 0;
  for (Index i = 0; i < outer; i += unit) {
    const Index w = numext::mini(unit, outer - i);
    end = numext::maxi(end, i * dst_stride + dst_offset * w + depth * w);
  }
  return end;
}

template <typename Scalar, bool PanelMode, bool Conjugate, typename MatrixType>
static void fill_lhs_ref(SmeVector<Scalar>& ref, const MatrixType& V, Index rows, Index depth, Index dst_stride,
                         Index dst_offset) {
  const Index MR = sme_mr<Scalar>();
  ref.setConstant(pack_sentinel<Scalar>());
  for (Index i = 0; i < rows; i += MR) {
    const Index w = numext::mini(MR, rows - i);
    const Index base = PanelMode ? i * dst_stride + dst_offset * w : i * depth;
    for (Index k = 0; k < depth; ++k)
      for (Index r = 0; r < w; ++r) {
        const Scalar v = V(i + r, k);
        set_packed(ref.data() + base, w, k, r, Conjugate ? numext::conj(v) : v);
      }
  }
}

template <typename Scalar, bool PanelMode, bool Conjugate, typename MatrixType>
static void fill_rhs_ref(SmeVector<Scalar>& ref, const MatrixType& V, Index cols, Index depth, Index dst_stride,
                         Index dst_offset) {
  const Index NR = sme_nr<Scalar>();
  ref.setConstant(pack_sentinel<Scalar>());
  for (Index j = 0; j < cols; j += NR) {
    const Index w = numext::mini(NR, cols - j);
    const Index base = PanelMode ? j * dst_stride + dst_offset * w : j * depth;
    for (Index k = 0; k < depth; ++k)
      for (Index c = 0; c < w; ++c) {
        const Scalar v = V(k, j + c);
        set_packed(ref.data() + base, w, k, c, Conjugate ? numext::conj(v) : v);
      }
  }
}

// Every packer check below is the same nine lines -- fill a sentinel buffer,
// build the scalar reference, run the packer, compare exactly -- over a
// different mapper. The mapper is what each check is really about, so it is
// built by the caller and the rest lives here once.
template <typename Scalar, int StorageOrder, bool PanelMode, bool Conjugate, typename Mapper, typename MatrixType>
static void check_packed_lhs(const Mapper& mapper, const MatrixType& V, Index rows, Index depth) {
  const Index dst_stride = PanelMode ? depth + 5 : 0;
  const Index dst_offset = PanelMode ? 3 : 0;
  const Index len = packed_len<PanelMode>(rows, depth, sme_mr<Scalar>(), dst_stride, dst_offset);
  SmeVector<Scalar> packed = SmeVector<Scalar>::Constant(len, pack_sentinel<Scalar>());
  SmeVector<Scalar> ref(len);
  fill_lhs_ref<Scalar, PanelMode, Conjugate>(ref, V, rows, depth, dst_stride, dst_offset);

  internal::gemm_pack_lhs<Scalar, Index, Mapper, sme_mr<Scalar>(), 1, typename internal::packet_traits<Scalar>::type,
                          StorageOrder, Conjugate, PanelMode>
      pack;
  pack(packed.data(), mapper, depth, rows, dst_stride, dst_offset);
  VERIFY_IS_EQUAL(packed, ref);
}

template <typename Scalar, int StorageOrder, bool PanelMode, bool Conjugate, typename Mapper, typename MatrixType>
static void check_packed_rhs(const Mapper& mapper, const MatrixType& V, Index depth, Index cols) {
  const Index dst_stride = PanelMode ? depth + 5 : 0;
  const Index dst_offset = PanelMode ? 3 : 0;
  const Index len = packed_len<PanelMode>(cols, depth, sme_nr<Scalar>(), dst_stride, dst_offset);
  SmeVector<Scalar> packed = SmeVector<Scalar>::Constant(len, pack_sentinel<Scalar>());
  SmeVector<Scalar> ref(len);
  fill_rhs_ref<Scalar, PanelMode, Conjugate>(ref, V, cols, depth, dst_stride, dst_offset);

  internal::gemm_pack_rhs<Scalar, Index, Mapper, sme_nr<Scalar>(), StorageOrder, Conjugate, PanelMode> pack;
  pack(packed.data(), mapper, depth, cols, dst_stride, dst_offset);
  VERIFY_IS_EQUAL(packed, ref);
}

// The raw pointer + stride packers, taken whenever the mapper grants direct
// unit-inner-stride access.
template <typename Scalar, int StorageOrder, bool PanelMode, bool Conjugate>
static void verify_pack_lhs_direct(Index rows, Index depth) {
  using MatrixType = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>;
  using Mapper = internal::const_blas_data_mapper<Scalar, Index, StorageOrder>;
  MatrixType V = MatrixType::Random(rows, depth);
  check_packed_lhs<Scalar, StorageOrder, PanelMode, Conjugate>(Mapper(V.data(), V.outerStride()), V, rows, depth);
}

template <typename Scalar, int StorageOrder, bool PanelMode, bool Conjugate>
static void verify_pack_rhs_direct(Index depth, Index cols) {
  using MatrixType = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>;
  using Mapper = internal::const_blas_data_mapper<Scalar, Index, StorageOrder>;
  MatrixType V = MatrixType::Random(depth, cols);
  check_packed_rhs<Scalar, StorageOrder, PanelMode, Conjugate>(Mapper(V.data(), V.outerStride()), V, depth, cols);
}

template <bool Conjugate, typename Scalar>
EIGEN_DONT_INLINE __arm_locally_streaming static void sme_copy_range(Scalar* dst, const Scalar* src, Index stride,
                                                                     Index k0, Index k1, int width) {
  internal::sve_copy_panel_range<Conjugate>(dst, src, stride, k0, k1, width);
}

// Every partial-vector width, including copies ending with an unshifted tail and exact vector boundaries.
template <typename Scalar>
static void test_copy_panel_range() {
  const int svl = internal::sme_packet_traits<typename NumTraits<Scalar>::Real>::size();
  for (int width = 1; width <= 2 * svl; ++width) {
    for (Index k0 : {0, 1, svl - 1}) {
      for (Index depth : {0, 1, 2, 3, svl - 1, svl, svl + 1, 2 * svl + 1}) {
        const Index k1 = k0 + depth;
        SmeRowMajorMat<Scalar> src = SmeRowMajorMat<Scalar>::Random(k1, width + 3);
        SmeVector<Scalar> packed = SmeVector<Scalar>::Constant(k1 * width + 2, pack_sentinel<Scalar>());
        SmeVector<Scalar> ref = packed;
        sme_copy_range<false>(packed.data() + 1, src.data(), src.outerStride(), k0, k1, width);
        for (Index k = k0; k < k1; ++k)
          for (Index c = 0; c < width; ++c) set_packed(ref.data() + 1, Index(width), k, c, Scalar(src(k, c)));
        VERIFY_IS_EQUAL(packed, ref);
        packed.setConstant(pack_sentinel<Scalar>());
        sme_copy_range<true>(packed.data() + 1, src.data(), src.outerStride(), k0, k1, width);
        for (Index k = k0; k < k1; ++k)
          for (Index c = 0; c < width; ++c)
            set_packed(ref.data() + 1, Index(width), k, c, Scalar(numext::conj(src(k, c))));
        VERIFY_IS_EQUAL(packed, ref);
      }
    }
  }
}

// Inner-strided blas mapper LHS: element(i, k) laid out with inner stride
// `incr`.  ColMajor takes the vectorised gather path; RowMajor takes the scalar
// path (its packets would run along depth, not rows).
template <typename Scalar, int StorageOrder, bool PanelMode, bool Conjugate>
static void verify_fallback_lhs_strided(Index rows, Index depth, Index incr) {
  using Mapper = internal::blas_data_mapper<Scalar, Index, StorageOrder, Unaligned, Dynamic>;
  Matrix<Scalar, Dynamic, Dynamic> V = Matrix<Scalar, Dynamic, Dynamic>::Random(rows, depth);
  const Index mstride = (StorageOrder == ColMajor ? rows : depth) * incr;
  SmeVector<Scalar> buf = SmeVector<Scalar>::Zero((StorageOrder == ColMajor ? depth : rows) * mstride + incr);
  for (Index k = 0; k < depth; ++k)
    for (Index i = 0; i < rows; ++i)
      buf[StorageOrder == ColMajor ? i * incr + k * mstride : k * incr + i * mstride] = V(i, k);
  check_packed_lhs<Scalar, StorageOrder, PanelMode, Conjugate>(Mapper(buf.data(), mstride, incr), V, rows, depth);
}

// By-value LHS mappers exercise both packet directions. RowMajor must stay
// scalar because its packets advance depth rather than rows.
template <typename Scalar, int StorageOrder, bool PanelMode, bool Conjugate>
static void verify_fallback_lhs_byvalue(Index rows, Index depth) {
  using MatrixType = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>;
  using Mapper = typename std::conditional<StorageOrder == ColMajor, ByValueColMajorLhsMapper<Scalar>,
                                           ByValueRowMajorLhsMapper<Scalar>>::type;
  MatrixType V = MatrixType::Random(rows, depth);
  check_packed_lhs<Scalar, StorageOrder, PanelMode, Conjugate>(Mapper{V.data(), V.outerStride()}, V, rows, depth);
}

// Inner-strided blas mapper RHS: element(k, col) with inner stride `incr`.
// ColMajor takes the vectorised transpose path; RowMajor takes the scalar path.
template <typename Scalar, int StorageOrder, bool PanelMode, bool Conjugate>
static void verify_fallback_rhs_strided(Index depth, Index cols, Index incr) {
  using Mapper = internal::blas_data_mapper<Scalar, Index, StorageOrder, Unaligned, Dynamic>;
  Matrix<Scalar, Dynamic, Dynamic> V = Matrix<Scalar, Dynamic, Dynamic>::Random(depth, cols);
  const Index mstride = (StorageOrder == ColMajor ? depth : cols) * incr;
  SmeVector<Scalar> buf = SmeVector<Scalar>::Zero((StorageOrder == ColMajor ? cols : depth) * mstride + incr);
  for (Index col = 0; col < cols; ++col)
    for (Index k = 0; k < depth; ++k)
      buf[StorageOrder == ColMajor ? k * incr + col * mstride : col * incr + k * mstride] = V(k, col);
  check_packed_rhs<Scalar, StorageOrder, PanelMode, Conjugate>(Mapper(buf.data(), mstride, incr), V, depth, cols);
}

// By-value RHS mappers likewise cover both packet directions. RowMajor packets
// advance columns, so the depth-oriented transpose fallback must stay scalar.
template <typename Scalar, int StorageOrder, bool PanelMode, bool Conjugate>
static void verify_fallback_rhs_byvalue(Index depth, Index cols) {
  using MatrixType = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>;
  using Mapper = typename std::conditional<StorageOrder == ColMajor, ByValueColMajorRhsMapper<Scalar>,
                                           ByValueRowMajorRhsMapper<Scalar>>::type;
  MatrixType V = MatrixType::Random(depth, cols);
  check_packed_rhs<Scalar, StorageOrder, PanelMode, Conjugate>(Mapper{V.data(), V.outerStride()}, V, depth, cols);
}

// Each packer is swept over all four (PanelMode, Conjugate) combinations.
// Conjugate=true is the identity on real scalars; for complex ones it negates
// the packed imaginary half, and reaches the direct packers from the SYMM
// above-diagonal transposed pack and the fallbacks from an inner-strided
// selfadjoint operand.
template <typename Scalar, int SO>
static void sweep_pack_direct(Index n, Index depth) {
  verify_pack_lhs_direct<Scalar, SO, false, false>(n, depth);
  verify_pack_lhs_direct<Scalar, SO, false, true>(n, depth);
  verify_pack_lhs_direct<Scalar, SO, true, false>(n, depth);
  verify_pack_lhs_direct<Scalar, SO, true, true>(n, depth);
  verify_pack_rhs_direct<Scalar, SO, false, false>(depth, n);
  verify_pack_rhs_direct<Scalar, SO, false, true>(depth, n);
  verify_pack_rhs_direct<Scalar, SO, true, false>(depth, n);
  verify_pack_rhs_direct<Scalar, SO, true, true>(depth, n);
}

template <typename Scalar, int SO>
static void sweep_fallback_strided(Index n, Index depth, Index incr) {
  verify_fallback_lhs_strided<Scalar, SO, false, false>(n, depth, incr);
  verify_fallback_lhs_strided<Scalar, SO, false, true>(n, depth, incr);
  verify_fallback_lhs_strided<Scalar, SO, true, false>(n, depth, incr);
  verify_fallback_lhs_strided<Scalar, SO, true, true>(n, depth, incr);
  verify_fallback_rhs_strided<Scalar, SO, false, false>(depth, n, incr);
  verify_fallback_rhs_strided<Scalar, SO, false, true>(depth, n, incr);
  verify_fallback_rhs_strided<Scalar, SO, true, false>(depth, n, incr);
  verify_fallback_rhs_strided<Scalar, SO, true, true>(depth, n, incr);
}

template <typename Scalar, int SO>
static void sweep_fallback_byvalue(Index n, Index depth) {
  verify_fallback_lhs_byvalue<Scalar, SO, false, false>(n, depth);
  verify_fallback_lhs_byvalue<Scalar, SO, false, true>(n, depth);
  verify_fallback_lhs_byvalue<Scalar, SO, true, false>(n, depth);
  verify_fallback_lhs_byvalue<Scalar, SO, true, true>(n, depth);
  verify_fallback_rhs_byvalue<Scalar, SO, false, false>(depth, n);
  verify_fallback_rhs_byvalue<Scalar, SO, false, true>(depth, n);
  verify_fallback_rhs_byvalue<Scalar, SO, true, false>(depth, n);
  verify_fallback_rhs_byvalue<Scalar, SO, true, true>(depth, n);
}

template <typename Scalar>
static void test_pack_direct() {
  const int TILE = sme_tile<Scalar>();
  const int MR = sme_mr<Scalar>();
  const int NR = sme_nr<Scalar>();
  // Widths around the tile side and both panel widths, which differ for
  // complex scalars.
  const int widths[] = {1, TILE - 1, TILE, TILE + 1, MR, MR + 1, NR, NR + 1, 2 * NR + 1};
  // Depths up to sme_neon_max_depth (24 / 16 / 16 / 8 by scalar) reach pack_neon,
  // the deeper ones the streaming pack_direct and its predicated tails.
  const int depths[] = {1, 3, 8, 9, 17, 25, 27, 35};
  for (int d : depths) {
    for (int n : widths) {
      sweep_pack_direct<Scalar, ColMajor>(n, d);
      sweep_pack_direct<Scalar, RowMajor>(n, d);
    }
  }
  // Four-panel LHS copies with leftover panels and tails, over several transposing depth chunks.
  for (int d : {35, 67, 97}) {
    for (int n : {4 * MR, 5 * MR + 1, 8 * MR + 3}) {
      sweep_pack_direct<Scalar, ColMajor>(n, d);
      sweep_pack_direct<Scalar, RowMajor>(n, d);
    }
  }
  // Depths at and one past whole multiples of four depth chunks, where the full-panel transposer's loops hand over to
  // the next and to the single-tile tail; a complex panel transposes a real view of twice the depth.
  const int svl = internal::sme_packet_traits<typename NumTraits<Scalar>::Real>::size();
  const int deep = (NumTraits<Scalar>::IsComplex ? 16 : 32) * svl;
  std::vector<int> deep_depths = {deep - 1, deep + 1};
  for (int j : {2, 3, 8, 9}) {
    deep_depths.push_back(4 * svl * j);
    deep_depths.push_back(4 * svl * j + 1);
  }
  // Deduplicate deep + 1 == 4*svl*8 + 1 for real scalars.
  std::sort(deep_depths.begin(), deep_depths.end());
  deep_depths.erase(std::unique(deep_depths.begin(), deep_depths.end()), deep_depths.end());
  // The paired transposer takes only unconjugated 2*svl-wide panels. At SVL 512 that is NR for every scalar and MR for
  // all but complex<float>. Below it, eligible real tails reach the paired transposer through tail_pack once deep
  // enough, and complex<float> has MR == 2*svl at SVL 256. Complex tails take the NEON tail_transpose_pack, so
  // complex<double> at SVL 256 reaches the four-chunk loop only through test_symm_pack's 2*svl-wide direct regions.
  // Real tails between svl and 2*svl wide store whole packed lines, through the copy or the ZA transposer.
  std::vector<int> deep_widths = {MR, MR + 1, NR, 2 * NR + 1, 2 * svl, svl + 1, 2 * svl - 1};
  std::sort(deep_widths.begin(), deep_widths.end());
  deep_widths.erase(std::unique(deep_widths.begin(), deep_widths.end()), deep_widths.end());
  for (int d : deep_depths) {
    for (int n : deep_widths) {
      sweep_pack_direct<Scalar, ColMajor>(n, d);
      sweep_pack_direct<Scalar, RowMajor>(n, d);
    }
  }
}

template <typename Scalar>
static void test_mapper_fallback() {
  const int MR = sme_mr<Scalar>();
  const int widths[] = {4, 5, MR, MR + 1, 2 * MR + 1};  // rows/cols around 4 and the panel width
  const int depths[] = {1, 3, 8, 35};                   // depth remainders 1..3 and larger
  for (int n : widths) {
    for (int d : depths) {
      // RowMajor mappers take the scalar path: their packets advance the index
      // the fallback does not transpose.
      for (int incr : {2, 3}) {
        sweep_fallback_strided<Scalar, ColMajor>(n, d, incr);
        sweep_fallback_strided<Scalar, RowMajor>(n, d, incr);
      }
      sweep_fallback_byvalue<Scalar, ColMajor>(n, d);
      sweep_fallback_byvalue<Scalar, RowMajor>(n, d);
    }
  }
}

// ---------------------------------------------------------------------------
// Product-level coverage, swept relative to the scalar's own block width.
// ---------------------------------------------------------------------------

// Sizes that land just on and off the block tails and the intra-block ZA-tile
// splits, for both panel widths (they differ for complex scalars) and the tile
// side at the SVL=512 design point.
template <typename Scalar>
static std::vector<int> sme_edge_sizes() {
  const int T = sme_tile<Scalar>();
  const int MR = sme_mr<Scalar>();
  const int NR = sme_nr<Scalar>();
  std::vector<int> sizes = {T / 2,  T + T / 2, 1,      T - 1,      T,      T + 1,      MR - 1,     MR,     MR + 1,
                            NR - 1, NR,        NR + 1, 2 * MR - 1, 2 * MR, 2 * MR + 1, 2 * NR - 1, 2 * NR, 2 * NR + 1};
  std::sort(sizes.begin(), sizes.end());
  sizes.erase(std::unique(sizes.begin(), sizes.end()), sizes.end());
  return sizes;
}

// The dimension-sum and output-area bounds in GeneralProduct.h route a small
// square product to the coeff-based evaluator, so the n x n sweeps above stop
// reaching the SME kernel once n gets small -- which is exactly where the
// intra-block ZA tile-split edges (T-1, T, T+1) live for double and
// complex<double>.  Repeat the same widths over a depth that clears both
// bounds, so every m and n edge still reaches the kernel, and check the
// conjugated forms too: at those sizes test_conjugated_products compares a
// lazyProduct against a lazyProduct and cannot fail.
template <typename Scalar>
static void test_edge_sizes_deep_k() {
  const int depth = 4 * sme_mr<Scalar>();
  for (int n : sme_edge_sizes<Scalar>()) {
    if (n < 2) continue;
    const SmeColMajorMat<Scalar> A = SmeColMajorMat<Scalar>::Random(n, depth);
    const SmeColMajorMat<Scalar> B = SmeColMajorMat<Scalar>::Random(depth, n);
    const SmeColMajorMat<Scalar> Aa = A.adjoint().eval();
    const SmeColMajorMat<Scalar> Ba = B.adjoint().eval();
    const SmeColMajorMat<Scalar> c_before = SmeColMajorMat<Scalar>::Random(n, n);
    const SmeColMajorMat<Scalar> expected = c_before + A.lazyProduct(B);

    SmeColMajorMat<Scalar> C = c_before;
    C.noalias() += A * B;
    VERIFY_IS_APPROX(C, expected);

    C = c_before;
    C.noalias() += Aa.adjoint() * B;
    VERIFY_IS_APPROX(C, expected);

    C = c_before;
    C.noalias() += A * Ba.adjoint();
    VERIFY_IS_APPROX(C, expected);
  }
}

template <typename Scalar>
static void test_conjugated_products() {
  for (int n : sme_edge_sizes<Scalar>()) {
    verify_conjugated_products<Scalar, SmeColMajorMat<Scalar>, SmeColMajorMat<Scalar>, SmeColMajorMat<Scalar>>(n);
    verify_conjugated_products<Scalar, SmeRowMajorMat<Scalar>, SmeColMajorMat<Scalar>, SmeColMajorMat<Scalar>>(n);
    verify_conjugated_products<Scalar, SmeColMajorMat<Scalar>, SmeRowMajorMat<Scalar>, SmeRowMajorMat<Scalar>>(n);
    verify_conjugated_products<Scalar, SmeRowMajorMat<Scalar>, SmeRowMajorMat<Scalar>, SmeColMajorMat<Scalar>>(n);
  }
}

template <typename Scalar>
static void test_products() {
  const int MR = sme_mr<Scalar>();

  // Square edge cases around the block and tile boundaries.
  for (int n : sme_edge_sizes<Scalar>()) product(SmeColMajorMat<Scalar>(n, n));

  // Thin / wide rectangular cases (M x 1, 1 x N) and non-square cases that
  // exercise tail paths for both M and N.
  product(SmeColMajorMat<Scalar>(MR, 1));
  product(SmeColMajorMat<Scalar>(1, MR));
  product(SmeColMajorMat<Scalar>(1, 2 * MR));
  product(SmeColMajorMat<Scalar>(2 * MR, 1));
  product(SmeColMajorMat<Scalar>(MR + 1, 2 * MR + 1));
  product(SmeColMajorMat<Scalar>(2 * MR + 1, MR + 1));
  product(SmeColMajorMat<Scalar>(MR - 1, 2 * MR - 1));
  product(SmeColMajorMat<Scalar>(MR + 1, 7));
  product(SmeColMajorMat<Scalar>(7, MR + 1));
  product(SmeColMajorMat<Scalar>(4 * MR, 3));
  product(SmeColMajorMat<Scalar>(3, 4 * MR));

  test_deep_k_split<Scalar>();
  test_edge_sizes_deep_k<Scalar>();

  // Random sizes
  for (int i = 0; i < g_repeat; i++) {
    product(SmeColMajorMat<Scalar>(internal::random<int>(1, EIGEN_TEST_MAX_SIZE),
                                   internal::random<int>(1, EIGEN_TEST_MAX_SIZE)));
  }

  // Exercise the RowMajor packers and RowMajor result path.  When the input
  // MatrixType is RowMajor, product() instantiates m1/m2/m3/res in RowMajor,
  // so every matrix-matrix product in the suite flows through:
  //   - the RowMajor LHS packer (gemm_pack_lhs<..., RowMajor>)
  //   - the RowMajor RHS packer (gemm_pack_rhs<..., RowMajor>)
  //   - the RowMajor-C dispatch in GeneralMatrixMatrix.h (which transposes
  //     the computation: C^T = B^T * A^T).
  for (int n : sme_edge_sizes<Scalar>()) {
    if (n > 1) product(SmeRowMajorMat<Scalar>(n, n));
  }
  product(SmeRowMajorMat<Scalar>(MR + 1, 2 * MR + 1));
  product(SmeRowMajorMat<Scalar>(2 * MR + 1, MR + 1));
  for (int i = 0; i < g_repeat; i++) {
    product(SmeRowMajorMat<Scalar>(internal::random<int>(1, EIGEN_TEST_MAX_SIZE),
                                   internal::random<int>(1, EIGEN_TEST_MAX_SIZE)));
  }

  // Exercise the general-stride branch of sme_store_za_tile: fires when both
  // C_stride_row != 1 and C_stride_col != 1, e.g. a Map<Matrix> with an
  // explicit non-unit inner stride.  product.h never builds such a result, so
  // without this subtest the scalar-scatter path is effectively untested.
  for (int n : sme_edge_sizes<Scalar>()) {
    if (n < 2) continue;
    test_general_strided_result<Scalar, SmeColMajorMat<Scalar>>(n);
    test_general_strided_result<Scalar, SmeRowMajorMat<Scalar>>(n);
    test_rowmajor_strided_result<Scalar>(n);
  }

  // Row-LHS x Row-RHS -> Col-C: the one LHS/RHS/C storage combination that
  // product.h's transpose-style expressions never build directly (it always
  // flips one side of the multiplication).  The code paths are the same as
  // other combinations via Eigen's dispatch, but exercise them explicitly.
  for (int n : sme_edge_sizes<Scalar>()) {
    if (n < 2) continue;
    SmeRowMajorMat<Scalar> A = SmeRowMajorMat<Scalar>::Random(n, n);
    SmeRowMajorMat<Scalar> B = SmeRowMajorMat<Scalar>::Random(n, n);
    SmeColMajorMat<Scalar> C = SmeColMajorMat<Scalar>::Zero(n, n);
    C.noalias() += A * B;
    VERIFY_IS_APPROX(C, (A.lazyProduct(B)).eval());
  }
}

// Folding C into the complex alpha chain, (c_re + re*ar) - im*ai, overflows
// for acc = (h,h), alpha = (1,1), C = (3h,0) although the exact result
// (3h,2h) is finite (#3139 review).  Widths straddle the one-vector slice
// width, which selects the store form; a real alpha takes the fused store.
// RowMajor GEMM transposes the computation, so the packed kernel is also
// called directly to reach the horizontal slice store.
template <typename Real, int Order>
static void test_complex_alpha_range() {
  using Scalar = std::complex<Real>;
  using Mat = Matrix<Scalar, Dynamic, Dynamic, Order>;
  const int svl = internal::sme_packet_traits<Real>::size();
  const Real h = std::ldexp(Real(1), std::numeric_limits<Real>::max_exponent - 2);
  const Scalar initial(3 * h, 0);
  const Scalar alphas[] = {Scalar(1, 1), Scalar(-1, 0)};
  const Scalar expected[] = {Scalar(3 * h, 2 * h), Scalar(2 * h, -h)};
  // At depth 8 the products with 5 < n <= sme_neon_shallow_width run on the NEON kernel (smaller ones are
  // coefficient-based), at 65 on ZA.
  for (Index depth : {Index(8), Index(65)}) {
    for (int n : {svl / 2 - 1, svl / 2, svl / 2 + 1, svl, svl + 1}) {
      if (n < 2) continue;
      Mat a = Mat::Zero(n, depth), b = Mat::Zero(depth, n);
      a.col(0).setConstant(Scalar(h, h));
      b.row(0).setOnes();
      const Mat ac = a.conjugate(), bc = b.conjugate();
      for (int t = 0; t < 2; ++t) {
        const Scalar alpha = alphas[t];
        for (int conjugation = 0; conjugation < 4; ++conjugation) {
          Mat c = Mat::Constant(n, n, initial);
          if (conjugation == 0) c.noalias() += alpha * (a * b);
          if (conjugation == 1) c.noalias() += alpha * (ac.conjugate() * b);
          if (conjugation == 2) c.noalias() += alpha * (a * bc.conjugate());
          if (conjugation == 3) c.noalias() += alpha * (ac.conjugate() * bc.conjugate());
          for (Index j = 0; j < n; ++j) {
            for (Index i = 0; i < n; ++i) {
              VERIFY((std::isfinite)(c(i, j).real()) && (std::isfinite)(c(i, j).imag()));
              VERIFY(c(i, j) == expected[t]);
            }
          }
        }
        // One packed panel through each kernel at every depth; a larger runtime SVL may exceed its static width.
        if (n <= numext::mini(sme_mr<Scalar>(), sme_nr<Scalar>())) {
          std::vector<Scalar> pa(n * depth), pb(n * depth);
          for (Index k = 0; k < depth; ++k) {
            for (Index i = 0; i < n; ++i) {
              set_packed(pa.data(), Index(n), k, i, a(i, k));
              set_packed(pb.data(), Index(n), k, i, b(k, i));
            }
          }
          const Index rs = Order == RowMajor ? n : 1, cs = Order == RowMajor ? 1 : n;
          for (int kernel = 0; kernel < 2; ++kernel) {
            Mat c = Mat::Constant(n, n, initial);
            if (kernel == 0) {
              internal::sme_gebp_dispatch<Scalar, false, false>::run(c.data(), rs, cs, pa.data(), pb.data(), Index(n),
                                                                     depth, Index(n), alpha, depth, depth, Index(0),
                                                                     Index(0));
            } else {
              internal::sme_gebp_neon<Scalar, false, false>(c.data(), rs, cs, pa.data(), pb.data(), Index(n), depth,
                                                            Index(n), alpha, depth, depth, Index(0), Index(0));
            }
            for (Index j = 0; j < n; ++j) {
              for (Index i = 0; i < n; ++i) {
                VERIFY((std::isfinite)(c(i, j).real()) && (std::isfinite)(c(i, j).imag()));
                VERIFY(c(i, j) == expected[t]);
              }
            }
          }
        }
      }
    }
  }
}

// With a real alpha, an overflowing imaginary part leaves the real part finite
// in every store of the NEON kernel: 5 - (2 + inf i) = 3 - inf i, where
// ar*re - ai*im would give NaN. Rows cover the packet and scalar-row tiles,
// rs = 2 the strided store.
template <typename Real>
static void test_neon_real_alpha_overflow() {
  using Scalar = std::complex<Real>;
  const Index PS = Index(internal::packet_traits<Real>::size);
  const Index MR = sme_mr<Scalar>();
  const Index cols = numext::mini(Index(5), Index(sme_nr<Scalar>()));
  const Index depth = 2;
  const Real big = Real(0.75) * NumTraits<Real>::highest();
  for (Index rows : {Index(1), PS, PS + 1, 2 * PS + 1}) {
    if (rows > MR) continue;
    std::vector<Scalar> pa(rows * depth), pb(cols * depth);
    for (Index k = 0; k < depth; ++k) {
      for (Index i = 0; i < rows; ++i) set_packed(pa.data(), rows, k, i, Scalar(1, big));
      for (Index j = 0; j < cols; ++j) set_packed(pb.data(), cols, k, j, Scalar(1, 0));
    }
    for (Index rs : {Index(1), Index(2)}) {
      const Index cs = rs * rows;
      std::vector<Scalar> c(cs * cols, Scalar(5, 0));
      internal::sme_gebp_neon<Scalar, false, false>(c.data(), rs, cs, pa.data(), pb.data(), rows, depth, cols,
                                                    Scalar(-1, 0), depth, depth, Index(0), Index(0));
      for (Index j = 0; j < cols; ++j) {
        for (Index i = 0; i < rows; ++i) {
          VERIFY_IS_EQUAL(c[j * cs + i * rs].real(), Real(3));
          VERIFY(c[j * cs + i * rs].imag() == -NumTraits<Real>::infinity());
        }
      }
    }
  }
}

// GCC 14 must keep the full-grid store's ZA tuples intact across C loads. Exact integer results and sentinel
// padding also cover the single-tile stores used for ragged rows and columns.
template <typename Scalar>
EIGEN_DONT_INLINE __arm_locally_streaming __arm_new("za") void sme_test_store_grid(Scalar* out, Index ldc,
                                                                                   const Scalar* a, const Scalar* b,
                                                                                   int rows, int cols) {
  using Traits = internal::sme_packet_traits<Scalar>;
  const int svl = Traits::size();
  const auto pg = Traits::ptrue();
  svzero_za();
  internal::outer_product_2x2<Scalar>(internal::ploadu(pg, a), internal::ploadu(pg, a + svl), internal::ploadu(pg, b),
                                      internal::ploadu(pg, b + svl));
  internal::sme_store_2x2_grid(out, Index(1), ldc, Scalar(2), Index(2), svl, rows - svl, Index(1),
                               internal::sme_min(cols, svl), cols > svl ? cols - svl : 0);
}

template <typename Scalar>
static void test_tile_pair_stores() {
  const int svl = internal::sme_packet_traits<Scalar>::size();
  const Index ldc = 2 * svl + 7;
  std::vector<Scalar> a(2 * svl), b(2 * svl);
  for (int i = 0; i < 2 * svl; ++i) {
    a[i] = Scalar(i % 7 + 1);
    b[i] = Scalar(i % 5 + 2);
  }
  const int widths[] = {1, 2, 3, 4, 5, svl - 1, svl, svl + 1, svl + 2, 2 * svl - 1, 2 * svl};
  std::vector<Scalar> out(ldc * (2 * svl + 3)), expected(out.size());
  for (int rows = svl + 1; rows <= 2 * svl; ++rows) {
    for (int cols : widths) {
      if (cols > 2 * svl) continue;
      std::fill(out.begin(), out.end(), Scalar(-9999));
      expected = out;
      for (int j = 0; j < cols; ++j)
        for (int i = 0; i < rows; ++i) expected[2 + i + (1 + j) * ldc] += Scalar(2) * a[i] * b[j];
      internal::sme_fpsr_guard guard;
      sme_test_store_grid(out.data(), ldc, a.data(), b.data(), rows, cols);
      for (Index i = 0; i < Index(out.size()); ++i) VERIFY_IS_EQUAL(out[i], expected[i]);
    }
  }
}

EIGEN_DECLARE_TEST(product_sme) {
  CALL_SUBTEST_1(test_whilelt_operand_types<float>());
  CALL_SUBTEST_1(test_packet_ops<float>());
  CALL_SUBTEST_1(test_products<float>());
  CALL_SUBTEST_1(test_conjugated_products<float>());
  CALL_SUBTEST_1(test_symm_pack<float>());
  CALL_SUBTEST_1(test_pack_direct<float>());
  CALL_SUBTEST_1(test_copy_panel_range<float>());
  CALL_SUBTEST_1(test_mapper_fallback<float>());
  CALL_SUBTEST_1(test_neon_small_blocks<float>());
  CALL_SUBTEST_1(test_small_products<float>());
  CALL_SUBTEST_1(test_direct_lhs<float>());
  CALL_SUBTEST_1(test_small_k_and_single_panel_rhs<float>());
  CALL_SUBTEST_1(test_deep_tail_panels<float>());
  CALL_SUBTEST_1(test_tiny_results<float>());
  CALL_SUBTEST_1(test_fp_flags<float>());
  CALL_SUBTEST_1(test_disjoint_parts<float>());

  // double reaches the SME kernel and packers only with FEAT_SME_F64F64; the
  // product sweep is meaningful either way, but the packed-layout tests name
  // specializations that only exist when it is available.
  CALL_SUBTEST_2(test_products<double>());
  CALL_SUBTEST_2(test_conjugated_products<double>());
#ifdef EIGEN_VECTORIZE_SME_F64F64
  CALL_SUBTEST_2(test_whilelt_operand_types<double>());
  CALL_SUBTEST_2(test_packet_ops<double>());
  CALL_SUBTEST_2(test_symm_pack<double>());
  CALL_SUBTEST_2(test_pack_direct<double>());
  CALL_SUBTEST_2(test_copy_panel_range<double>());
  CALL_SUBTEST_2(test_mapper_fallback<double>());
  CALL_SUBTEST_2(test_neon_small_blocks<double>());
  CALL_SUBTEST_2(test_small_products<double>());
  CALL_SUBTEST_2(test_direct_lhs<double>());
  CALL_SUBTEST_2(test_small_k_and_single_panel_rhs<double>());
  CALL_SUBTEST_2(test_deep_tail_panels<double>());
  CALL_SUBTEST_2(test_tiny_results<double>());
  CALL_SUBTEST_2(test_fp_flags<double>());
  CALL_SUBTEST_2(test_disjoint_parts<double>());
#endif

  CALL_SUBTEST_3(test_products<std::complex<float>>());
  CALL_SUBTEST_3(test_conjugated_products<std::complex<float>>());
  CALL_SUBTEST_3(test_symm_pack<std::complex<float>>());
  CALL_SUBTEST_3(test_pack_direct<std::complex<float>>());
  CALL_SUBTEST_3(test_copy_panel_range<std::complex<float>>());
  CALL_SUBTEST_3(test_mapper_fallback<std::complex<float>>());
  CALL_SUBTEST_3(test_neon_small_blocks<std::complex<float>>());
  CALL_SUBTEST_3(test_small_products<std::complex<float>>());
  CALL_SUBTEST_3(test_fp_flags<std::complex<float>>());
  CALL_SUBTEST_3(test_disjoint_parts<std::complex<float>>());
  CALL_SUBTEST_3((test_complex_alpha_range<float, ColMajor>()));
  CALL_SUBTEST_3((test_complex_alpha_range<float, RowMajor>()));
  CALL_SUBTEST_3(test_neon_real_alpha_overflow<float>());

  // complex<double> accumulates into ZA.D tiles, so it needs FEAT_SME_F64F64
  // exactly as double does.
  CALL_SUBTEST_4(test_products<std::complex<double>>());
  CALL_SUBTEST_4(test_conjugated_products<std::complex<double>>());
#ifdef EIGEN_VECTORIZE_SME_F64F64
  CALL_SUBTEST_4(test_symm_pack<std::complex<double>>());
  CALL_SUBTEST_4(test_pack_direct<std::complex<double>>());
  CALL_SUBTEST_4(test_copy_panel_range<std::complex<double>>());
  CALL_SUBTEST_4(test_mapper_fallback<std::complex<double>>());
  CALL_SUBTEST_4(test_neon_small_blocks<std::complex<double>>());
  CALL_SUBTEST_4(test_small_products<std::complex<double>>());
  CALL_SUBTEST_4(test_fp_flags<std::complex<double>>());
  CALL_SUBTEST_4(test_disjoint_parts<std::complex<double>>());
  CALL_SUBTEST_4((test_complex_alpha_range<double, ColMajor>()));
  CALL_SUBTEST_4((test_complex_alpha_range<double, RowMajor>()));
  CALL_SUBTEST_4(test_neon_real_alpha_overflow<double>());
#endif

  // A scalar type SME does not specialize, proving it still routes through the
  // generic product path inside an SME build -- where packet traits, alignment
  // and cache blocking all differ from a plain NEON build.
  CALL_SUBTEST_5(product(Matrix<long double, Dynamic, Dynamic>(33, 17)));

  CALL_SUBTEST_6(test_small_product_overflow<int>());
  CALL_SUBTEST_6(test_small_product_overflow<int64_t>());
  CALL_SUBTEST_6(test_small_products<float>());

#if defined(EIGEN_TEST_PART_7) || defined(EIGEN_TEST_PART_8)
  CALL_SUBTEST(test_small_products<float>());
  CALL_SUBTEST(test_small_products<std::complex<float>>());
#ifdef EIGEN_VECTORIZE_SME_F64F64
  CALL_SUBTEST(test_small_products<double>());
  CALL_SUBTEST(test_small_products<std::complex<double>>());
#endif
#endif
  CALL_SUBTEST_9(test_tile_pair_stores<float>());
#ifdef EIGEN_VECTORIZE_SME_F64F64
  CALL_SUBTEST_9(test_tile_pair_stores<double>());
#endif
}