From f518e23902c7015628f1da5b51dc371417da37cf Mon Sep 17 00:00:00 2001 From: michael tesch Date: Wed, 21 Oct 2026 23:46:20 -0500 Subject: [PATCH 1/2] arm64: SME2 SGEMV and DGEMV kernels for Apple M4 (VORTEXM4) GEMV N keeps up to 64 vectors of y in ZA and adds four columns at a time with the indexed multi-vector FMLA, column after column as the reference BLAS. GEMV T accumulates 16 columns per block in ZA (one ZA vector per column, strided multi-vector loads) and sums them through a tile transpose. Both enter streaming mode once per call. Small problems, T column tails and CPUs without SME2 use NEON code (sgemv_n_neon.c for fp32 N). --- kernel/arm64/KERNEL.VORTEXM4 | 5 + kernel/arm64/gemv_n_sme2.c | 172 ++++++++++++++++++++++++++++++ kernel/arm64/gemv_t_sme2.c | 201 +++++++++++++++++++++++++++++++++++ kernel/arm64/sme2_gemv.h | 88 +++++++++++++++ 4 files changed, 466 insertions(+) create mode 100644 kernel/arm64/gemv_n_sme2.c create mode 100644 kernel/arm64/gemv_t_sme2.c create mode 100644 kernel/arm64/sme2_gemv.h diff --git a/kernel/arm64/KERNEL.VORTEXM4 b/kernel/arm64/KERNEL.VORTEXM4 index 13a0a38c59..8f50822af8 100644 --- a/kernel/arm64/KERNEL.VORTEXM4 +++ b/kernel/arm64/KERNEL.VORTEXM4 @@ -51,4 +51,9 @@ SSYMMUCOPY_M = symm_copy_sme2.c SSYMMLCOPY_M = symm_copy_sme2.c DSYMMUCOPY_M = symm_copy_sme2.c DSYMMLCOPY_M = symm_copy_sme2.c + +SGEMVNKERNEL = gemv_n_sme2.c +SGEMVTKERNEL = gemv_t_sme2.c +DGEMVNKERNEL = gemv_n_sme2.c +DGEMVTKERNEL = gemv_t_sme2.c endif diff --git a/kernel/arm64/gemv_n_sme2.c b/kernel/arm64/gemv_n_sme2.c new file mode 100644 index 0000000000..05d01e0df6 --- /dev/null +++ b/kernel/arm64/gemv_n_sme2.c @@ -0,0 +1,172 @@ +/*************************************************************************** +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. +*****************************************************************************/ + + +/* GEMV N, y += alpha A x: SME2 with up to 64 vectors of y updated in ZA; NEON for small problems and without SME2. */ + +#include "common.h" +#include "sme2_gemv.h" + +/* Smallest m and m * n for the SME2 path, including the streaming mode switch. */ +#ifndef GVN_SME2_MIN +#define GVN_SME2_MIN 1024 +#endif +#ifndef GVN_SME2_MINM +#define GVN_SME2_MINM 16 +#endif + +#ifndef DOUBLE +/* fp32 small problems and CPUs without SME2: the NEON kernel of other arm64 targets. */ +#pragma push_macro("CNAME") +#undef CNAME +#define CNAME sgemv_n_neon_vortexm4 +#include "sgemv_n_neon.c" +#pragma pop_macro("CNAME") +#else +/* y (contiguous) += alpha A x, four columns at a time. */ +static void gvn_neon(BLASLONG m, BLASLONG n, FLOAT alpha, const FLOAT *a, BLASLONG lda, const FLOAT *x, BLASLONG incx, + FLOAT *y) { + BLASLONG j = 0; + for (; j + 4 <= n; j += 4) { + const FLOAT *a0 = a + j * lda, *a1 = a0 + lda, *a2 = a1 + lda, *a3 = a2 + lda; + const FLOAT x0 = alpha * x[j * incx], x1 = alpha * x[(j + 1) * incx]; + const FLOAT x2 = alpha * x[(j + 2) * incx], x3 = alpha * x[(j + 3) * incx]; + BLASLONG i = 0; + for (; i + 2 * GV_NV <= m; i += 2 * GV_NV) { + gv_v y0 = GV_LD(y + i), y1 = GV_LD(y + i + GV_NV); + y0 = GV_FMAN(y0, GV_LD(a0 + i), x0); + y1 = GV_FMAN(y1, GV_LD(a0 + i + GV_NV), x0); + y0 = GV_FMAN(y0, GV_LD(a1 + i), x1); + y1 = GV_FMAN(y1, GV_LD(a1 + i + GV_NV), x1); + y0 = GV_FMAN(y0, GV_LD(a2 + i), x2); + y1 = GV_FMAN(y1, GV_LD(a2 + i + GV_NV), x2); + y0 = GV_FMAN(y0, GV_LD(a3 + i), x3); + y1 = GV_FMAN(y1, GV_LD(a3 + i + GV_NV), x3); + GV_ST(y + i, y0); + GV_ST(y + i + GV_NV, y1); + } + for (; i < m; i++) y[i] += a0[i] * x0 + a1[i] * x1 + a2[i] * x2 + a3[i] * x3; + } + for (; j < n; j++) { + const FLOAT *aj = a + j * lda; + const FLOAT xj = alpha * x[j * incx]; + for (BLASLONG i = 0; i < m; i++) y[i] += aj[i] * xj; + } +} +#endif + +#ifdef HAVE_SME2_GEMM +#pragma clang attribute push(__attribute__((target("sme2"))), apply_to = function) + +/* Columns per 128-bit segment of x, used by the indexed FMLA. */ +#ifdef DOUBLE +#define GVN_NQ 2 +#define GVN_LANES(g, pc, aj, lda, xq, r) \ + GV_MLAL4(g, S2_LD4(pc, aj), xq, 0); \ + if (r > 1) GV_MLAL4(g, S2_LD4(pc, (aj) + (lda)), xq, 1); +#else +#define GVN_NQ 4 +#define GVN_LANES(g, pc, aj, lda, xq, r) \ + GV_MLAL4(g, S2_LD4(pc, aj), xq, 0); \ + if (r > 1) GV_MLAL4(g, S2_LD4(pc, (aj) + (lda)), xq, 1); \ + if (r > 2) GV_MLAL4(g, S2_LD4(pc, (aj) + 2 * (lda)), xq, 2); \ + if (r > 3) GV_MLAL4(g, S2_LD4(pc, (aj) + 3 * (lda)), xq, 3); +#endif + +/* ZA groups 0..ng-1 += A (alpha x) for a block of mb rows, column after column as the reference BLAS; full: mb = GV_G * GV_GL. */ +S2_INL void gvn_block(int ng, int full, BLASLONG mb, BLASLONG n, FLOAT alpha, const FLOAT *a, BLASLONG lda, + const FLOAT *x) S2_S S2_ZA { + BLASLONG j = 0; + for (; j + GVN_NQ <= n; j += GVN_NQ) { + const s2_v xq = S2_MUL(GV_LDRQ(S2_PT(), x + j), alpha); + const FLOAT *aj = a + j * lda; + for (int g = 0; g < ng; g++) { + const svcount_t pc = full ? S2_CT() : S2_CW(g * GV_GL, mb); + GVN_LANES(g, pc, aj + g * GV_GL, lda, xq, GVN_NQ) + } + } + if (j < n) { + const int r = (int)(n - j); + const s2_v xq = S2_MUL(GV_LDRQ(S2_PW(j, n), x + j), alpha); + const FLOAT *aj = a + j * lda; + for (int g = 0; g < ng; g++) { + const svcount_t pc = full ? S2_CT() : S2_CW(g * GV_GL, mb); + GVN_LANES(g, pc, aj + g * GV_GL, lda, xq, r) + } + } +} + +__arm_new("za") __arm_locally_streaming static void gvn_sme2(BLASLONG m, BLASLONG n, FLOAT alpha, const FLOAT *a, + BLASLONG lda, const FLOAT *x, FLOAT *y) { + for (BLASLONG i0 = 0; i0 < m; i0 += GV_G * GV_GL) { + const BLASLONG mb = m - i0 < GV_G * GV_GL ? m - i0 : GV_G * GV_GL; + const int ng = (int)((mb + GV_GL - 1) / GV_GL); + for (int g = 0; g < ng; g++) GV_WRZ4(g, S2_LD4(S2_CW(g * GV_GL, mb), y + i0 + g * GV_GL)); + if (ng == GV_G && mb == GV_G * GV_GL) + gvn_block(GV_G, 1, mb, n, alpha, a + i0, lda, x); + else + gvn_block(ng, 0, mb, n, alpha, a + i0, lda, x); + for (int g = 0; g < ng; g++) S2_ST4(S2_CW(g * GV_GL, mb), y + i0 + g * GV_GL, GV_RDZ4(g)); + } +} + +#pragma clang attribute pop +#endif + +int CNAME(BLASLONG m, BLASLONG n, BLASLONG dummy1, FLOAT alpha, FLOAT *a, BLASLONG lda, FLOAT *x, BLASLONG inc_x, + FLOAT *y, BLASLONG inc_y, FLOAT *buffer) { + if (m < 1 || n < 1 || alpha == ZERO) return 0; + FLOAT *yy = y; +#ifdef HAVE_SME2_GEMM + if (m >= GVN_SME2_MINM && m * n >= GVN_SME2_MIN && s2_usable()) { + FLOAT *xx = x; + if (inc_x != 1) { + xx = buffer; + gv_gather(n, x, inc_x, xx); + buffer += n; + } + if (inc_y != 1) { + yy = buffer; + gv_gather(m, y, inc_y, yy); + } + gvn_sme2(m, n, alpha, a, lda, xx, yy); + if (inc_y != 1) gv_scatter(m, yy, y, inc_y); + return 0; + } +#endif +#ifndef DOUBLE + return sgemv_n_neon_vortexm4(m, n, dummy1, alpha, a, lda, x, inc_x, y, inc_y, buffer); +#else + if (inc_y != 1) { + yy = buffer; + gv_gather(m, y, inc_y, yy); + } + gvn_neon(m, n, alpha, a, lda, x, inc_x, yy); + if (inc_y != 1) gv_scatter(m, yy, y, inc_y); + return 0; +#endif +} diff --git a/kernel/arm64/gemv_t_sme2.c b/kernel/arm64/gemv_t_sme2.c new file mode 100644 index 0000000000..13a7f3c375 --- /dev/null +++ b/kernel/arm64/gemv_t_sme2.c @@ -0,0 +1,201 @@ +/*************************************************************************** +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. +*****************************************************************************/ + + +/* GEMV T, y += alpha A^T x: SME2 with 16 columns at a time accumulated in ZA; NEON for small problems and without SME2. */ + +#include "common.h" +#include "sme2_gemv.h" + +/* Smallest m and m * n for the SME2 path, which takes the columns in blocks of 16. */ +#ifndef GVT_SME2_MIN +#define GVT_SME2_MIN 4096 +#endif +#ifndef GVT_SME2_MINM +#define GVT_SME2_MINM 48 +#endif + +/* y += alpha A^T x with a contiguous x, four columns at a time. */ +static void gvt_neon(BLASLONG m, BLASLONG n, FLOAT alpha, const FLOAT *a, BLASLONG lda, const FLOAT *x, FLOAT *y, + BLASLONG incy) { + BLASLONG j = 0; + for (; j + 4 <= n; j += 4) { + const FLOAT *a0 = a + j * lda, *a1 = a0 + lda, *a2 = a1 + lda, *a3 = a2 + lda; + gv_v s0 = GV_ZERO(), s1 = s0, s2 = s0, s3 = s0, u0 = s0, u1 = s0, u2 = s0, u3 = s0; + BLASLONG i = 0; + for (; i + 2 * GV_NV <= m; i += 2 * GV_NV) { + const gv_v x0 = GV_LD(x + i), x1 = GV_LD(x + i + GV_NV); + s0 = GV_FMA(s0, GV_LD(a0 + i), x0); + u0 = GV_FMA(u0, GV_LD(a0 + i + GV_NV), x1); + s1 = GV_FMA(s1, GV_LD(a1 + i), x0); + u1 = GV_FMA(u1, GV_LD(a1 + i + GV_NV), x1); + s2 = GV_FMA(s2, GV_LD(a2 + i), x0); + u2 = GV_FMA(u2, GV_LD(a2 + i + GV_NV), x1); + s3 = GV_FMA(s3, GV_LD(a3 + i), x0); + u3 = GV_FMA(u3, GV_LD(a3 + i + GV_NV), x1); + } + FLOAT r0 = GV_ADDV(s0, u0), r1 = GV_ADDV(s1, u1), r2 = GV_ADDV(s2, u2), r3 = GV_ADDV(s3, u3); + for (; i < m; i++) { + r0 += a0[i] * x[i]; + r1 += a1[i] * x[i]; + r2 += a2[i] * x[i]; + r3 += a3[i] * x[i]; + } + y[j * incy] += alpha * r0; + y[(j + 1) * incy] += alpha * r1; + y[(j + 2) * incy] += alpha * r2; + y[(j + 3) * incy] += alpha * r3; + } + for (; j < n; j++) { + const FLOAT *aj = a + j * lda; + gv_v s0 = GV_ZERO(), u0 = s0; + BLASLONG i = 0; + for (; i + 2 * GV_NV <= m; i += 2 * GV_NV) { + s0 = GV_FMA(s0, GV_LD(aj + i), GV_LD(x + i)); + u0 = GV_FMA(u0, GV_LD(aj + i + GV_NV), GV_LD(x + i + GV_NV)); + } + FLOAT r0 = GV_ADDV(s0, u0); + for (; i < m; i++) r0 += aj[i] * x[i]; + y[j * incy] += alpha * r0; + } +} + +#ifdef HAVE_SME2_GEMM +#pragma clang attribute push(__attribute__((target("sme2"))), apply_to = function) + +#ifdef DOUBLE +#define GVT_T "d" +#define GVT_LD "ld1d" +#define GVT_SH "3" +#else +#define GVT_T "s" +#define GVT_LD "ld1w" +#define GVT_SH "2" +#endif + +/* Full block of 16 columns, ZA row s + 16 k accumulates column 4 s + k over rows i..i+GV_GL-1; strided loads feed the FMLA without moves. */ +S2_INL void gvt_chunk(BLASLONG i, BLASLONG m, const FLOAT *x, const FLOAT *a, BLASLONG ldb) S2_S S2_ZA { + __asm__ volatile( + "whilelt pn8." GVT_T ", %[i], %[m], vlx4\n" + GVT_LD " {z0." GVT_T "-z3." GVT_T "}, pn8/z, [%[x], %[i], lsl #" GVT_SH "]\n" + "mov w8, #0\n" + "mov x10, %[a]\n" + "1:\n" + "add x11, x10, %[ldb]\n" + "add x12, x11, %[ldb]\n" + "add x13, x12, %[ldb]\n" + GVT_LD " {z16." GVT_T ", z20." GVT_T ", z24." GVT_T ", z28." GVT_T "}, pn8/z, [x10, %[i], lsl #" GVT_SH "]\n" + GVT_LD " {z17." GVT_T ", z21." GVT_T ", z25." GVT_T ", z29." GVT_T "}, pn8/z, [x11, %[i], lsl #" GVT_SH "]\n" + GVT_LD " {z18." GVT_T ", z22." GVT_T ", z26." GVT_T ", z30." GVT_T "}, pn8/z, [x12, %[i], lsl #" GVT_SH "]\n" + GVT_LD " {z19." GVT_T ", z23." GVT_T ", z27." GVT_T ", z31." GVT_T "}, pn8/z, [x13, %[i], lsl #" GVT_SH "]\n" + "add x10, x13, %[ldb]\n" + "fmla za." GVT_T "[w8, 0, vgx4], {z16." GVT_T "-z19." GVT_T "}, z0." GVT_T "\n" + "fmla za." GVT_T "[w8, 0, vgx4], {z20." GVT_T "-z23." GVT_T "}, z1." GVT_T "\n" + "fmla za." GVT_T "[w8, 0, vgx4], {z24." GVT_T "-z27." GVT_T "}, z2." GVT_T "\n" + "fmla za." GVT_T "[w8, 0, vgx4], {z28." GVT_T "-z31." GVT_T "}, z3." GVT_T "\n" + "add w8, w8, #1\n" + "cmp w8, #4\n" + "b.lt 1b\n" + : + : [i] "r"(i), [m] "r"(m), [x] "r"(x), [a] "r"(a), [ldb] "r"(ldb) + : "x8", "x10", "x11", "x12", "x13", "p8", "z0", "z1", "z2", "z3", "z16", "z17", "z18", "z19", "z20", "z21", "z22", + "z23", "z24", "z25", "z26", "z27", "z28", "z29", "z30", "z31", "za", "cc", "memory"); +} + +#define GVT_SUM4(r) GV_ADD(GV_ADD(S2_GET4(r, 0), S2_GET4(r, 1)), GV_ADD(S2_GET4(r, 2), S2_GET4(r, 3))) +#define GVT_VSUM(t, s, acc) acc = GV_ADD(acc, GVT_SUM4(S2_RDV4(t, s))) + +/* y[0..16) += alpha * sums of the 16 column vectors in ZA groups 0..3 (column 4 s + k in vector k of group s). */ +S2_INL void gvt_reduce(FLOAT alpha, FLOAT *y) S2_S S2_ZA { + const s2_v4 r0 = GV_RDZ4(0), r1 = GV_RDZ4(1), r2 = GV_RDZ4(2), r3 = GV_RDZ4(3); + s2_v t0 = S2_ZEROV(); +#ifdef DOUBLE + s2_v t1 = S2_ZEROV(); + S2_WRH4(0, 0, r0); + S2_WRH4(0, 4, r1); + S2_WRH4(1, 0, r2); + S2_WRH4(1, 4, r3); + GVT_VSUM(0, 0, t0); + GVT_VSUM(0, 4, t0); + GVT_VSUM(1, 0, t1); + GVT_VSUM(1, 4, t1); + S2_ST1(S2_PT(), y + S2_VL, GV_MLAV(S2_LD1(S2_PT(), y + S2_VL), t1, alpha)); +#else + s2_v t1 = S2_ZEROV(); + S2_WRH4(0, 0, r0); + S2_WRH4(0, 4, r1); + S2_WRH4(0, 8, r2); + S2_WRH4(0, 12, r3); + GVT_VSUM(0, 0, t0); + GVT_VSUM(0, 4, t1); + GVT_VSUM(0, 8, t0); + GVT_VSUM(0, 12, t1); + t0 = GV_ADD(t0, t1); +#endif + S2_ST1(S2_PT(), y, GV_MLAV(S2_LD1(S2_PT(), y), t0, alpha)); +} + +/* n is a multiple of 16. */ +__arm_new("za") __arm_locally_streaming static void gvt_sme2(BLASLONG m, BLASLONG n, FLOAT alpha, const FLOAT *a, + BLASLONG lda, const FLOAT *x, FLOAT *y) { + for (BLASLONG j0 = 0; j0 < n; j0 += 16) { + svzero_za(); + for (BLASLONG i = 0; i < m; i += GV_GL) gvt_chunk(i, m, x, a + j0 * lda, lda * (BLASLONG)sizeof(FLOAT)); + gvt_reduce(alpha, y + j0); + } +} + +#pragma clang attribute pop +#endif + +int CNAME(BLASLONG m, BLASLONG n, BLASLONG dummy1, FLOAT alpha, FLOAT *a, BLASLONG lda, FLOAT *x, BLASLONG inc_x, + FLOAT *y, BLASLONG inc_y, FLOAT *buffer) { + if (m < 1 || n < 1 || alpha == ZERO) return 0; + FLOAT *xx = x; + if (inc_x != 1) { + xx = buffer; + gv_gather(m, x, inc_x, xx); + buffer += m; + } +#ifdef HAVE_SME2_GEMM + const BLASLONG n16 = n & ~(BLASLONG)15; + if (n16 && m >= GVT_SME2_MINM && m * n >= GVT_SME2_MIN && s2_usable()) { + FLOAT *yy = y; + if (inc_y != 1) { + yy = buffer; + gv_gather(n16, y, inc_y, yy); + } + gvt_sme2(m, n16, alpha, a, lda, xx, yy); + if (inc_y != 1) gv_scatter(n16, yy, y, inc_y); + a += n16 * lda; + y += n16 * inc_y; + n -= n16; + } +#endif + gvt_neon(m, n, alpha, a, lda, xx, y, inc_y); + return 0; +} diff --git a/kernel/arm64/sme2_gemv.h b/kernel/arm64/sme2_gemv.h new file mode 100644 index 0000000000..c84e81d46e --- /dev/null +++ b/kernel/arm64/sme2_gemv.h @@ -0,0 +1,88 @@ +/*************************************************************************** +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. +*****************************************************************************/ + + +/* Shared parts of the SME2 GEMV kernels: NEON helpers for small problems and CPUs without SME2, ZA macros. */ +#ifndef SME2_GEMV_H +#define SME2_GEMV_H + +#include +#include "sme2_gemm_detect.h" + +#ifdef DOUBLE +#define GV_NV 2 +typedef float64x2_t gv_v; +#define GV_LD(p) vld1q_f64(p) +#define GV_ST(p, v) vst1q_f64(p, v) +#define GV_FMA(c, a, b) vfmaq_f64(c, a, b) +#define GV_FMAN(c, a, s) vfmaq_n_f64(c, a, s) +#define GV_ZERO() vdupq_n_f64(0) +#define GV_ADDV(a, b) vaddvq_f64(vaddq_f64(a, b)) +#else +#define GV_NV 4 +typedef float32x4_t gv_v; +#define GV_LD(p) vld1q_f32(p) +#define GV_ST(p, v) vst1q_f32(p, v) +#define GV_FMA(c, a, b) vfmaq_f32(c, a, b) +#define GV_FMAN(c, a, s) vfmaq_n_f32(c, a, s) +#define GV_ZERO() vdupq_n_f32(0) +#define GV_ADDV(a, b) vaddvq_f32(vaddq_f32(a, b)) +#endif + +/* Copy n elements with stride inc (the pointer as the level-2 interface passes it) to a contiguous d, and back. */ +static inline void gv_gather(BLASLONG n, const FLOAT *s, BLASLONG inc, FLOAT *d) { + for (BLASLONG i = 0; i < n; i++) d[i] = s[i * inc]; +} +static inline void gv_scatter(BLASLONG n, const FLOAT *s, FLOAT *d, BLASLONG inc) { + for (BLASLONG i = 0; i < n; i++) d[i * inc] = s[i]; +} + +#ifdef HAVE_SME2_GEMM +#include "sme2_gemm_tile.h" + +#ifdef DOUBLE +#define GV_RDZ4(s) svread_za64_f64_vg1x4(s) +#define GV_WRZ4(s, v) svwrite_za64_f64_vg1x4(s, v) +#define GV_MLAL4(s, a, z, l) svmla_lane_za64_f64_vg1x4(s, a, z, l) +#define GV_LDRQ(p, q) svld1rq_f64(p, q) +#define GV_ADD(a, b) svadd_f64_x(svptrue_b64(), a, b) +#define GV_MLAV(c, v, s) svmla_n_f64_x(svptrue_b64(), c, v, s) +#else +#define GV_RDZ4(s) svread_za32_f32_vg1x4(s) +#define GV_WRZ4(s, v) svwrite_za32_f32_vg1x4(s, v) +#define GV_MLAL4(s, a, z, l) svmla_lane_za32_f32_vg1x4(s, a, z, l) +#define GV_LDRQ(p, q) svld1rq_f32(p, q) +#define GV_ADD(a, b) svadd_f32_x(svptrue_b32(), a, b) +#define GV_MLAV(c, v, s) svmla_n_f32_x(svptrue_b32(), c, v, s) +#endif + +/* ZA vg1x4 groups: 16 slice indices of 4 vectors each, 64 vectors in all. */ +#define GV_G 16 +#define GV_GL (4 * S2_VL) +#endif + +#endif From e385aa3b762ea8108fb41bc068681bf98a107e39 Mon Sep 17 00:00:00 2001 From: michael tesch Date: Wed, 21 Oct 2026 23:46:20 -0500 Subject: [PATCH 2/2] interface: one GEMV thread per SME unit on VORTEXM4 The SME2 GEMV kernels share the SME unit of a P-cluster, so more threads than units are slower. Use num_cpu_avail(3), which sme_level3_threads() caps at the number of units, above 8 MB of matrix, and one thread below. --- interface/gemv.c | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/interface/gemv.c b/interface/gemv.c index d718e17011..ac74f2a34e 100644 --- a/interface/gemv.c +++ b/interface/gemv.c @@ -104,6 +104,13 @@ static inline int get_gemv_optimal_nthreads_neoversev2(BLASLONG MN, int ncpu) { } #endif +#if defined(__clang__) && (defined(DYNAMIC_ARCH) || defined(VORTEXM4)) +/* Apple M4: ncpu from num_cpu_avail(3) is one thread per SME unit; one unit is faster below about 8 MB of matrix. */ +static inline int get_gemv_optimal_nthreads_vortexm4(BLASLONG MN, int ncpu) { + return MN * (BLASLONG)sizeof(FLOAT) <= (8L << 20) ? 1 : ncpu; +} +#endif + static inline int get_gemv_optimal_nthreads(BLASLONG MN) { int ncpu = num_cpu_avail(3); #if defined(_WIN64) && defined(_M_ARM64) @@ -115,6 +122,8 @@ static inline int get_gemv_optimal_nthreads(BLASLONG MN) { return get_gemv_optimal_nthreads_neoversev1(MN, ncpu); #elif defined(NEOVERSEV2) && !defined(COMPLEX) && !defined(DOUBLE) && !defined(BFLOAT16) return get_gemv_optimal_nthreads_neoversev2(MN, ncpu); +#elif defined(VORTEXM4) && defined(__clang__) && !defined(COMPLEX) && !defined(BFLOAT16) + return get_gemv_optimal_nthreads_vortexm4(MN, ncpu); #elif defined(DYNAMIC_ARCH) && !defined(COMPLEX) && !defined(BFLOAT16) if (strcmp(gotoblas_corename(), "neoversev1") == 0) { return get_gemv_optimal_nthreads_neoversev1(MN, ncpu); @@ -122,6 +131,11 @@ static inline int get_gemv_optimal_nthreads(BLASLONG MN) { if (strcmp(gotoblas_corename(), "neoversev2") == 0) { return get_gemv_optimal_nthreads_neoversev2(MN, ncpu); } +#ifdef __clang__ + if (strcmp(gotoblas_corename(), "vortexm4") == 0) { + return get_gemv_optimal_nthreads_vortexm4(MN, ncpu); + } +#endif #endif if ( MN < 115200L * GEMM_MULTITHREAD_THRESHOLD )