diff --git a/interface/gemm.c b/interface/gemm.c index b2e1815911..d6d43c635c 100644 --- a/interface/gemm.c +++ b/interface/gemm.c @@ -582,13 +582,21 @@ else if (order == CblasRowMajor && k==lda && n==ldb && n==ldc && TransA == Cblas #if !defined(BFLOAT16) && !defined(HFLOAT16) #if defined(ARCH_ARM64) && (defined(USE_SGEMM_KERNEL_DIRECT)||defined(DYNAMIC_ARCH)) +int sme_kernel = 1; #if defined(DYNAMIC_ARCH) -if (strcmp(gotoblas_corename(), "armv9sme") == 0 +sme_kernel = strcmp(gotoblas_corename(), "armv9sme") == 0 #if defined(__clang__) || strcmp(gotoblas_corename(), "vortexm4") == 0 #endif -) +; #endif //defined dynarch +#if USE_SMALL_MATRIX_OPT && !defined(COMPLEX) +/* a target's small-matrix kernels come before its SME kernels */ +if (sme_kernel && *(FLOAT *)args.alpha != ZERO && + GEMM_SMALL_MATRIX_PERMIT(transa, transb, args.m, args.n, args.k, *(FLOAT *)args.alpha, *(FLOAT *)args.beta)) + sme_kernel = 0; +#endif +if (sme_kernel) { char* TA,*TB; if (transa & 1) diff --git a/kernel/arm64/KERNEL.VORTEXM4 b/kernel/arm64/KERNEL.VORTEXM4 index 13a0a38c59..736b647ef5 100644 --- a/kernel/arm64/KERNEL.VORTEXM4 +++ b/kernel/arm64/KERNEL.VORTEXM4 @@ -51,4 +51,23 @@ SSYMMUCOPY_M = symm_copy_sme2.c SSYMMLCOPY_M = symm_copy_sme2.c DSYMMUCOPY_M = symm_copy_sme2.c DSYMMLCOPY_M = symm_copy_sme2.c + +SGEMM_SMALL_M_PERMIT = gemm_small_kernel_permit_vortexm4.c +SGEMM_SMALL_K_NN = gemm_small_kernel_nn_neon.c +SGEMM_SMALL_K_NT = gemm_small_kernel_nt_neon.c +SGEMM_SMALL_K_TN = gemm_small_kernel_tn_neon.c +SGEMM_SMALL_K_TT = gemm_small_kernel_tt_neon.c +SGEMM_SMALL_K_B0_NN = gemm_small_kernel_nn_neon.c +SGEMM_SMALL_K_B0_NT = gemm_small_kernel_nt_neon.c +SGEMM_SMALL_K_B0_TN = gemm_small_kernel_tn_neon.c +SGEMM_SMALL_K_B0_TT = gemm_small_kernel_tt_neon.c +DGEMM_SMALL_M_PERMIT = gemm_small_kernel_permit_vortexm4.c +DGEMM_SMALL_K_NN = gemm_small_kernel_nn_neon.c +DGEMM_SMALL_K_NT = gemm_small_kernel_nt_neon.c +DGEMM_SMALL_K_TN = gemm_small_kernel_tn_neon.c +DGEMM_SMALL_K_TT = gemm_small_kernel_tt_neon.c +DGEMM_SMALL_K_B0_NN = gemm_small_kernel_nn_neon.c +DGEMM_SMALL_K_B0_NT = gemm_small_kernel_nt_neon.c +DGEMM_SMALL_K_B0_TN = gemm_small_kernel_tn_neon.c +DGEMM_SMALL_K_B0_TT = gemm_small_kernel_tt_neon.c endif diff --git a/kernel/arm64/gemm_small_kernel_neon.h b/kernel/arm64/gemm_small_kernel_neon.h new file mode 100644 index 0000000000..18e825817a --- /dev/null +++ b/kernel/arm64/gemm_small_kernel_neon.h @@ -0,0 +1,238 @@ +/*************************************************************************** +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. +*****************************************************************************/ + + +/* Register-blocked NEON GEMM for small problems (SMALL_MATRIX_OPT): blocks of C in registers, no packing. */ + +#ifndef GEMM_SMALL_KERNEL_NEON_H +#define GEMM_SMALL_KERNEL_NEON_H + +#include +#ifdef DOUBLE +#define SK_NV 2 +typedef float64x2_t sk_v; +#define SK_LD vld1q_f64 +#define SK_ST vst1q_f64 +#define SK_DUP vdupq_n_f64 +#define SK_FMA_L(acc, a, b, l) vfmaq_laneq_f64(acc, a, b, l) +#define SK_MUL vmulq_f64 +#define SK_FMA vfmaq_f64 +#else +#define SK_NV 4 +typedef float32x4_t sk_v; +#define SK_LD vld1q_f32 +#define SK_ST vst1q_f32 +#define SK_DUP vdupq_n_f32 +#define SK_FMA_L(acc, a, b, l) vfmaq_laneq_f32(acc, a, b, l) +#define SK_MUL vmulq_f32 +#define SK_FMA vfmaq_f32 +#endif +#define SK_MR (2 * SK_NV) +#define SK_NR 4 + +#define SK_INL static inline __attribute__((always_inline)) + +SK_INL sk_v sk_fmal(sk_v acc, sk_v a, sk_v b, const int t) { + switch (t) { + case 0: return SK_FMA_L(acc, a, b, 0); +#if SK_NV > 2 + case 2: return SK_FMA_L(acc, a, b, 2); + case 3: return SK_FMA_L(acc, a, b, 3); +#endif + default: return SK_FMA_L(acc, a, b, 1); + } +} + +/* C block of MV vectors x NR columns at (a, b, c); the first skip rows of the block are left unchanged in C. */ +/* As in the reference GEMM for A op(B), alpha scales op(B) before the products. */ +SK_INL void sk_blk(const int TB, const int MV, const int NR, int skip, BLASLONG k, const FLOAT *a, BLASLONG lda, + const FLOAT *b, BLASLONG ldb, FLOAT alpha, FLOAT beta, FLOAT *c, BLASLONG ldc) { + const sk_v va = SK_DUP(alpha), vb = SK_DUP(beta); + sk_v acc[4][2]; + for (int j = 0; j < NR; ++j) + for (int h = 0; h < MV; ++h) acc[j][h] = SK_DUP(0); + BLASLONG l = 0; + if (!TB) /* B(l, j) at b[l + j * ldb]: SK_NV depth steps per load of each column */ + for (; l + SK_NV <= k; l += SK_NV) { + sk_v bv[4]; + for (int j = 0; j < NR; ++j) bv[j] = SK_MUL(SK_LD(b + l + j * ldb), va); +#pragma clang loop unroll(full) + for (int t = 0; t < SK_NV; ++t) { + sk_v av[2]; + for (int h = 0; h < MV; ++h) av[h] = SK_LD(a + (l + t) * lda + h * SK_NV); +#pragma clang loop unroll(full) + for (int j = 0; j < NR; ++j) + for (int h = 0; h < MV; ++h) acc[j][h] = sk_fmal(acc[j][h], av[h], bv[j], t); + } + } + else if (NR == 4) /* B(l, j) at b[j + l * ldb]: the four columns of one depth step in 4 / SK_NV loads */ + for (; l < k; ++l) { + sk_v av[2], bv[4 / SK_NV]; + for (int h = 0; h < MV; ++h) av[h] = SK_LD(a + l * lda + h * SK_NV); + for (int q = 0; q < 4 / SK_NV; ++q) bv[q] = SK_MUL(SK_LD(b + l * ldb + q * SK_NV), va); +#pragma clang loop unroll(full) + for (int j = 0; j < 4; ++j) + for (int h = 0; h < MV; ++h) acc[j][h] = sk_fmal(acc[j][h], av[h], bv[j / SK_NV], j % SK_NV); + } + const BLASLONG sl = TB ? ldb : 1, sj = TB ? 1 : ldb; + for (; l < k; ++l) { + sk_v av[2]; + for (int h = 0; h < MV; ++h) av[h] = SK_LD(a + l * lda + h * SK_NV); + for (int j = 0; j < NR; ++j) { + const sk_v bj = SK_DUP(alpha * b[l * sl + j * sj]); + for (int h = 0; h < MV; ++h) acc[j][h] = SK_FMA(acc[j][h], av[h], bj); + } + } + for (int j = 0; j < NR; ++j) { + FLOAT *cj = c + j * ldc; + for (int h = 0; h < MV; ++h) { + FLOAT *ch = cj + h * SK_NV; + const sk_v r = beta == (FLOAT)0 ? acc[j][h] : SK_FMA(acc[j][h], SK_LD(ch), vb); + if (skip == 0) SK_ST(ch, r); + else { + FLOAT tmp[SK_NV]; + SK_ST(tmp, r); + for (int i = skip; i < SK_NV; ++i) ch[i] = tmp[i]; + } + } + } +} + +/* Rows fewer than one vector: plain loops. */ +static inline void sk_rows(int tb, BLASLONG m, BLASLONG n, BLASLONG k, const FLOAT *a, BLASLONG lda, const FLOAT *b, + BLASLONG ldb, FLOAT alpha, FLOAT beta, FLOAT *c, BLASLONG ldc) { + const BLASLONG sl = tb ? ldb : 1, sj = tb ? 1 : ldb; + for (BLASLONG j = 0; j < n; ++j) + for (BLASLONG i = 0; i < m; ++i) { + FLOAT s = 0; + for (BLASLONG l = 0; l < k; ++l) s += a[i + l * lda] * (alpha * b[l * sl + j * sj]); + c[i + j * ldc] = s + (beta == (FLOAT)0 ? (FLOAT)0 : beta * c[i + j * ldc]); + } +} + +/* One block column of NR columns: two-vector row blocks, one-vector blocks, and the row tail as an overlapped block. */ +SK_INL void sk_col(const int TB, const int NR, BLASLONG m, BLASLONG k, const FLOAT *a, BLASLONG lda, const FLOAT *b, + BLASLONG ldb, FLOAT alpha, FLOAT beta, FLOAT *c, BLASLONG ldc) { + BLASLONG i = 0; + for (; i + 2 * SK_NV <= m; i += 2 * SK_NV) sk_blk(TB, 2, NR, 0, k, a + i, lda, b, ldb, alpha, beta, c + i, ldc); + if (i + SK_NV <= m) { + sk_blk(TB, 1, NR, 0, k, a + i, lda, b, ldb, alpha, beta, c + i, ldc); + i += SK_NV; + } + if (i < m) { + if (m >= SK_NV) + sk_blk(TB, 1, NR, (int)(SK_NV - (m - i)), k, a + m - SK_NV, lda, b, ldb, alpha, beta, c + m - SK_NV, ldc); + else + sk_rows(TB, m, NR, k, a, lda, b, ldb, alpha, beta, c, ldc); + } +} + +/* C = alpha A op(B) + beta C with A column-major; tb: B transposed. */ +static inline void sk_gemm_n(int tb, BLASLONG m, BLASLONG n, BLASLONG k, FLOAT alpha, const FLOAT *a, BLASLONG lda, + const FLOAT *b, BLASLONG ldb, FLOAT beta, FLOAT *c, BLASLONG ldc) { + const BLASLONG sj = tb ? 1 : ldb; + BLASLONG j = 0; + if (tb) + for (; j + 4 <= n; j += 4) sk_col(1, 4, m, k, a, lda, b + j * sj, ldb, alpha, beta, c + j * ldc, ldc); + else + for (; j + 4 <= n; j += 4) sk_col(0, 4, m, k, a, lda, b + j * sj, ldb, alpha, beta, c + j * ldc, ldc); +#define SK_TAIL(NR) \ + if (tb) sk_col(1, NR, m, k, a, lda, b + j * sj, ldb, alpha, beta, c + j * ldc, ldc); \ + else sk_col(0, NR, m, k, a, lda, b + j * sj, ldb, alpha, beta, c + j * ldc, ldc); + switch (n - j) { + case 3: SK_TAIL(3) break; + case 2: SK_TAIL(2) break; + case 1: SK_TAIL(1) break; + } +#undef SK_TAIL +} + +#ifdef DOUBLE +#define SK_ADDV vaddvq_f64 +#else +#define SK_ADDV vaddvq_f32 +#endif + +/* C block of R rows x NR columns for C = alpha A^T B + beta C: dot products along the contiguous depth. */ +SK_INL void sk_dot(const int R, const int NR, BLASLONG k, const FLOAT *a, BLASLONG lda, const FLOAT *b, BLASLONG ldb, + FLOAT alpha, FLOAT beta, FLOAT *c, BLASLONG ldc) { + sk_v acc[4][4]; + for (int r = 0; r < R; ++r) + for (int j = 0; j < NR; ++j) acc[r][j] = SK_DUP(0); + BLASLONG l = 0; + for (; l + SK_NV <= k; l += SK_NV) { + sk_v av[4], bv[4]; + for (int r = 0; r < R; ++r) av[r] = SK_LD(a + r * lda + l); + for (int j = 0; j < NR; ++j) bv[j] = SK_LD(b + j * ldb + l); + for (int r = 0; r < R; ++r) + for (int j = 0; j < NR; ++j) acc[r][j] = SK_FMA(acc[r][j], av[r], bv[j]); + } + for (int r = 0; r < R; ++r) + for (int j = 0; j < NR; ++j) { + FLOAT s = SK_ADDV(acc[r][j]); + for (BLASLONG q = l; q < k; ++q) s += a[r * lda + q] * b[j * ldb + q]; + FLOAT *cc = c + r + j * ldc; + *cc = alpha * s + (beta == (FLOAT)0 ? (FLOAT)0 : beta * *cc); + } +} + +SK_INL void sk_dot_rows(const int NR, BLASLONG m, BLASLONG k, const FLOAT *a, BLASLONG lda, const FLOAT *b, + BLASLONG ldb, FLOAT alpha, FLOAT beta, FLOAT *c, BLASLONG ldc) { + BLASLONG i = 0; + for (; i + 4 <= m; i += 4) sk_dot(4, NR, k, a + i * lda, lda, b, ldb, alpha, beta, c + i, ldc); + switch (m - i) { + case 3: sk_dot(3, NR, k, a + i * lda, lda, b, ldb, alpha, beta, c + i, ldc); break; + case 2: sk_dot(2, NR, k, a + i * lda, lda, b, ldb, alpha, beta, c + i, ldc); break; + case 1: sk_dot(1, NR, k, a + i * lda, lda, b, ldb, alpha, beta, c + i, ldc); break; + } +} + +/* C = alpha A^T B + beta C. */ +static inline void sk_gemm_tn(BLASLONG m, BLASLONG n, BLASLONG k, FLOAT alpha, const FLOAT *a, BLASLONG lda, + const FLOAT *b, BLASLONG ldb, FLOAT beta, FLOAT *c, BLASLONG ldc) { + BLASLONG j = 0; + for (; j + 4 <= n; j += 4) sk_dot_rows(4, m, k, a, lda, b + j * ldb, ldb, alpha, beta, c + j * ldc, ldc); + switch (n - j) { + case 3: sk_dot_rows(3, m, k, a, lda, b + j * ldb, ldb, alpha, beta, c + j * ldc, ldc); break; + case 2: sk_dot_rows(2, m, k, a, lda, b + j * ldb, ldb, alpha, beta, c + j * ldc, ldc); break; + case 1: sk_dot_rows(1, m, k, a, lda, b + j * ldb, ldb, alpha, beta, c + j * ldc, ldc); break; + } +} + +/* C = alpha A^T B^T + beta C: plain loops. */ +static inline void sk_gemm_tt(BLASLONG m, BLASLONG n, BLASLONG k, FLOAT alpha, const FLOAT *a, BLASLONG lda, + const FLOAT *b, BLASLONG ldb, FLOAT beta, FLOAT *c, BLASLONG ldc) { + for (BLASLONG j = 0; j < n; ++j) + for (BLASLONG i = 0; i < m; ++i) { + FLOAT s = 0; + for (BLASLONG l = 0; l < k; ++l) s += a[l + i * lda] * b[j + l * ldb]; + FLOAT *cc = c + i + j * ldc; + *cc = alpha * s + (beta == (FLOAT)0 ? (FLOAT)0 : beta * *cc); + } +} + +#endif diff --git a/kernel/arm64/gemm_small_kernel_nn_neon.c b/kernel/arm64/gemm_small_kernel_nn_neon.c new file mode 100644 index 0000000000..151c6e9125 --- /dev/null +++ b/kernel/arm64/gemm_small_kernel_nn_neon.c @@ -0,0 +1,46 @@ +/*************************************************************************** +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. +*****************************************************************************/ + +/* Small-matrix GEMM C = alpha A B + beta C (gemm_small_kernel_neon.h). */ + +#include "common.h" +#include "gemm_small_kernel_neon.h" + +#ifdef B0 +int CNAME(BLASLONG M, BLASLONG N, BLASLONG K, IFLOAT *A, BLASLONG lda, FLOAT alpha, IFLOAT *B, BLASLONG ldb, FLOAT *C, + BLASLONG ldc) +#else +int CNAME(BLASLONG M, BLASLONG N, BLASLONG K, IFLOAT *A, BLASLONG lda, FLOAT alpha, IFLOAT *B, BLASLONG ldb, FLOAT beta, + FLOAT *C, BLASLONG ldc) +#endif +{ +#ifdef B0 + const FLOAT beta = ZERO; +#endif + sk_gemm_n(0, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc); + return 0; +} diff --git a/kernel/arm64/gemm_small_kernel_nt_neon.c b/kernel/arm64/gemm_small_kernel_nt_neon.c new file mode 100644 index 0000000000..6fd437f0ff --- /dev/null +++ b/kernel/arm64/gemm_small_kernel_nt_neon.c @@ -0,0 +1,46 @@ +/*************************************************************************** +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. +*****************************************************************************/ + +/* Small-matrix GEMM C = alpha A B^T + beta C (gemm_small_kernel_neon.h). */ + +#include "common.h" +#include "gemm_small_kernel_neon.h" + +#ifdef B0 +int CNAME(BLASLONG M, BLASLONG N, BLASLONG K, IFLOAT *A, BLASLONG lda, FLOAT alpha, IFLOAT *B, BLASLONG ldb, FLOAT *C, + BLASLONG ldc) +#else +int CNAME(BLASLONG M, BLASLONG N, BLASLONG K, IFLOAT *A, BLASLONG lda, FLOAT alpha, IFLOAT *B, BLASLONG ldb, FLOAT beta, + FLOAT *C, BLASLONG ldc) +#endif +{ +#ifdef B0 + const FLOAT beta = ZERO; +#endif + sk_gemm_n(1, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc); + return 0; +} diff --git a/kernel/arm64/gemm_small_kernel_permit_vortexm4.c b/kernel/arm64/gemm_small_kernel_permit_vortexm4.c new file mode 100644 index 0000000000..85fc204c3f --- /dev/null +++ b/kernel/arm64/gemm_small_kernel_permit_vortexm4.c @@ -0,0 +1,39 @@ +/*************************************************************************** +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. +*****************************************************************************/ + +/* Small-matrix GEMM on NEON below the multiply-adds where the SME2 kernels become faster (Apple M4). */ + +#include "common.h" + +int CNAME(int transa, int transb, BLASLONG M, BLASLONG N, BLASLONG K, FLOAT alpha, FLOAT beta) { + if (M == 1 || N == 1) return 0; /* vector shapes stay on the SME2 kernels; forwarded to GEMV, LAPACK SHS test 11 fails */ +#ifdef DOUBLE + return (double)M * (double)N * (double)K < 12000.; +#else + return (double)M * (double)N * (double)K < 16000.; +#endif +} diff --git a/kernel/arm64/gemm_small_kernel_tn_neon.c b/kernel/arm64/gemm_small_kernel_tn_neon.c new file mode 100644 index 0000000000..569276c7c1 --- /dev/null +++ b/kernel/arm64/gemm_small_kernel_tn_neon.c @@ -0,0 +1,46 @@ +/*************************************************************************** +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. +*****************************************************************************/ + +/* Small-matrix GEMM C = alpha A^T B + beta C (gemm_small_kernel_neon.h). */ + +#include "common.h" +#include "gemm_small_kernel_neon.h" + +#ifdef B0 +int CNAME(BLASLONG M, BLASLONG N, BLASLONG K, IFLOAT *A, BLASLONG lda, FLOAT alpha, IFLOAT *B, BLASLONG ldb, FLOAT *C, + BLASLONG ldc) +#else +int CNAME(BLASLONG M, BLASLONG N, BLASLONG K, IFLOAT *A, BLASLONG lda, FLOAT alpha, IFLOAT *B, BLASLONG ldb, FLOAT beta, + FLOAT *C, BLASLONG ldc) +#endif +{ +#ifdef B0 + const FLOAT beta = ZERO; +#endif + sk_gemm_tn(M, N, K, alpha, A, lda, B, ldb, beta, C, ldc); + return 0; +} diff --git a/kernel/arm64/gemm_small_kernel_tt_neon.c b/kernel/arm64/gemm_small_kernel_tt_neon.c new file mode 100644 index 0000000000..0f21a9d2af --- /dev/null +++ b/kernel/arm64/gemm_small_kernel_tt_neon.c @@ -0,0 +1,46 @@ +/*************************************************************************** +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. +*****************************************************************************/ + +/* Small-matrix GEMM C = alpha A^T B^T + beta C (gemm_small_kernel_neon.h). */ + +#include "common.h" +#include "gemm_small_kernel_neon.h" + +#ifdef B0 +int CNAME(BLASLONG M, BLASLONG N, BLASLONG K, IFLOAT *A, BLASLONG lda, FLOAT alpha, IFLOAT *B, BLASLONG ldb, FLOAT *C, + BLASLONG ldc) +#else +int CNAME(BLASLONG M, BLASLONG N, BLASLONG K, IFLOAT *A, BLASLONG lda, FLOAT alpha, IFLOAT *B, BLASLONG ldb, FLOAT beta, + FLOAT *C, BLASLONG ldc) +#endif +{ +#ifdef B0 + const FLOAT beta = ZERO; +#endif + sk_gemm_tt(M, N, K, alpha, A, lda, B, ldb, beta, C, ldc); + return 0; +} diff --git a/kernel/arm64/sgemm_direct_performant.c b/kernel/arm64/sgemm_direct_performant.c index 0fc7252971..a5270f86e2 100644 --- a/kernel/arm64/sgemm_direct_performant.c +++ b/kernel/arm64/sgemm_direct_performant.c @@ -9,8 +9,8 @@ int CNAME(BLASLONG M, BLASLONG N, BLASLONG K) 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()) + /* the SME2 kernel behind SME_SGEMM_KERNEL is faster from about 64^3, and its NEON loops below about 22^3 (Apple M4) */ + if ((mnk >= 64ULL * 64ULL * 64ULL || mnk < 10000ULL) && s2_usable()) return 0; #endif /* benchmark performance on M4 peaks around 512 and crosses the graph of the NEON SGEMM at about 3100 */ diff --git a/kernel/arm64/sme_dgemm_kernel.c b/kernel/arm64/sme_dgemm_kernel.c index 6edd82fde0..9517a421b2 100644 --- a/kernel/arm64/sme_dgemm_kernel.c +++ b/kernel/arm64/sme_dgemm_kernel.c @@ -1190,8 +1190,7 @@ void CNAME(const char *transa, const char *transb, const BLASLONG m, const BLASL bool trans_a = (*transa == 'T' || *transa == 't' || *transa == 'C' || *transa == 'c'); bool trans_b = (*transb == 'T' || *transb == 't' || *transb == 'C' || *transb == 'c'); #if defined(HAVE_SME2_GEMM) - /* below 5000 multiply-adds the existing kernel is faster when it needs no edge tiles (M and N multiples of 16) */ - if (s2_usable() && S2_FITS(m, n, k) && !((double)m * (double)n * (double)k < 5000. && m % 16 == 0 && n % 16 == 0)) { + if (s2_usable() && S2_FITS(m, n, k)) { s2_gemm(trans_a, trans_b, m, n, k, *alpha, a, lda, b, ldb, *beta, c, ldc); return; }