arm64 SME2: faster SGEMM/DGEMM, and SYMM/SYRK/SYR2K/TRMM/TRSM on the SME GEMM kernel

CONTRIBUTORS.md

@@ -288,3 +288,6 @@ hheei <hheei@users.noreply.github.com>
 
 * Hugo Meiland <hugo@meiland.nl>
   * [2026-08-09] Add Cortex-A72 DGEMM 6x8 microkernel and blocking
+
+* Michael Tesch <https://github.com/tesch1>
+  * [2026-09-29] SME2 SGEMM/DGEMM kernel for Apple M4; SYMM, SYRK, SYR2K, TRMM and TRSM on the SME GEMM kernel

README.md

@@ -187,7 +187,7 @@ Please read `GotoBLAS_01Readme.txt` for older CPU models already supported by th
 - **Neoverse V3**: preliminary support
 - **Neoverse V3AE**: preliminary support
 - **Apple Vortex**: preliminary support based on ThunderX2/3
-- **Apple VortexM4**: preliminary support based on ThunderX2/3, SME kernels for SGEMM,SSYMM,STRMM,SSYRK,SSYR2K
+- **Apple VortexM4**: preliminary support based on ThunderX2/3, SME kernels for SGEMM,SSYMM,STRMM,SSYRK,SSYR2K; on SME2 cores with a 512-bit streaming vector length (tested on M4) an SME2 kernel for SGEMM and DGEMM, which S/D SYMM, SYRK, SYR2K, TRMM and TRSM also use
 - **A64FX**:  preliminary support, optimized Level-3 BLAS
 - **ARMV8SVE**: any ARMV8 cpu with SVE extensions 
 - **ARMV9SME**: any ARMV9 cpu with SVE and SME extensions 

cpuid_arm64.c

@@ -408,7 +408,13 @@ int detect(void)
 	if (value64 == 3660830781) return CPU_VORTEX; //A15/M2
         if (value64 == 2271604202) return CPU_VORTEX; //A16/M3
         if (value64 == 1867590060) return CPU_VORTEXM4; //M4
+        if (value64 == 399882554) return CPU_VORTEXM4; //M4 Pro/Max
 	    if (value64 == 492472296) return CPU_VORTEXM4; //M5
+	{ /* later Apple cores with SME: same kernels until they get their own entry */
+	int sme = 0;
+	size_t len = sizeof(sme);
+	if (sysctlbyname("hw.optional.arm.FEAT_SME", &sme, &len, NULL, 0) == 0 && sme) return CPU_VORTEXM4;
+	}
 #else
 #ifdef OS_WINDOWS
 	HKEY reghandle;

interface/gemm.c

@@ -586,6 +586,12 @@ if (strcmp(gotoblas_corename(), "armv9sme") == 0
 #endif
 )
 #endif //defined dynarch
