From 382d1863216b314fbb6916c9bd78b7e801e714ba Mon Sep 17 00:00:00 2001 From: michael tesch Date: Wed, 21 Oct 2026 23:43:36 -0500 Subject: [PATCH 1/2] arm64: SME2 complex level-3 kernels for Apple M4 (VORTEXM4) zgemm_kernel_sme2.c is the CGEMM/ZGEMM and TRMM kernel of the level-3 driver. The packed panels keep the generic interleaved layout; FMOPA on the interleaved vectors accumulates re*re, im*re, re*im and im*im in separate ZA rows, and FCADD/FCMLA form the complex result and apply alpha on the store. All four conjugation variants come from the same code. CGEMM unroll is 16x16, ZGEMM 16x8 (ZA full, B one streaming vector); the copies are the generic ones, with new generic zgemm_tcopy_16.c and zlaswp_ncopy_16.c. ztrsm_kernel_sme2.c does the TRSM update in the GEMM kernel and the solve in plain or NEON code. SYMM, HEMM, SYRK, HERK, SYR2K, HER2K, TRMM, TRSM and the complex LAPACK routines on the level-3 driver use these kernels. Complex GEMM itself still goes to the SME1 direct kernels from interface/gemm.c. --- cmake/prebuild.cmake | 12 +- kernel/arm64/KERNEL.VORTEXM4 | 25 ++- kernel/arm64/zgemm_kernel_sme2.c | 297 +++++++++++++++++++++++++++++++ kernel/arm64/ztrsm_kernel_sme2.c | 251 ++++++++++++++++++++++++++ kernel/generic/zgemm_tcopy_16.c | 44 +++++ kernel/generic/zlaswp_ncopy_16.c | 59 ++++++ param.h | 13 +- 7 files changed, 690 insertions(+), 11 deletions(-) create mode 100644 kernel/arm64/zgemm_kernel_sme2.c create mode 100644 kernel/arm64/ztrsm_kernel_sme2.c create mode 100644 kernel/generic/zgemm_tcopy_16.c create mode 100644 kernel/generic/zlaswp_ncopy_16.c diff --git a/cmake/prebuild.cmake b/cmake/prebuild.cmake index 61ed48c000..ba2c42c103 100644 --- a/cmake/prebuild.cmake +++ b/cmake/prebuild.cmake @@ -1278,16 +1278,20 @@ endif () set(SGEMM_UNROLL_N 16) set(DGEMM_UNROLL_M 64) set(DGEMM_UNROLL_N 8) + set(CGEMM_UNROLL_M 16) + set(CGEMM_UNROLL_N 16) + set(ZGEMM_UNROLL_M 16) + set(ZGEMM_UNROLL_N 8) else () set(SGEMM_UNROLL_M 16) set(SGEMM_UNROLL_N 4) set(DGEMM_UNROLL_M 8) set(DGEMM_UNROLL_N 4) + set(CGEMM_UNROLL_M 8) + set(CGEMM_UNROLL_N 4) + set(ZGEMM_UNROLL_M 4) + set(ZGEMM_UNROLL_N 4) endif () - set(CGEMM_UNROLL_M 8) - set(CGEMM_UNROLL_N 4) - set(ZGEMM_UNROLL_M 4) - set(ZGEMM_UNROLL_N 4) set(SYMV_P 16) elseif ("${TCORE}" STREQUAL "A64FX") file(APPEND ${TARGET_CONF_TEMP} diff --git a/kernel/arm64/KERNEL.VORTEXM4 b/kernel/arm64/KERNEL.VORTEXM4 index 13a0a38c59..084534eed9 100644 --- a/kernel/arm64/KERNEL.VORTEXM4 +++ b/kernel/arm64/KERNEL.VORTEXM4 @@ -1,6 +1,6 @@ include $(KERNELDIR)/KERNEL.NEOVERSEN1 -# SME2 level-3 kernels (GEMM_UNROLL_M 64, GEMM_UNROLL_N one streaming vector); GCC builds keep the NEON kernels +# SME2 level-3 kernels (GEMM_UNROLL_M 64 real, 16 complex; GEMM_UNROLL_N one streaming vector); GCC keeps NEON ifneq ($(C_COMPILER), GCC) SGEMMKERNEL = gemm_kernel_sme2.c STRMMKERNEL = gemm_kernel_sme2.c @@ -51,4 +51,27 @@ SSYMMUCOPY_M = symm_copy_sme2.c SSYMMLCOPY_M = symm_copy_sme2.c DSYMMUCOPY_M = symm_copy_sme2.c DSYMMLCOPY_M = symm_copy_sme2.c + +CGEMMKERNEL = zgemm_kernel_sme2.c +CTRMMKERNEL = zgemm_kernel_sme2.c +CGEMMINCOPY = ../generic/zgemm_ncopy_$(CGEMM_UNROLL_M).c +CGEMMITCOPY = ../generic/zgemm_tcopy_$(CGEMM_UNROLL_M).c +CGEMMONCOPY = ../generic/zgemm_ncopy_$(CGEMM_UNROLL_N).c +CGEMMOTCOPY = ../generic/zgemm_tcopy_$(CGEMM_UNROLL_N).c + +ZGEMMKERNEL = zgemm_kernel_sme2.c +ZTRMMKERNEL = zgemm_kernel_sme2.c +ZGEMMINCOPY = ../generic/zgemm_ncopy_$(ZGEMM_UNROLL_M).c +ZGEMMITCOPY = ../generic/zgemm_tcopy_$(ZGEMM_UNROLL_M).c +ZGEMMONCOPY = ../generic/zgemm_ncopy_$(ZGEMM_UNROLL_N).c +ZGEMMOTCOPY = ../generic/zgemm_tcopy_$(ZGEMM_UNROLL_N).c + +CTRSMKERNEL_LN = ztrsm_kernel_sme2.c +CTRSMKERNEL_LT = ztrsm_kernel_sme2.c +CTRSMKERNEL_RN = ztrsm_kernel_sme2.c +CTRSMKERNEL_RT = ztrsm_kernel_sme2.c +ZTRSMKERNEL_LN = ztrsm_kernel_sme2.c +ZTRSMKERNEL_LT = ztrsm_kernel_sme2.c +ZTRSMKERNEL_RN = ztrsm_kernel_sme2.c +ZTRSMKERNEL_RT = ztrsm_kernel_sme2.c endif diff --git a/kernel/arm64/zgemm_kernel_sme2.c b/kernel/arm64/zgemm_kernel_sme2.c new file mode 100644 index 0000000000..12639d23d7 --- /dev/null +++ b/kernel/arm64/zgemm_kernel_sme2.c @@ -0,0 +1,297 @@ +/*************************************************************************** +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. +*****************************************************************************/ + +/* SME2 complex GEMM kernel of the level-3 driver, the TRMM kernel with TRMMKERNEL; plain loops without SME2. */ + +#include "common.h" +#include "sme2_gemm_detect.h" + +/* Conjugation of A (bit 1) and of B (bit 0), as the generic complex kernels. */ +#if defined(NN) || defined(NT) || defined(TN) || defined(TT) +#define S2Z_CJ 0 +#elif defined(NR) || defined(NC) || defined(TR) || defined(TC) +#define S2Z_CJ 1 +#elif defined(RN) || defined(RT) || defined(CN) || defined(CT) +#define S2Z_CJ 2 +#else +#define S2Z_CJ 3 +#endif + +#ifdef HAVE_SME2_GEMM +#define S2Z_SC __arm_streaming_compatible +#else +#define S2Z_SC +#endif + +/* Width of the panel at j0 of n as the generic copies lay it out: full panels of u, then 8, 4, 2, 1. */ +static inline int s2z_w(BLASLONG n, BLASLONG j0, int u) S2Z_SC { + BLASLONG r = n - j0; + if (r >= u) return u; + int w = 1; + while (2 * w <= r) w *= 2; + return w; +} + +/* Depth range of the block at (i0, j0): the TRMM kernel skips the zero part of the triangle as trmmkernel_4x4.c. */ +static inline void s2z_krange(BLASLONG i0, BLASLONG j0, int wm, int wn, BLASLONG k, BLASLONG offset, BLASLONG *ks, + BLASLONG *kl) S2Z_SC { +#ifdef TRMMKERNEL +#ifdef LEFT + const BLASLONG off = offset + i0; +#else + const BLASLONG off = j0 - offset; +#endif +#if (defined(LEFT) && !defined(TRANSA)) || (!defined(LEFT) && defined(TRANSA)) + BLASLONG s = off, e = k; +#elif defined(LEFT) + BLASLONG s = 0, e = off + wm; +#else + BLASLONG s = 0, e = off + wn; +#endif + if (s < 0) s = 0; + if (s > k) s = k; + if (e > k) e = k; + if (e < s) e = s; + *ks = s; + *kl = e - s; +#else + *ks = 0; + *kl = k; +#endif +} + +#ifdef HAVE_SME2_GEMM +#include "sme2_gemm_tile.h" + +/* A panels of 2 (fp32) or 4 (fp64) vectors, B panels of 2 vectors; complex numbers stay interleaved. */ +#define S2Z_TA (2 * GEMM_DEFAULT_UNROLL_M / S2_VL) +#define S2Z_TB (2 * GEMM_DEFAULT_UNROLL_N / S2_VL) +#if S2Z_TB != 2 || S2Z_TA * S2Z_TB != S2_NT +#error "the SME2 complex kernel needs GEMM_UNROLL_N of one streaming vector and UNROLL_M x UNROLL_N to fill ZA" +#endif + +#pragma clang attribute push(__attribute__((target("sme2"))), apply_to = function) + +#ifdef DOUBLE +#define S2Z_CADD(x, y, r) svcadd_f64_x(svptrue_b64(), x, y, r) +#define S2Z_CMLA(c, x, y, r) svcmla_f64_x(svptrue_b64(), c, x, y, r) +#define S2Z_DUPQ(r, i) svdupq_n_f64(r, i) +#else +#define S2Z_CADD(x, y, r) svcadd_f32_x(svptrue_b32(), x, y, r) +#define S2Z_CMLA(c, x, y, r) svcmla_f32_x(svptrue_b32(), c, x, y, r) +#define S2Z_DUPQ(r, i) svdupq_n_f32(r, i, r, i) +#endif + +/* Tile ib * TA + ia += B vector ib (rows) x A vector ia; ZA rows 2q, 2q + 1: A re b_q, A im b_q (interleaved). */ +S2_INL void s2z_mac(const int TA, const int FULL, BLASLONG kb, const FLOAT *a, int ma, const FLOAT *b, + int nb) S2_S S2_ZA { + const svbool_t pt = S2_PT(); + BLASLONG l = 0; + if (FULL) { + for (; l + 2 <= kb; l += 2) { + const s2_v4 bb = s2_ld4(b + l * 2 * S2_VL); + s2_v4 a0, a1; + if (TA == 4) { + a0 = s2_ld4(a + l * 4 * S2_VL); + a1 = s2_ld4(a + (l + 1) * 4 * S2_VL); + } else { + a0 = s2_ld4(a + l * 2 * S2_VL); + a1 = a0; + } + S2_UNROLL for (int U = 0; U < 2; ++U) { + S2_UNROLL for (int ib = 0; ib < 2; ++ib) { + const s2_v bv = s2_get4(bb, 2 * U + ib); + S2_UNROLL for (int ia = 0; ia < TA; ++ia) { + const s2_v av = TA == 4 ? s2_get4(U ? a1 : a0, ia) : s2_get4(a0, 2 * U + ia); + s2_mopa(ib * TA + ia, pt, bv, av); + } + } + } + } + } + const svbool_t pb0 = S2_PW(0, nb), pb1 = S2_PW(S2_VL, nb); + const int two = nb > S2_VL; + for (; l < kb; ++l) { + const FLOAT *ar = a + l * ma, *br = b + l * nb; + const s2_v b0 = S2_LD1(pb0, br), b1 = two ? S2_LD1(pb1, br + S2_VL) : b0; + S2_UNROLL for (int ia = 0; ia < TA; ++ia) { + const s2_v av = S2_LD1(S2_PW(ia * S2_VL, ma), ar + ia * S2_VL); + s2_mopa(ia, pt, b0, av); + if (two) s2_mopa(TA + ia, pt, b1, av); + } + } +} + +/* C (+)= alpha T for the two ZA rows x, y of one column of op(A) op(B) (T = x + i y with the conjugations). */ +S2_INL void s2z_put(FLOAT *c, svbool_t pg, s2_v x, s2_v y, s2_v al, const int add) S2_S { +#if S2Z_CJ == 0 || S2Z_CJ == 3 + const s2_v t = S2Z_CADD(x, y, 90); +#else + const s2_v t = S2Z_CADD(x, y, 270); +#endif + s2_v v = add ? S2_LD1(pg, c) : S2_ZEROV(); + v = S2Z_CMLA(v, t, al, 0); +#if S2Z_CJ < 2 + v = S2Z_CMLA(v, t, al, 90); +#else + v = S2Z_CMLA(v, t, al, 270); +#endif + S2_ST1(pg, c, v); +} + +/* C (+)= alpha op(A) op(B) for wm x wn complex numbers from the ZA layout of s2z_mac; FULL: whole block. */ +S2_INL void s2z_store(const int TA, const int FULL, FLOAT *C, BLASLONG ldc, int wm, int wn, const s2_v al, + const int add) S2_S S2_ZA { + const svbool_t pt = S2_PT(); + if (FULL) { /* four ZA rows (two columns) per move, all tiles at once for independent chains */ + S2_UNROLL for (int ib = 0; ib < 2; ++ib) + S2_UNROLL for (int r = 0; r < S2_VL; r += 4) { + FLOAT *c = C + (BLASLONG)(ib * S2_VL + r) * ldc; + S2_UNROLL for (int ia = 0; ia < TA; ++ia) { + const s2_v4 z = s2_rdh4(ib * TA + ia, r); + s2z_put(c + ia * S2_VL, pt, s2_get4(z, 0), s2_get4(z, 1), al, add); + s2z_put(c + 2 * ldc + ia * S2_VL, pt, s2_get4(z, 2), s2_get4(z, 3), al, add); + } + } + return; + } + const int ma = 2 * wm; + for (int q = 0; q < wn; ++q) { + const int ib = 2 * q / S2_VL, r = 2 * q % S2_VL; + S2_UNROLL for (int ia = 0; ia < TA; ++ia) + s2z_put(C + 2 * q * ldc + ia * S2_VL, S2_PW(ia * S2_VL, ma), s2_rdh(ib * TA + ia, r), s2_rdh(ib * TA + ia, r + 1), + al, add); + } +} + +typedef void (*s2z_bfn)(BLASLONG, const FLOAT *, int, const FLOAT *, int, FLOAT *, BLASLONG, FLOAT, FLOAT) S2_S S2_ZA; + +#ifdef TRMMKERNEL +#define S2Z_ADD 0 +#else +#define S2Z_ADD 1 +#endif + +/* Block of wm x wn complex numbers of C, TA A vectors, FULL: both panels at full width. */ +#define S2Z_DEF(TA, FULL) \ + S2_NOINL void s2z_blk_##TA##_##FULL(BLASLONG kb, const FLOAT *a, int wm, const FLOAT *b, int wn, FLOAT *C, \ + BLASLONG ldc, FLOAT ar, FLOAT ai) S2_S S2_ZA { \ + svzero_za(); \ + s2z_mac(TA, FULL, kb, a, 2 * wm, b, 2 * wn); \ + s2z_store(TA, FULL, C, ldc, wm, wn, S2Z_DUPQ(ar, ai), S2Z_ADD); \ + } +S2Z_DEF(1, 0) +S2Z_DEF(2, 0) +#ifdef DOUBLE +S2Z_DEF(4, 0) +S2Z_DEF(4, 1) +static const s2z_bfn s2z_blk[4] = {s2z_blk_1_0, s2z_blk_2_0, s2z_blk_4_0, s2z_blk_4_0}; +#define S2Z_FULL s2z_blk_4_1 +#else +S2Z_DEF(2, 1) +static const s2z_bfn s2z_blk[2] = {s2z_blk_1_0, s2z_blk_2_0}; +#define S2Z_FULL s2z_blk_2_1 +#endif + +__arm_new("za") __arm_locally_streaming static void s2z_drive(BLASLONG m, BLASLONG n, BLASLONG k, FLOAT ar, + FLOAT ai, const FLOAT *sa, const FLOAT *sb, FLOAT *C, + BLASLONG ldc, BLASLONG offset) { + for (BLASLONG j0 = 0; j0 < n;) { + const int wn = s2z_w(n, j0, GEMM_DEFAULT_UNROLL_N); + for (BLASLONG i0 = 0; i0 < m;) { + const int wm = s2z_w(m, i0, GEMM_DEFAULT_UNROLL_M); + BLASLONG ks, kl; + s2z_krange(i0, j0, wm, wn, k, offset, &ks, &kl); + const FLOAT *a = sa + 2 * (i0 * k + ks * wm), *b = sb + 2 * (j0 * k + ks * wn); + FLOAT *c = C + 2 * (i0 + j0 * ldc); + if (wm == GEMM_DEFAULT_UNROLL_M && wn == GEMM_DEFAULT_UNROLL_N) + S2Z_FULL(kl, a, wm, b, wn, c, ldc, ar, ai); + else + s2z_blk[(2 * wm + S2_VL - 1) / S2_VL - 1](kl, a, wm, b, wn, c, ldc, ar, ai); + i0 += wm; + } + j0 += wn; + } +} + +#pragma clang attribute pop +#endif + +static void s2z_loops(BLASLONG m, BLASLONG n, BLASLONG k, FLOAT alr, FLOAT ali, const FLOAT *sa, const FLOAT *sb, + FLOAT *C, BLASLONG ldc, BLASLONG offset) { + const FLOAT sa_i = (S2Z_CJ & 2) ? -1 : 1, sb_i = (S2Z_CJ & 1) ? -1 : 1; + for (BLASLONG j0 = 0; j0 < n;) { + const int wn = s2z_w(n, j0, GEMM_DEFAULT_UNROLL_N); + for (BLASLONG i0 = 0; i0 < m;) { + const int wm = s2z_w(m, i0, GEMM_DEFAULT_UNROLL_M); + BLASLONG ks, kl; + s2z_krange(i0, j0, wm, wn, k, offset, &ks, &kl); + const FLOAT *a = sa + 2 * (i0 * k + ks * wm), *b = sb + 2 * (j0 * k + ks * wn); + for (int j = 0; j < wn; ++j) + for (int i = 0; i < wm; ++i) { + FLOAT sr = 0, si = 0; + for (BLASLONG l = 0; l < kl; ++l) { + const FLOAT xr = a[2 * (l * wm + i)], xi = sa_i * a[2 * (l * wm + i) + 1]; + const FLOAT yr = b[2 * (l * wn + j)], yi = sb_i * b[2 * (l * wn + j) + 1]; + sr += xr * yr - xi * yi; + si += xr * yi + xi * yr; + } + FLOAT *c = C + 2 * (i0 + i + (j0 + j) * ldc); +#ifdef TRMMKERNEL + c[0] = alr * sr - ali * si; + c[1] = alr * si + ali * sr; +#else + c[0] += alr * sr - ali * si; + c[1] += alr * si + ali * sr; +#endif + } + i0 += wm; + } + j0 += wn; + } +} + +int CNAME(BLASLONG bm, BLASLONG bn, BLASLONG bk, FLOAT alphar, FLOAT alphai, FLOAT *ba, FLOAT *bb, FLOAT *C, + BLASLONG ldc +#ifdef TRMMKERNEL + , BLASLONG offset +#endif +) { +#ifndef TRMMKERNEL + const BLASLONG offset = 0; + if (bk <= 0) return 0; +#endif + if (bm <= 0 || bn <= 0) return 0; +#ifdef HAVE_SME2_GEMM + if (s2_usable()) { + s2z_drive(bm, bn, bk, alphar, alphai, ba, bb, C, ldc, offset); + return 0; + } +#endif + s2z_loops(bm, bn, bk, alphar, alphai, ba, bb, C, ldc, offset); + return 0; +} diff --git a/kernel/arm64/ztrsm_kernel_sme2.c b/kernel/arm64/ztrsm_kernel_sme2.c new file mode 100644 index 0000000000..25192d45c2 --- /dev/null +++ b/kernel/arm64/ztrsm_kernel_sme2.c @@ -0,0 +1,251 @@ +/*************************************************************************** +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. +*****************************************************************************/ + +/* Complex TRSM kernels for the SME2 complex GEMM kernel: the update in GEMM_KERNEL, the solve in plain or NEON code. */ + +#include "common.h" +#include + +#if defined(LN) || defined(LT) +#define ZT_LEFT 1 +#endif +#if defined(LN) || defined(RT) +#define ZT_BACK 1 +#endif +#ifdef CONJ +#define ZT_CJ 1 +#else +#define ZT_CJ 0 +#endif +#ifdef ZT_LEFT +#define ZT_GEMM (ZT_CJ ? GEMM_KERNEL_L : GEMM_KERNEL_N) +#else +#define ZT_GEMM (ZT_CJ ? GEMM_KERNEL_R : GEMM_KERNEL_N) +#endif + +#ifdef DOUBLE +typedef float64x2_t zt_v; +#define ZT_N 1 +#define ZT_LD(p) vld1q_f64(p) +#define ZT_ST(p, v) vst1q_f64(p, v) +#define ZT_DUP(r, i) ((float64x2_t){r, i}) +#define ZT_ZERO() vdupq_n_f64(0) +#define ZT_C0(c, n, m) vcmlaq_f64(c, n, m) +#define ZT_C90(c, n, m) vcmlaq_rot90_f64(c, n, m) +#else +typedef float32x4_t zt_v; +#define ZT_N 2 +#define ZT_LD(p) vld1q_f32(p) +#define ZT_ST(p, v) vst1q_f32(p, v) +#define ZT_DUP(r, i) ((float32x4_t){r, i, r, i}) +#define ZT_ZERO() vdupq_n_f32(0) +#define ZT_C0(c, n, m) vcmlaq_f32(c, n, m) +#define ZT_C90(c, n, m) vcmlaq_rot90_f32(c, n, m) +#endif + +/* Panel width at j0 of n as the generic copies lay it out: full panels of u, then 8, 4, 2, 1. */ +static inline int zt_w(BLASLONG n, BLASLONG j0, int u) { + BLASLONG r = n - j0; + if (r >= u) return u; + int w = 1; + while (2 * w <= r) w *= 2; + return w; +} + +#ifndef ZT_LEFT +/* c[0..len) += v[0..len) s for complex numbers. */ +static inline void zt_axpy(int len, FLOAT *c, const FLOAT *v, FLOAT sr, FLOAT si) { + const zt_v m = ZT_DUP(sr, si); + int r = 0; + for (; r + ZT_N <= len; r += ZT_N) { + const zt_v n = ZT_LD(v + 2 * r); + ZT_ST(c + 2 * r, ZT_C90(ZT_C0(ZT_LD(c + 2 * r), n, m), n, m)); + } + for (; r < len; ++r) { + const FLOAT vr = v[2 * r], vi = v[2 * r + 1]; + c[2 * r] += vr * sr - vi * si; + c[2 * r + 1] += vr * si + vi * sr; + } +} + +/* x = c = c s for len complex numbers. */ +static inline void zt_scale(int len, FLOAT *c, FLOAT *x, FLOAT sr, FLOAT si) { + const zt_v m = ZT_DUP(sr, si); + int r = 0; + for (; r + ZT_N <= len; r += ZT_N) { + const zt_v n = ZT_LD(c + 2 * r); + const zt_v v = ZT_C90(ZT_C0(ZT_ZERO(), n, m), n, m); + ZT_ST(c + 2 * r, v); + ZT_ST(x + 2 * r, v); + } + for (; r < len; ++r) { + const FLOAT cr = c[2 * r], ci = c[2 * r + 1]; + x[2 * r] = c[2 * r] = cr * sr - ci * si; + x[2 * r + 1] = c[2 * r + 1] = cr * si + ci * sr; + } +} +#endif + +#ifdef ZT_LEFT +/* Left solve on real and imaginary planes of the block with rows contiguous: plain vectorizable loops along wn. */ +static inline void zt_solve_l(int wm, const int wn, const FLOAT *t, FLOAT *x, FLOAT *c, BLASLONG ldc) { + const FLOAT sg = ZT_CJ ? -1 : 1; + FLOAT pr[GEMM_DEFAULT_UNROLL_M * GEMM_DEFAULT_UNROLL_N], pi[GEMM_DEFAULT_UNROLL_M * GEMM_DEFAULT_UNROLL_N]; + for (int j = 0; j < wn; ++j) + for (int i = 0; i < wm; ++i) { + pr[i * wn + j] = c[2 * (i + j * ldc)]; + pi[i * wn + j] = c[2 * (i + j * ldc) + 1]; + } + for (int u = 0; u < wm; ++u) { +#ifdef ZT_BACK + const int i = wm - 1 - u, r0 = 0, r1 = i; +#else + const int i = u, r0 = i + 1, r1 = wm; +#endif + const FLOAT *ti = t + 2 * (BLASLONG)i * wm; + const FLOAT dr = ti[2 * i], di = sg * ti[2 * i + 1]; + FLOAT *restrict xr = pr + i * wn, *restrict xi = pi + i * wn, *restrict xo = x + 2 * (BLASLONG)i * wn; + for (int j = 0; j < wn; ++j) { + const FLOAT a = xr[j], b = xi[j]; + xr[j] = a * dr - b * di; + xi[j] = a * di + b * dr; + xo[2 * j] = xr[j]; + xo[2 * j + 1] = xi[j]; + } + for (int r = r0; r < r1; ++r) { + const FLOAT tr = ti[2 * r], tm = sg * ti[2 * r + 1]; + FLOAT *restrict yr = pr + r * wn, *restrict yi = pi + r * wn; + for (int j = 0; j < wn; ++j) { + yr[j] -= tr * xr[j] - tm * xi[j]; + yi[j] -= tr * xi[j] + tm * xr[j]; + } + } + } + for (int j = 0; j < wn; ++j) + for (int i = 0; i < wm; ++i) { + c[2 * (i + j * ldc)] = pr[i * wn + j]; + c[2 * (i + j * ldc) + 1] = pi[i * wn + j]; + } +} +#endif + +/* Solve of a wm x wn block against the triangle t (inverted diagonal); the solution also goes to x. */ +static void zt_solve(int wm, int wn, const FLOAT *t, FLOAT *x, FLOAT *c, BLASLONG ldc) { +#ifdef ZT_LEFT + if (wn == GEMM_DEFAULT_UNROLL_N) zt_solve_l(wm, GEMM_DEFAULT_UNROLL_N, t, x, c, ldc); + else zt_solve_l(wm, wn, t, x, c, ldc); +#else + const FLOAT sg = ZT_CJ ? -1 : 1; + for (int u = 0; u < wn; ++u) { +#ifdef ZT_BACK + const int i = wn - 1 - u, r0 = 0, r1 = i; +#else + const int i = u, r0 = i + 1, r1 = wn; +#endif + const FLOAT *ti = t + 2 * (BLASLONG)i * wn; + FLOAT *xi = x + 2 * (BLASLONG)i * wm; + zt_scale(wm, c + 2 * i * ldc, xi, ti[2 * i], sg * ti[2 * i + 1]); + for (int r = r0; r < r1; ++r) zt_axpy(wm, c + 2 * r * ldc, xi, -ti[2 * r], -sg * ti[2 * r + 1]); + } +#endif +} + +/* C block -= au bu over ku steps, then the solve. */ +static void zt_block(BLASLONG ku, FLOAT *au, FLOAT *bu, const FLOAT *t, FLOAT *x, int wm, int wn, FLOAT *c, + BLASLONG ldc) { + if (ku > 0) ZT_GEMM(wm, wn, ku, -ONE, ZERO, au, bu, c, ldc); + zt_solve(wm, wn, t, x, c, ldc); +} + +int CNAME(BLASLONG m, BLASLONG n, BLASLONG k, FLOAT dummy1, FLOAT dummy2, FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, + BLASLONG offset) { + const int um = GEMM_UNROLL_M, un = GEMM_UNROLL_N; + if (m <= 0 || n <= 0) return 0; +#ifdef ZT_LEFT + for (BLASLONG j0 = 0; j0 < n;) { + const int wn = zt_w(n, j0, un); + FLOAT *bj = b + 2 * j0 * k, *cj = c + 2 * j0 * ldc; +#ifdef LT + BLASLONG kk = offset; + for (BLASLONG i0 = 0; i0 < m;) { + const int wm = zt_w(m, i0, um); + FLOAT *aa = a + 2 * i0 * k; + zt_block(kk, aa, bj, aa + 2 * kk * wm, bj + 2 * kk * wn, wm, wn, cj + 2 * i0, ldc); + kk += wm; + i0 += wm; + } +#else + /* A panels from the last one back: the tails (1, 2, 4, ... rows), then the full panels */ + BLASLONG kk = m + offset, i1 = m; + for (int wm = 1; i1 > 0; wm = wm < um ? 2 * wm : um) { + if (wm < um && !(m & wm)) continue; + const BLASLONG i0 = i1 - wm; + FLOAT *aa = a + 2 * i0 * k; + zt_block(k - kk, aa + 2 * wm * kk, bj + 2 * wn * kk, aa + 2 * (kk - wm) * wm, bj + 2 * (kk - wm) * wn, wm, wn, + cj + 2 * i0, ldc); + kk -= wm; + i1 = i0; + } +#endif + j0 += wn; + } +#else +#ifdef RN + BLASLONG kk = -offset; + for (BLASLONG j0 = 0; j0 < n;) { + const int wn = zt_w(n, j0, un); + FLOAT *bj = b + 2 * j0 * k; + for (BLASLONG i0 = 0; i0 < m;) { + const int wm = zt_w(m, i0, um); + FLOAT *aa = a + 2 * i0 * k; + zt_block(kk, aa, bj, bj + 2 * kk * wn, aa + 2 * kk * wm, wm, wn, c + 2 * (i0 + j0 * ldc), ldc); + i0 += wm; + } + kk += wn; + j0 += wn; + } +#else + /* B panels from the last one back: the tails (1, 2, 4, ... columns), then the full panels */ + BLASLONG kk = n - offset, j1 = n; + for (int wn = 1; j1 > 0; wn = wn < un ? 2 * wn : un) { + if (wn < un && !(n & wn)) continue; + const BLASLONG j0 = j1 - wn; + FLOAT *bj = b + 2 * j0 * k; + for (BLASLONG i0 = 0; i0 < m;) { + const int wm = zt_w(m, i0, um); + FLOAT *aa = a + 2 * i0 * k; + zt_block(k - kk, aa + 2 * wm * kk, bj + 2 * wn * kk, bj + 2 * (kk - wn) * wn, aa + 2 * (kk - wn) * wm, wm, wn, + c + 2 * (i0 + j0 * ldc), ldc); + i0 += wm; + } + kk -= wn; + j1 = j0; + } +#endif +#endif + return 0; +} diff --git a/kernel/generic/zgemm_tcopy_16.c b/kernel/generic/zgemm_tcopy_16.c new file mode 100644 index 0000000000..1779540840 --- /dev/null +++ b/kernel/generic/zgemm_tcopy_16.c @@ -0,0 +1,44 @@ +/*************************************************************************** +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. +*****************************************************************************/ + +/* Panels of 16 complex columns of op(A) with contiguous rows, then 8, 4, 2, 1 (the layout of zgemm_ncopy_16.c). */ + +#include "common.h" + +int CNAME(BLASLONG m, BLASLONG n, FLOAT *a, BLASLONG lda, FLOAT *b) { + for (BLASLONG j0 = 0; j0 < n;) { + BLASLONG w = 16; + while (w > n - j0) w /= 2; + for (BLASLONG l = 0; l < m; ++l) { + const FLOAT *src = a + 2 * (j0 + l * lda); + for (BLASLONG c = 0; c < 2 * w; ++c) b[c] = src[c]; + b += 2 * w; + } + j0 += w; + } + return 0; +} diff --git a/kernel/generic/zlaswp_ncopy_16.c b/kernel/generic/zlaswp_ncopy_16.c new file mode 100644 index 0000000000..6631b45774 --- /dev/null +++ b/kernel/generic/zlaswp_ncopy_16.c @@ -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. +*****************************************************************************/ + +/* Row interchanges of LASWP on n complex columns, the rows k1..k2 copied into panels of 16 columns, then 8, 4, 2, 1. */ + +#include "common.h" + +int CNAME(BLASLONG n, BLASLONG k1, BLASLONG k2, FLOAT *a, BLASLONG lda, blasint *ipiv, FLOAT *buffer) { + if (n <= 0 || k1 > k2) return 0; + a -= 2; + ipiv += k1 - 1; + for (BLASLONG j0 = 0; j0 < n;) { + BLASLONG w = 16; + while (w > n - j0) w /= 2; + for (BLASLONG i = k1; i <= k2; ++i) { + const BLASLONG ip = ipiv[i - k1]; + for (BLASLONG c = 0; c < w; ++c) { + FLOAT *col = a + 2 * (j0 + c) * lda; + const FLOAT xr = col[2 * i], xi = col[2 * i + 1], yr = col[2 * ip], yi = col[2 * ip + 1]; + if (ip != i) { + col[2 * ip] = xr; + col[2 * ip + 1] = xi; + buffer[2 * c] = yr; + buffer[2 * c + 1] = yi; + } else { + buffer[2 * c] = xr; + buffer[2 * c + 1] = xi; + } + } + buffer += 2 * w; + } + j0 += w; + } + return 0; +} diff --git a/param.h b/param.h index 26d7bdc279..f86eedbecd 100644 --- a/param.h +++ b/param.h @@ -3573,16 +3573,17 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #define DGEMM_DEFAULT_UNROLL_M 64 #define DGEMM_DEFAULT_UNROLL_N 8 -#define CGEMM_DEFAULT_UNROLL_M 8 -#define CGEMM_DEFAULT_UNROLL_N 4 +/* kernel/arm64/zgemm_kernel_sme2.c: interleaved complex numbers, A in 2 (fp32) or 4 (fp64) vectors, B in 2 */ +#define CGEMM_DEFAULT_UNROLL_M 16 +#define CGEMM_DEFAULT_UNROLL_N 16 -#define ZGEMM_DEFAULT_UNROLL_M 4 -#define ZGEMM_DEFAULT_UNROLL_N 4 +#define ZGEMM_DEFAULT_UNROLL_M 16 +#define ZGEMM_DEFAULT_UNROLL_N 8 #define SGEMM_DEFAULT_P 1024 #define DGEMM_DEFAULT_P 512 -#define CGEMM_DEFAULT_P 256 -#define ZGEMM_DEFAULT_P 128 +#define CGEMM_DEFAULT_P 512 +#define ZGEMM_DEFAULT_P 256 #define SGEMM_DEFAULT_Q 1024 #define DGEMM_DEFAULT_Q 512 From f54ee8a94efca83e4e8c0811dd069b0696a620a8 Mon Sep 17 00:00:00 2001 From: michael tesch Date: Wed, 21 Oct 2026 23:43:41 -0500 Subject: [PATCH 2/2] arm64: complex GEMM on the level-3 driver and SME thread cap for VORTEXM4 On VORTEXM4, CGEMM/ZGEMM now skip the SME1 direct kernels and run on the level-3 driver with the SME2 complex kernel, which gives them the driver's threads and removes the static packing buffers of the SME1 kernels from this path. ARMV9SME keeps the SME1 path. The one-thread-per-SME-unit cap of num_cpu_avail() now covers complex level-3 BLAS and LAPACK too. --- common_thread.h | 4 ++-- interface/gemm.c | 5 +++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/common_thread.h b/common_thread.h index 8317d07941..2adec94e14 100644 --- a/common_thread.h +++ b/common_thread.h @@ -157,12 +157,12 @@ extern int blas_omp_threads_local; extern int blas_is_num_threads_set_explicitly; #if defined(ARCH_ARM64) && defined(OS_DARWIN) && defined(__clang__) && (defined(VORTEXM4) || defined(DYNAMIC_ARCH)) && \ - !defined(COMPLEX) && !defined(BFLOAT16) && !defined(HFLOAT16) + !defined(BFLOAT16) && !defined(HFLOAT16) #include #ifdef DYNAMIC_ARCH extern char *gotoblas_corename(void); #endif -/* Apple M4: real level-3 BLAS and LAPACK run on one SME unit per P-cluster; measured on M4 Pro only, not on M5+. */ +/* Apple M4: level-3 BLAS and LAPACK run on one SME unit per P-cluster; measured on M4 Pro only, not on M5+. */ static __inline int sme_level3_threads(int nthreads) { static int units = 0; #ifdef DYNAMIC_ARCH diff --git a/interface/gemm.c b/interface/gemm.c index b2e1815911..401fe6ffd3 100644 --- a/interface/gemm.c +++ b/interface/gemm.c @@ -581,10 +581,11 @@ else if (order == CblasRowMajor && k==lda && n==ldb && n==ldc && TransA == Cblas #endif //defined CBLAS #if !defined(BFLOAT16) && !defined(HFLOAT16) -#if defined(ARCH_ARM64) && (defined(USE_SGEMM_KERNEL_DIRECT)||defined(DYNAMIC_ARCH)) +/* complex GEMM on VORTEXM4 takes the level-3 driver with its SME2 kernels */ +#if defined(ARCH_ARM64) && (defined(USE_SGEMM_KERNEL_DIRECT)||defined(DYNAMIC_ARCH)) && !(defined(COMPLEX) && defined(VORTEXM4)) #if defined(DYNAMIC_ARCH) if (strcmp(gotoblas_corename(), "armv9sme") == 0 -#if defined(__clang__) +#if defined(__clang__) && !defined(COMPLEX) || strcmp(gotoblas_corename(), "vortexm4") == 0 #endif )