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/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 ) 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