+/* below this the NEON paths win over the cost of entering streaming mode (Apple M4: about 20^3, 16^3 in fp64) */
+#if defined(DOUBLE) && !defined(COMPLEX)
+if ((double)args.m * (double)args.n * (double)args.k >= 4000.)
+#else
+if ((double)args.m * (double)args.n * (double)args.k >= 8000.)
+#endif
 {
 char* TA,*TB;
       if (transa & 1) 

interface/symm.c

@@ -121,6 +121,10 @@ extern char* gotoblas_corename(void);
 #endif
 #endif
 
+#if defined(ARCH_ARM64)
+#include "../kernel/arm64/sme_level3.h"
+#endif
+
 static int (*symm[])(blas_arg_t *, BLASLONG *, BLASLONG *, FLOAT *, FLOAT *, BLASLONG) = {
 #ifndef GEMM3M
 #ifndef HEMM
@@ -398,6 +402,10 @@ if (strcmp(gotoblas_corename(), "armv9sme") == 0
 #endif
 #endif
 
+#endif
+
+#ifdef SME_LEVEL3 /* arm64 SME: recursive blocking on the SME GEMM kernel */
+  if (s3_symm_hook(side, uplo, &args)) return;
 #endif
 
   IDEBUG_START;

interface/syr2k.c

@@ -71,6 +71,10 @@
 #endif
 #endif
 
+#if defined(ARCH_ARM64)
+#include "../kernel/arm64/sme_level3.h"
+#endif
+
 static int (*syr2k[])(blas_arg_t *, BLASLONG *, BLASLONG *, FLOAT *, FLOAT *, BLASLONG) = {
 #ifndef HEMM
   SYR2K_UN, SYR2K_UC, SYR2K_LN, SYR2K_LC,
@@ -388,6 +392,10 @@ if (strcmp(gotoblas_corename(), "armv9sme") == 0
 #endif
 
 
+#ifdef SME_LEVEL3 /* arm64 SME: recursive blocking on the SME GEMM kernel */
+  if (s3_syrk_hook(1, uplo, trans, &args)) return;
+#endif
+
   IDEBUG_START;
 
   FUNCTION_PROFILE_START();

interface/syrk.c

@@ -77,6 +77,10 @@
 #define GEMM_MULTITHREAD_THRESHOLD 4
 #endif
 
+#if defined(ARCH_ARM64)
+#include "../kernel/arm64/sme_level3.h"
+#endif
+
 static int (*syrk[])(blas_arg_t *, BLASLONG *, BLASLONG *, FLOAT *, FLOAT *, BLASLONG) = {
 #ifndef HEMM
   SYRK_UN, SYRK_UC, SYRK_LN, SYRK_LC,
@@ -371,6 +375,10 @@ if (strcmp(gotoblas_corename(), "armv9sme") == 0
 #endif
 
 
+#ifdef SME_LEVEL3 /* arm64 SME: recursive blocking on the SME GEMM kernel */
+  if (s3_syrk_hook(0, uplo, trans, &args)) return;
+#endif
+
   IDEBUG_START;
 
   FUNCTION_PROFILE_START();

interface/trsm.c

@@ -91,6 +91,10 @@
 extern char* gotoblas_corename(void);                                
 #endif                                                               
 
+#if defined(ARCH_ARM64)
+#include "../kernel/arm64/sme_level3.h"
+#endif
+
 static int (*trsm[])(blas_arg_t *, BLASLONG *, BLASLONG *, FLOAT *, FLOAT *, BLASLONG) = {
 #ifndef TRMM
   TRSM_LNUU, TRSM_LNUN, TRSM_LNLU, TRSM_LNLN,
@@ -392,6 +396,10 @@ if (strcmp(gotoblas_corename(), "armv9sme") == 0
 #endif
 
   if ((args.m == 0) || (args.n == 0)) return;
+#ifdef SME_LEVEL3 /* arm64 SME: recursive blocking on the SME GEMM kernel */
+  if (s3_trxm_hook(side, uplo, trans, unit, &args)) return;
+#endif
+
   IDEBUG_START;
 
   FUNCTION_PROFILE_START();

kernel/arm64/sgemm_direct_performant.c

@@ -1,4 +1,10 @@
 #include "common.h"
+#if defined(__clang__) && defined(__ARM_FEATURE_SME)
+#if (defined(__apple_build_version__) && __clang_major__ >= 17) || (!defined(__apple_build_version__) && __clang_major__ >= 18)
+#define HAVE_SME2_GEMM 1
+#include "sme2_gemm_detect.h"
+#endif
+#endif
 /* helper for the direct sgemm code adapted from Arjan van der Ven's x86_64 version */
 
 int CNAME(BLASLONG M, BLASLONG N, BLASLONG K)
@@ -7,6 +13,11 @@ int CNAME(BLASLONG M, BLASLONG N, BLASLONG K)
 		return 0;
 
 	unsigned long long mnk = (unsigned long long)M * (unsigned long long)N * (unsigned long long)K;
+#ifdef HAVE_SME2_GEMM
+	/* the SME2 kernel behind SME_SGEMM_KERNEL is faster from about 64^3 (Apple M4) */
+	if (mnk >= 64ULL * 64ULL * 64ULL && s2_usable())
+		return 0;
+#endif
 	/* benchmark performance on M4 peaks around 512 and crosses the graph of the NEON SGEMM at about 3100  */
 	if (mnk >= 3100ULL * 3100ULL * 3100ULL)
 		return 0;

kernel/arm64/sme2_gemm_detect.h

@@ -0,0 +1,59 @@
+/***************************************************************************
+Copyright (c) 2026, The OpenBLAS Project
+All rights reserved.
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in
+the documentation and/or other materials provided with the
+distribution.
+3. Neither the name of the OpenBLAS project nor the names of
+its contributors may be used to endorse or promote products
+derived from this software without specific prior written permission.
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE OPENBLAS PROJECT OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
+GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
+THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*****************************************************************************/
+
+/* Runtime check for the SME2 GEMM of sme2_gemm_impl.h; also used to keep the SME1 direct sgemm out of its way. */
+#ifndef SME2_GEMM_DETECT_H
+#define SME2_GEMM_DETECT_H
+
+#include <arm_sme.h>
+#if defined(__APPLE__)
+#include <sys/sysctl.h>
+#elif defined(__linux__)
+#include <sys/auxv.h>
+#endif
+
+/* SME2 with a 512-bit streaming vector length (the kernels assume 16 fp32 lanes). */
+static inline int s2_usable(void) {
+  static int ok = -1;
+  if (ok < 0) {
+    int sme2 = 0;
+#if defined(__APPLE__)
+    int v = 0;
+    size_t len = sizeof(v);
+    sme2 = sysctlbyname("hw.optional.arm.FEAT_SME2", &v, &len, NULL, 0) == 0 && v;
+#elif defined(__linux__)
+    sme2 = (getauxval(AT_HWCAP2) & (1UL << 37)) != 0; /* HWCAP2_SME2 */
+#endif
+    ok = sme2 && svcntsw() == 16;
+  }
+  return ok;
+}
+
+/* The kernel keeps M, N and K in int; larger problems (INTERFACE64) stay on the other kernels. */
+#define S2_FITS(m, n, k) ((m) < (1L << 30) && (n) < (1L << 30) && (k) < (1L << 30))
+
+#endif