SME: Add a double-precision GEMM kernel, and fix four streaming-mode pessimizations

libeigen/eigen!2873

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
This commit is contained in:
Rasmus Munk Larsen
2026-08-21 09:43:52 -07:00
co-authored by Rasmus Munk Larsen
parent 28d5232c8f
commit a474202651
9 changed files with 996 additions and 533 deletions
File diff suppressed because it is too large Load Diff
@@ -127,6 +127,19 @@ inline void manage_caching_sizes(Action action, std::ptrdiff_t* l1, std::ptrdiff
* \sa setCpuCacheSizes */
#ifdef EIGEN_VECTORIZE_SME
// True for the scalar pairs the SME gebp_kernel specializes (see
// arch/SME/GeneralBlockPanelKernel.h, which static_asserts that it agrees with
// this list); every other pair keeps Eigen's generic kernel, packers, cache
// blocking and GEMM loop order.
template <typename LhsScalar, typename RhsScalar>
struct sme_has_gebp_kernel : std::false_type {};
template <>
struct sme_has_gebp_kernel<float, float> : std::true_type {};
#ifdef EIGEN_VECTORIZE_SME_F64F64
template <>
struct sme_has_gebp_kernel<double, double> : std::true_type {};
#endif
// Overridable SME packed-panel budgets. The defaults are empirically tuned
// fp32 working-set limits for Apple M4 — heuristic budgets, not generic ARM64
// cache defaults; redefine them to retune for other SME implementations.
@@ -162,8 +175,11 @@ void evaluateProductBlockingSizesHeuristicForSme(Index& k, Index& m, Index& n) {
#endif
// Keep kc large enough to amortize SME setup and accumulation, but cap very
// deep products to avoid too many result store passes.
k = (numext::mini)(k, sme_max_kc);
// deep products to avoid too many result store passes. The cap is a scalar
// count tuned for fp32; scale it by the scalar width so every element type
// gets the same packed-panel byte budget.
const Index max_kc = (numext::maxi)(Index(1), sme_max_kc * Index(sizeof(float)) / Index(sizeof(LhsScalar)));
k = (numext::mini)(k, max_kc);
// Bound the packed RHS strip so very wide matrices do not allocate an
// unbounded blockB panel.
@@ -267,9 +283,9 @@ void evaluateProductBlockingSizesHeuristic(Index& k, Index& m, Index& n, Index n
if ((numext::maxi)(k, (numext::maxi)(m, n)) < 48) return;
#ifdef EIGEN_VECTORIZE_SME
// Only float×float uses the SME kernel; other scalar pairs run the generic
// kernel below and would thrash L1/L2 with the SME-sized budgets.
if (std::is_same<LhsScalar, float>::value && std::is_same<RhsScalar, float>::value) {
// Only the scalar pairs the SME kernel specializes use the SME budgets;
// the others run the generic kernel below and would thrash L1/L2 with them.
EIGEN_IF_CONSTEXPR ((sme_has_gebp_kernel<LhsScalar, RhsScalar>::value)) {
evaluateProductBlockingSizesHeuristicForSme<LhsScalar, RhsScalar>(k, m, n);
return;
}
@@ -228,11 +228,13 @@ struct general_matrix_matrix_product<Index, LhsScalar, LhsStorageOrder, Conjugat
ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, sizeA, blocking.blockA());
ei_declare_aligned_stack_constructed_variable(RhsScalar, blockB, sizeB, blocking.blockB());
// SME uses RHS-first order so consecutive gebp calls stream through
// adjacent row panels of a ColMajor result. Other kernels keep Eigen's
// The SME kernel uses RHS-first order so consecutive gebp calls stream
// through adjacent row panels of a ColMajor result. Other kernels --
// including the scalar pairs SME does not specialize -- keep Eigen's
// default LHS-first order.
#ifdef EIGEN_VECTORIZE_SME
using SequentialGemmLoop = gemm_pack_rhs_first_loop_policy;
using SequentialGemmLoop = std::conditional_t<sme_has_gebp_kernel<LhsScalar, RhsScalar>::value,
gemm_pack_rhs_first_loop_policy, gemm_pack_lhs_first_loop_policy>;
#else
using SequentialGemmLoop = gemm_pack_lhs_first_loop_policy;
#endif
@@ -462,6 +462,16 @@ extern "C" {
"EIGEN_ARM64_USE_SME must be built without -msve-vector-bits (scalable/VLA mode): a fixed SVE vector length pins the kernel to one runtime streaming SVL and silently miscomputes at any other."
#endif
// Double-precision outer products (FMOPA into a ZA.D tile) need the optional
// FEAT_SME_F64F64, which each compiler reports differently: GCC defines the ACLE
// macro, clang defines no macro but gates the builtin on the target feature.
// Both halves are needed -- clang otherwise accepts svmopa_za64_f64_m without the
// feature, so a missed gate faults at run time rather than at build time.
#if !defined(EIGEN_ARM64_NO_SME_F64F64) && \
(defined(__ARM_FEATURE_SME_F64F64) || EIGEN_HAS_BUILTIN(__builtin_sme_svmopa_za64_f64_m))
#define EIGEN_VECTORIZE_SME_F64F64
#endif
#elif EIGEN_ARCH_RISCV
#if defined(__riscv_zfh)
+35 -3
View File
@@ -429,9 +429,12 @@ build:linux:cross:arm:clang-14:default:
# line (ci/scripts/build.linux.script.sh), so each -D value must be a single
# shell word: multi-token flags use CMake's ';' list separator, with no
# spaces or embedded quotes (cf. the arm / generic-vector jobs).
# +sme-f64f64 (FEAT_SME_F64F64) enables the double-precision GEMM kernel:
# GCC refuses to compile its fp64 FMOPA without it, and qemu's `max` CPU
# implements the extension, so the test jobs below exercise that kernel.
EIGEN_CI_ADDITIONAL_ARGS: >-
-DEIGEN_TEST_SME=ON
-DEIGEN_TEST_CUSTOM_CXX_FLAGS=-march=armv9.2-a+sme2
-DEIGEN_TEST_CUSTOM_CXX_FLAGS=-march=armv9.2-a+sme2+sme-f64f64
-DCMAKE_CROSSCOMPILING_EMULATOR=qemu-aarch64;-L;/usr/aarch64-linux-gnu
# SME ACLE (arm_sme.h, __arm_streaming, +sme2) requires GCC >= 14.
# The build stage only cross-compiles and never invokes the emulator, so
@@ -440,8 +443,9 @@ build:linux:cross:arm:clang-14:default:
# product_symm/product_trsolve/product_trmm/product_syrk/mixingtypes exercise
# the SME packers through their non-dense-GEMM consumers (selfadjoint,
# triangular-solve, triangular-product, rank-update, mixed-scalar) that
# product_sme can't reach. Explicit list, not `buildtests`: cross-building
# every test under the SME flags overruns the runner timeout.
# product_sme can't reach; each has double subtests, so they cover both SME
# scalar types. Explicit list, not `buildtests`: cross-building every test
# under the SME flags overruns the runner timeout.
EIGEN_CI_BUILD_TARGET: product_sme product_symm product_trsolve product_trmm product_syrk mixingtypes
# Inherit the 2xlarge runner from .build:linux:cross: GitLab's open-source
# program bills small and large runners the same, so don't downsize.
@@ -475,6 +479,28 @@ build:linux:cross:sme:gcc-14:full:
EIGEN_CI_BUILD_TARGET: buildtests
timeout: 2h
# FEAT_SME_F64F64 is optional, and every job above builds with it, so nothing
# else compiles the configuration where it is absent: the #else static_assert
# that keeps double off the GEMM driver, and double on the generic kernel inside
# an SME build. GCC hard-errors on the fp64 FMOPA without the feature, so this
# is also what catches an fp64 SME intrinsic escaping the
# EIGEN_VECTORIZE_SME_F64F64 gate. The EIGEN_ARM64_NO_SME_F64F64 opt-out reaches
# the same branch, so one job covers both ways in. product_sme alone: its double
# subtest is written to run in either configuration, and one target keeps the
# job cheap.
build:linux:cross:sme:gcc-14:no-f64f64:
extends: .build:linux:cross:sme
variables:
EIGEN_CI_C_COMPILER: gcc-14
EIGEN_CI_CXX_COMPILER: g++-14
EIGEN_CI_CROSS_C_COMPILER: aarch64-linux-gnu-gcc-14
EIGEN_CI_CROSS_CXX_COMPILER: aarch64-linux-gnu-g++-14
# The anchor's arguments, less +sme-f64f64.
EIGEN_CI_ADDITIONAL_ARGS: >-
-DEIGEN_TEST_SME=ON
-DEIGEN_TEST_CUSTOM_CXX_FLAGS=-march=armv9.2-a+sme2
-DCMAKE_CROSSCOMPILING_EMULATOR=qemu-aarch64;-L;/usr/aarch64-linux-gnu
EIGEN_CI_BUILD_TARGET: product_sme
build:linux:cross:aarch64:gcc-10:default:
@@ -728,3 +754,9 @@ build:linux:riscv64:gcc-15:default:affected:
build:linux:cross:sme:gcc-14:full:affected:
extends: build:linux:cross:sme:gcc-14:full
rules: !reference [.rules:libeigen:affected-tests:sve-sme, rules]
# One product_sme build, so a change to the gate is caught by the affected tier
# rather than only by the nightly run.
build:linux:cross:sme:gcc-14:no-f64f64:affected:
extends: build:linux:cross:sme:gcc-14:no-f64f64
rules: !reference [.rules:libeigen:affected-tests:sve-sme, rules]
+12
View File
@@ -470,6 +470,18 @@ test:linux:sme2048:gcc-14:default:official:
variables:
QEMU_CPU: max,sme=on,sme2048=on,sme128=off,sme256=off,sme512=off,sme1024=off
# The no-FEAT_SME_F64F64 build, at one SVL: double has to reach the generic
# kernel and still be correct inside an SME build, where packet traits,
# alignment and cache blocking all differ from a plain NEON build. Its subtest
# skips the SME-only symm and mapper cases there, so the run is short. One SVL
# is enough -- the generic kernel it exercises is SVL-independent.
test:linux:sme512:gcc-14:no-f64f64:official:
extends: .test:linux:sme:gcc-14:default
needs: [ build:linux:cross:sme:gcc-14:no-f64f64 ]
variables:
QEMU_CPU: max,sme=on,sme512=on,sme128=off,sme256=off,sme1024=off,sme2048=off
EIGEN_CI_CTEST_REGEX: product_sme
# Full nightly SME test run: the complete Official + Unsupported ctest suite
# against the full build above, at a single representative SVL (512: the logical
# block matches the ZA tile grid exactly, exercising the multi-vector-load fast
+6 -2
View File
@@ -124,8 +124,12 @@ run time. However, these assertions do cost time and can thus be turned off.
- \b \c EIGEN_DEFAULT_L1_CACHE_SIZE - Sets the default L1 cache size that is used in Eigen's GEBP kernel when the correct cache size cannot be determined at runtime.
- \b \c EIGEN_DEFAULT_L2_CACHE_SIZE - Sets the default L2 cache size that is used in Eigen's GEBP kernel when the correct cache size cannot be determined at runtime.
- \b \c EIGEN_DEFAULT_L3_CACHE_SIZE - Sets the default L3 cache size that is used in Eigen's GEBP kernel when the correct cache size cannot be determined at runtime.
- \b \c EIGEN_SME_MAX_KC - Maximum depth (k) blocking size used by the ARM SME GEMM kernel. The default (2048) is
empirically tuned for Apple M4; override to retune for other SME implementations.
- \b \c EIGEN_SME_MAX_KC - Maximum depth (k) blocking size used by the ARM SME GEMM kernel, expressed in \c float
elements and scaled by the scalar width for wider types. The default (2048) is empirically tuned for Apple M4;
override to retune for other SME implementations.
- \b \c EIGEN_ARM64_NO_SME_F64F64 - Disables the double-precision ARM SME GEMM kernel, which needs the optional
FEAT_SME_F64F64 extension. %Eigen enables it whenever the compiler reports that extension; define this macro when
the build target is wider than the run target, so that \c double keeps the generic kernel.
- \b \c EIGEN_SME_PACKED_RHS_BUDGET_BYTES - Byte budget for the packed RHS panel in the ARM SME GEMM kernel.
The default (32 MB) is empirically tuned for Apple M4.
- \b \c EIGEN_SME_LHS_WORKING_SET_BUDGET_BYTES - Byte budget for the LHS working set in the ARM SME GEMM kernel.
+5 -2
View File
@@ -20,8 +20,11 @@ On CPUs, %Eigen provides vectorized kernels for the following instruction sets:
<tr class="alt"><td>ARM / AArch64</td><td>NEON; SVE and SME as opt-in backends</td>
<td>SVE requires \c EIGEN_ARM64_USE_SVE and a fixed vector length
(<tt>-msve-vector-bits=N</tt>); the SME backend, enabled with \c EIGEN_ARM64_USE_SME,
accelerates matrix products and must be built \em without <tt>-msve-vector-bits</tt>, since
a fixed length would pin the kernels to one runtime streaming vector length</td></tr>
accelerates \c float and \c double matrix products and must be built \em without
<tt>-msve-vector-bits</tt>, since a fixed length would pin the kernels to one runtime
streaming vector length. Double precision additionally needs the optional FEAT_SME_F64F64
extension (<tt>+sme-f64f64</tt>, or a <tt>-mcpu</tt> that implies it); without it \c double
keeps the generic kernel</td></tr>
<tr><td>PowerPC</td><td>AltiVec, VSX, MMA</td><td></td></tr>
<tr class="alt"><td>IBM Z (s390x)</td><td>ZVector</td><td></td></tr>
<tr><td>MIPS</td><td>MSA</td><td></td></tr>
+229 -180
View File
@@ -9,6 +9,10 @@
// SME GEMM kernel tests.
// Requires compiler flags: -march=armv9.2-a+sme2 and -DEIGEN_ARM64_USE_SME.
// Double precision additionally needs FEAT_SME_F64F64 (+sme-f64f64, or a -mcpu
// that implies it); without it EIGEN_VECTORIZE_SME_F64F64 is undefined and
// double keeps the generic kernel, so the double subtest packs its cases
// through that path instead.
#include "product.h"
@@ -23,10 +27,28 @@
"the typical CMake invocation)."
#endif
using SmeColMajorMatF = Matrix<float, Dynamic, Dynamic, ColMajor>;
using SmeRowMajorMatF = Matrix<float, Dynamic, Dynamic, RowMajor>;
using SmeColMajorStridedMatF = Map<SmeColMajorMatF, 0, Stride<Dynamic, Dynamic>>;
using SmeRowMajorStridedMatF = Map<SmeRowMajorMatF, 0, Stride<Dynamic, Dynamic>>;
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 width for Scalar: kSmeMr for float, kSmeMrD
// for double. Sizes below are expressed in terms of it so each scalar sweeps
// its own block and ZA-tile boundaries.
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 InputMat, typename ResultMat, typename ResultMap>
static void verify_strided_result(int n, ResultMat& storage, const Stride<Dynamic, Dynamic>& stride) {
@@ -43,44 +65,48 @@ static void verify_strided_result(int n, ResultMat& storage, const Stride<Dynami
VERIFY_IS_APPROX(got, ref);
}
template <typename InputMat>
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.
SmeColMajorMatF storage = SmeColMajorMatF::Zero(2 * n, n);
verify_strided_result<InputMat, SmeColMajorMatF, SmeColMajorStridedMatF>(
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) == float(0));
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.
SmeRowMajorMatF storage = SmeRowMajorMatF::Zero(n, 2 * n);
verify_strided_result<SmeRowMajorMatF, SmeRowMajorMatF, SmeRowMajorStridedMatF>(
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) == float(0));
VERIFY(storage(i, j) == Scalar(0));
}
}
}
// 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;
constexpr int depth = 2050;
const int depth = int(2 * (Index(EIGEN_SME_MAX_KC) * Index(sizeof(float)) / Index(sizeof(Scalar)))) + 2;
constexpr int cols = 64;
SmeColMajorMatF A = SmeColMajorMatF::Random(rows, depth);
SmeColMajorMatF B = SmeColMajorMatF::Random(depth, cols);
SmeColMajorMatF C = SmeColMajorMatF::Random(rows, cols);
SmeColMajorMatF c_before = C;
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;
@@ -102,81 +128,86 @@ static void test_deep_k_split() {
// 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.
static constexpr float kPackSentinel = 98765.0f;
template <typename Scalar>
static Scalar pack_sentinel() {
return Scalar(98765);
}
// Lower-triangular n x n operand plus the dense symmetric reference the packer
// must emit. The unused triangle is filled with kPackSentinel so a packer that
// 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 stored(col,row) above it.
template <int StorageOrder>
static void make_lower_stored_symmetric(Index n, Matrix<float, Dynamic, Dynamic, StorageOrder>& stored,
Matrix<float, Dynamic, Dynamic, StorageOrder>& full) {
full = Matrix<float, Dynamic, Dynamic, StorageOrder>::Random(n, n);
full = ((full + full.transpose()) * 0.5f).eval();
stored = Matrix<float, Dynamic, Dynamic, StorageOrder>::Constant(n, n, kPackSentinel);
template <typename Scalar, int StorageOrder>
static void make_lower_stored_symmetric(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.transpose()) * 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 kSmeMr-wide depth-major panels. Reference: full(i+r, k).
template <int StorageOrder>
// 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 = internal::kSmeMr;
Matrix<float, Dynamic, Dynamic, StorageOrder> stored, full;
make_lower_stored_symmetric<StorageOrder>(n, stored, full);
const Index MR = sme_mr<Scalar>();
Matrix<Scalar, Dynamic, Dynamic, StorageOrder> stored, full;
make_lower_stored_symmetric<Scalar, StorageOrder>(n, stored, full);
VectorXf packed = VectorXf::Constant(n * n, kPackSentinel);
VectorXf ref = VectorXf::Constant(n * n, kPackSentinel);
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) ref[i * n + k * w + r] = full(i + r, k);
}
internal::symm_pack_lhs<float, Index, internal::kSmeMr, 1, StorageOrder> pack;
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 kSmeNr-wide depth-major panels. Reference:
// 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 <int StorageOrder>
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 = internal::kSmeNr;
Matrix<float, Dynamic, Dynamic, StorageOrder> stored, full;
make_lower_stored_symmetric<StorageOrder>(N, stored, full);
const Index NR = sme_nr<Scalar>();
Matrix<Scalar, Dynamic, Dynamic, StorageOrder> stored, full;
make_lower_stored_symmetric<Scalar, StorageOrder>(N, stored, full);
VectorXf packed = VectorXf::Constant(cols * depth, kPackSentinel);
VectorXf ref = VectorXf::Constant(cols * depth, kPackSentinel);
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) ref[j * depth + k * w + c] = full(k2 + k, j + c);
}
internal::symm_pack_rhs<float, Index, internal::kSmeNr, StorageOrder> pack;
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 (svlw, 2*svlw). The spread below
// hits that for svlw in {4, 8, 16, 32, 64} (SVL 128..2048).
// 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<ColMajor>(n);
verify_symm_pack_lhs<RowMajor>(n);
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<ColMajor>(n, n, n, 0);
verify_symm_pack_rhs<RowMajor>(n, n, n, 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
@@ -189,8 +220,8 @@ static void test_symm_pack() {
{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<ColMajor>(c.N, c.depth, c.cols, c.k2);
verify_symm_pack_rhs<RowMajor>(c.N, c.depth, c.cols, c.k2);
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);
}
}
@@ -207,53 +238,57 @@ static void test_symm_pack() {
// 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 float* data;
const Scalar* data;
Index stride; // element(i, k) = data[i + k * stride], contiguous in i
float operator()(Index i, Index k) const { return data[i + k * stride]; }
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 float* data;
const Scalar* data;
Index stride; // element(i, k) = data[i * stride + k], contiguous in k
float operator()(Index i, Index k) const { return data[i * stride + 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 float* p; // &element(0, col); contiguous in depth
float operator()(Index k) const { return p[k]; }
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 float* data;
const Scalar* data;
Index stride; // element(k, col) = data[k + col * stride], contiguous in k
float operator()(Index k, Index col) const { return data[k + col * stride]; }
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 float* p; // &element(row, col); packet offsets advance columns
float operator()(Index offset) const { return p[offset]; }
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 float* data;
const Scalar* data;
Index stride; // element(k, col) = data[k * stride + col], contiguous in col
float operator()(Index k, Index col) const { return data[k * stride + 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}; }
};
@@ -261,7 +296,7 @@ struct ByValueRowMajorRhsMapper {
// 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 kSmeMr or kSmeNr.
// `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) {
@@ -271,11 +306,11 @@ static Index packed_len(Index outer, Index depth, Index unit, Index dst_stride,
return end;
}
template <bool PanelMode, typename MatrixType>
static void fill_lhs_ref(VectorXf& ref, const MatrixType& V, Index rows, Index depth, Index dst_stride,
template <typename Scalar, bool PanelMode, 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 = internal::kSmeMr;
ref.setConstant(kPackSentinel);
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;
@@ -284,11 +319,11 @@ static void fill_lhs_ref(VectorXf& ref, const MatrixType& V, Index rows, Index d
}
}
template <bool PanelMode, typename MatrixType>
static void fill_rhs_ref(VectorXf& ref, const MatrixType& V, Index cols, Index depth, Index dst_stride,
template <typename Scalar, bool PanelMode, 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 = internal::kSmeNr;
ref.setConstant(kPackSentinel);
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;
@@ -300,12 +335,12 @@ static void fill_rhs_ref(VectorXf& ref, const MatrixType& V, Index cols, Index d
// 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 <int StorageOrder, bool PanelMode>
template <typename Scalar, int StorageOrder, bool PanelMode>
static void verify_fallback_lhs_strided(Index rows, Index depth, Index incr) {
typedef internal::blas_data_mapper<float, Index, StorageOrder, Unaligned, Dynamic> Mapper;
MatrixXf V = MatrixXf::Random(rows, depth);
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;
VectorXf buf = VectorXf::Zero((StorageOrder == ColMajor ? depth : rows) * mstride + 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);
@@ -313,12 +348,12 @@ static void verify_fallback_lhs_strided(Index rows, Index depth, Index incr) {
const Index dst_stride = PanelMode ? depth + 5 : 0;
const Index dst_offset = PanelMode ? 3 : 0;
const Index len = packed_len<PanelMode>(rows, depth, internal::kSmeMr, dst_stride, dst_offset);
VectorXf packed = VectorXf::Constant(len, kPackSentinel);
VectorXf ref(len);
fill_lhs_ref<PanelMode>(ref, V, rows, depth, dst_stride, dst_offset);
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>(ref, V, rows, depth, dst_stride, dst_offset);
internal::gemm_pack_lhs<float, Index, Mapper, internal::kSmeMr, 1, typename internal::packet_traits<float>::type,
internal::gemm_pack_lhs<Scalar, Index, Mapper, sme_mr<Scalar>(), 1, typename internal::packet_traits<Scalar>::type,
StorageOrder, false, PanelMode>
pack;
pack(packed.data(), mapper, depth, rows, dst_stride, dst_offset);
@@ -327,22 +362,22 @@ static void verify_fallback_lhs_strided(Index rows, Index depth, Index incr) {
// By-value LHS mappers exercise both packet directions. RowMajor must stay
// scalar because its packets advance depth rather than rows.
template <int StorageOrder, bool PanelMode>
template <typename Scalar, int StorageOrder, bool PanelMode>
static void verify_fallback_lhs_byvalue(Index rows, Index depth) {
typedef Matrix<float, Dynamic, Dynamic, StorageOrder> MatrixType;
typedef typename std::conditional<StorageOrder == ColMajor, ByValueColMajorLhsMapper, ByValueRowMajorLhsMapper>::type
Mapper;
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);
Mapper mapper{V.data(), V.outerStride()};
const Index dst_stride = PanelMode ? depth + 5 : 0;
const Index dst_offset = PanelMode ? 3 : 0;
const Index len = packed_len<PanelMode>(rows, depth, internal::kSmeMr, dst_stride, dst_offset);
VectorXf packed = VectorXf::Constant(len, kPackSentinel);
VectorXf ref(len);
fill_lhs_ref<PanelMode>(ref, V, rows, depth, dst_stride, dst_offset);
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>(ref, V, rows, depth, dst_stride, dst_offset);
internal::gemm_pack_lhs<float, Index, Mapper, internal::kSmeMr, 1, typename internal::packet_traits<float>::type,
internal::gemm_pack_lhs<Scalar, Index, Mapper, sme_mr<Scalar>(), 1, typename internal::packet_traits<Scalar>::type,
StorageOrder, false, PanelMode>
pack;
pack(packed.data(), mapper, depth, rows, dst_stride, dst_offset);
@@ -351,12 +386,12 @@ static void verify_fallback_lhs_byvalue(Index rows, Index 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 <int StorageOrder, bool PanelMode>
template <typename Scalar, int StorageOrder, bool PanelMode>
static void verify_fallback_rhs_strided(Index depth, Index cols, Index incr) {
typedef internal::blas_data_mapper<float, Index, StorageOrder, Unaligned, Dynamic> Mapper;
MatrixXf V = MatrixXf::Random(depth, cols);
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;
VectorXf buf = VectorXf::Zero((StorageOrder == ColMajor ? cols : depth) * mstride + 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);
@@ -364,109 +399,106 @@ static void verify_fallback_rhs_strided(Index depth, Index cols, Index incr) {
const Index dst_stride = PanelMode ? depth + 5 : 0;
const Index dst_offset = PanelMode ? 3 : 0;
const Index len = packed_len<PanelMode>(cols, depth, internal::kSmeNr, dst_stride, dst_offset);
VectorXf packed = VectorXf::Constant(len, kPackSentinel);
VectorXf ref(len);
fill_rhs_ref<PanelMode>(ref, V, cols, depth, dst_stride, dst_offset);
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>(ref, V, cols, depth, dst_stride, dst_offset);
internal::gemm_pack_rhs<float, Index, Mapper, internal::kSmeNr, StorageOrder, false, PanelMode> pack;
internal::gemm_pack_rhs<Scalar, Index, Mapper, sme_nr<Scalar>(), StorageOrder, false, PanelMode> pack;
pack(packed.data(), mapper, depth, cols, dst_stride, dst_offset);
VERIFY_IS_EQUAL(packed, ref);
}
// By-value RHS mappers likewise cover both packet directions. RowMajor packets
// advance columns, so the depth-oriented transpose fallback must stay scalar.
template <int StorageOrder, bool PanelMode>
template <typename Scalar, int StorageOrder, bool PanelMode>
static void verify_fallback_rhs_byvalue(Index depth, Index cols) {
typedef Matrix<float, Dynamic, Dynamic, StorageOrder> MatrixType;
typedef typename std::conditional<StorageOrder == ColMajor, ByValueColMajorRhsMapper, ByValueRowMajorRhsMapper>::type
Mapper;
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);
Mapper mapper{V.data(), V.outerStride()};
const Index dst_stride = PanelMode ? depth + 5 : 0;
const Index dst_offset = PanelMode ? 3 : 0;
const Index len = packed_len<PanelMode>(cols, depth, internal::kSmeNr, dst_stride, dst_offset);
VectorXf packed = VectorXf::Constant(len, kPackSentinel);
VectorXf ref(len);
fill_rhs_ref<PanelMode>(ref, V, cols, depth, dst_stride, dst_offset);
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>(ref, V, cols, depth, dst_stride, dst_offset);
internal::gemm_pack_rhs<float, Index, Mapper, internal::kSmeNr, StorageOrder, false, PanelMode> pack;
internal::gemm_pack_rhs<Scalar, Index, Mapper, sme_nr<Scalar>(), StorageOrder, false, PanelMode> pack;
pack(packed.data(), mapper, depth, cols, dst_stride, dst_offset);
VERIFY_IS_EQUAL(packed, ref);
}
template <typename Scalar>
static void test_mapper_fallback() {
const int widths[] = {4, 5, 32, 33, 65}; // rows/cols around 4 and 32
const int depths[] = {1, 3, 8, 35}; // depth remainders 1..3 and larger
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) {
for (int incr : {2, 3}) {
verify_fallback_lhs_strided<ColMajor, false>(n, d, incr);
verify_fallback_lhs_strided<ColMajor, true>(n, d, incr);
verify_fallback_lhs_strided<RowMajor, false>(n, d, incr); // scalar path
verify_fallback_lhs_strided<RowMajor, true>(n, d, incr);
verify_fallback_rhs_strided<ColMajor, false>(d, n, incr);
verify_fallback_rhs_strided<ColMajor, true>(d, n, incr);
verify_fallback_rhs_strided<RowMajor, false>(d, n, incr); // scalar path
verify_fallback_rhs_strided<RowMajor, true>(d, n, incr);
verify_fallback_lhs_strided<Scalar, ColMajor, false>(n, d, incr);
verify_fallback_lhs_strided<Scalar, ColMajor, true>(n, d, incr);
verify_fallback_lhs_strided<Scalar, RowMajor, false>(n, d, incr); // scalar path
verify_fallback_lhs_strided<Scalar, RowMajor, true>(n, d, incr);
verify_fallback_rhs_strided<Scalar, ColMajor, false>(d, n, incr);
verify_fallback_rhs_strided<Scalar, ColMajor, true>(d, n, incr);
verify_fallback_rhs_strided<Scalar, RowMajor, false>(d, n, incr); // scalar path
verify_fallback_rhs_strided<Scalar, RowMajor, true>(d, n, incr);
}
verify_fallback_lhs_byvalue<ColMajor, false>(n, d);
verify_fallback_lhs_byvalue<ColMajor, true>(n, d);
verify_fallback_lhs_byvalue<RowMajor, false>(n, d);
verify_fallback_lhs_byvalue<RowMajor, true>(n, d);
verify_fallback_rhs_byvalue<ColMajor, false>(d, n);
verify_fallback_rhs_byvalue<ColMajor, true>(d, n);
verify_fallback_rhs_byvalue<RowMajor, false>(d, n);
verify_fallback_rhs_byvalue<RowMajor, true>(d, n);
verify_fallback_lhs_byvalue<Scalar, ColMajor, false>(n, d);
verify_fallback_lhs_byvalue<Scalar, ColMajor, true>(n, d);
verify_fallback_lhs_byvalue<Scalar, RowMajor, false>(n, d);
verify_fallback_lhs_byvalue<Scalar, RowMajor, true>(n, d);
verify_fallback_rhs_byvalue<Scalar, ColMajor, false>(d, n);
verify_fallback_rhs_byvalue<Scalar, ColMajor, true>(d, n);
verify_fallback_rhs_byvalue<Scalar, RowMajor, false>(d, n);
verify_fallback_rhs_byvalue<Scalar, RowMajor, true>(d, n);
}
}
}
EIGEN_DECLARE_TEST(product_sme) {
// Square edge cases around the block and tile boundaries (the block is
// kSmeMr x kSmeNr and a ZA tile is svlw x svlw, so the sizes below land
// just on/off the intra-tile splits and the block tails at SVL=512).
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(1, 1)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(15, 15)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(16, 16)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(17, 17)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(31, 31)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(33, 33)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(63, 63)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(64, 64)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(65, 65)));
// ---------------------------------------------------------------------------
// Product-level coverage, swept relative to the scalar's own block width.
// ---------------------------------------------------------------------------
// Thin / wide rectangular cases (M x 1, 1 x N)
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(32, 1)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(1, 32)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(1, 64)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(64, 1)));
// Sizes that land just on and off the block tails and the intra-block ZA-tile
// splits. MR/2 is the tile side at the SVL=512 design point.
template <typename Scalar>
static std::vector<int> sme_edge_sizes() {
const int MR = sme_mr<Scalar>();
return {1, MR / 2 - 1, MR / 2, MR / 2 + 1, MR - 1, MR, MR + 1, 2 * MR - 1, 2 * MR, 2 * MR + 1};
}
// Non-float scalar smoke tests: SME only specializes fp32, so these prove
// unsupported scalar types still route through the generic product path.
CALL_SUBTEST_2(product(Matrix<double, Dynamic, Dynamic>(33, 17)));
CALL_SUBTEST_3(product(Matrix<std::complex<float>, Dynamic, Dynamic>(33, 17)));
template <typename Scalar>
static void test_products() {
const int MR = sme_mr<Scalar>();
// Non-square cases that exercise tail paths for both M and N
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(17, 65)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(65, 17)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(15, 63)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(33, 7)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(7, 33)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(128, 3)));
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(3, 128)));
// Square edge cases around the block and tile boundaries.
for (int n : sme_edge_sizes<Scalar>()) product(SmeColMajorMat<Scalar>(n, n));
// Exercise the kc split path just above the SME blocking heuristic's depth
// cap (sme_max_kc in GeneralBlockPanelKernel.h).
test_deep_k_split();
test_symm_pack();
test_mapper_fallback();
// 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>();
// Random sizes
for (int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1(product(Matrix<float, Dynamic, Dynamic>(internal::random<int>(1, EIGEN_TEST_MAX_SIZE),
internal::random<int>(1, EIGEN_TEST_MAX_SIZE))));
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
@@ -476,40 +508,57 @@ EIGEN_DECLARE_TEST(product_sme) {
// - 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).
CALL_SUBTEST_1(product(SmeRowMajorMatF(15, 15)));
CALL_SUBTEST_1(product(SmeRowMajorMatF(16, 16)));
CALL_SUBTEST_1(product(SmeRowMajorMatF(17, 17)));
CALL_SUBTEST_1(product(SmeRowMajorMatF(31, 31)));
CALL_SUBTEST_1(product(SmeRowMajorMatF(32, 32)));
CALL_SUBTEST_1(product(SmeRowMajorMatF(33, 33)));
CALL_SUBTEST_1(product(SmeRowMajorMatF(64, 64)));
CALL_SUBTEST_1(product(SmeRowMajorMatF(65, 65)));
CALL_SUBTEST_1(product(SmeRowMajorMatF(17, 65)));
CALL_SUBTEST_1(product(SmeRowMajorMatF(65, 17)));
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++) {
CALL_SUBTEST_1(product(
SmeRowMajorMatF(internal::random<int>(1, EIGEN_TEST_MAX_SIZE), internal::random<int>(1, EIGEN_TEST_MAX_SIZE))));
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 : {15, 16, 17, 31, 32, 33, 63, 64, 65}) {
test_general_strided_result<SmeColMajorMatF>(n);
test_general_strided_result<SmeRowMajorMatF>(n);
test_rowmajor_strided_result(n);
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 : {15, 16, 17, 31, 32, 33, 63, 64, 65}) {
Matrix<float, Dynamic, Dynamic, RowMajor> A = Matrix<float, Dynamic, Dynamic, RowMajor>::Random(n, n);
Matrix<float, Dynamic, Dynamic, RowMajor> B = Matrix<float, Dynamic, Dynamic, RowMajor>::Random(n, n);
SmeColMajorMatF C = SmeColMajorMatF::Zero(n, n);
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());
}
}
EIGEN_DECLARE_TEST(product_sme) {
CALL_SUBTEST_1(test_products<float>());
CALL_SUBTEST_1(test_symm_pack<float>());
CALL_SUBTEST_1(test_mapper_fallback<float>());
// double only reaches the SME kernel and packers 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>());
#ifdef EIGEN_VECTORIZE_SME_F64F64
CALL_SUBTEST_2(test_symm_pack<double>());
CALL_SUBTEST_2(test_mapper_fallback<double>());
#endif
// Scalar types SME does not specialize: these prove they still route through
// the generic product path.
CALL_SUBTEST_3(product(Matrix<std::complex<float>, Dynamic, Dynamic>(33, 17)));
CALL_SUBTEST_3(product(Matrix<std::complex<double>, Dynamic, Dynamic>(33, 17)));
}