#include <x86intrin.h>
#define PB_PT 1
/* lane m120 -- SI1: single induction register in the mikro k-loop.
base = lane_m111/arms/csp2s2pf.cpp (row artifact, sid 113912, 38.385 ms). */
#define FUSE67 1
#define GC_NT 0
#define MC 6
#define NC 1024
#define KC 512
/* m33: the pad-band cell. The record's board table has +6 at 38.563 and +12 at 38.568
-- 0.005 ms apart, INSIDE the 0.023 % board channel -- measured on an OLDER artifact.
This is the same cell on the BOARDED binary (tail2, which ships KCP = KC+6).
KCP is the A-panel row stride; KCP > KC is required (packA writes k <= KC per row). */
#define KCP (KC+12)
#define SCUT 512
// mmmd (double) : C = A*B -- Strassen-1 + packed AVX2/FMA 6x8 GEMM.
// Judge: g++-9 -static -O2 (no -march) => per-function target pragma.
#include <immintrin.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
typedef double T;
#ifndef MR
#define MR 6
#endif
#ifndef NR
#define NR 8
#endif
#ifndef KC
#define KC 192
#endif
#ifndef MC
#define MC 192
#endif
#ifndef NC
#define NC 1024
#endif
#ifndef SCUT
#define SCUT 640
#endif
#pragma GCC optimize("O3")
#pragma GCC push_options
#pragma GCC target("avx2,fma")
// ★ Under the judge's plain `g++-9 -O2` (no -march) GCC defaults to
// -mavx256-split-unaligned-load/-store: EVERY _mm256_loadu_*/_mm256_storeu_*
// becomes TWO instructions (xmm load + vinserti128 / vextracti128 + xmm store),
// and the insert/extract is a port-5 uop. `#pragma GCC target("avx2")` does NOT
// turn that off. Inline asm bypasses the lowering: one instruction, no port 5.
// AT&T operand order is src,dst for BOTH forms -> "vmovupd %1, %0" in each.
#define AVX2T __attribute__((target("avx2,fma")))
AVX2T static inline __m256d LD(const double *p) {
__m256d r; __asm__("vmovupd %1, %0" : "=x"(r) : "m"(*(const __m256d *)p)); return r;
}
AVX2T static inline void ST(double *p, __m256d v) {
__asm__ volatile("vmovupd %1, %0" : "=m"(*(__m256d *)p) : "x"(v));
}
// ---- E8NT: non-temporal stores for PURE (non-accumulating) stores. A pure
// store to a cold line pays a read-for-ownership of the whole line; vmovntpd
// skips it. Legal only when the destination is write-only AND 32B aligned.
AVX2T static inline void STNT(double *p, __m256d v) {
__asm__ volatile("vmovntpd %1, %0" : "=m"(*(__m256d *)p) : "x"(v));
}
// Packed A: panel p holds rows [p*MR, p*MR+MR): dst[p*MR*KC + r*KC + k] = S(row r, k)
// Packed B: panel g holds cols [g*NR, g*NR+NR): dst[g*KC*NR + k*NR + j] = S(k, col j)
static T g_par[(size_t)((MC + MR - 1) / MR) * MR * KCP + 64 + 448] __attribute__((aligned(64)));
static T g_pbr[(size_t)((NC + NR - 1) / NR) * NR * KC + 64 + 64] __attribute__((aligned(64)));
#define Apack (g_par + 265)
#define Bpack (g_pbr + 488)
// 8-wide op: dst[0..7] = s1 or s1+s2 or s1-s2
static inline void vop8(T *d, const T *s1, const T *s2, int op) {
if (op == 0) {
ST(d, LD(s1));
ST(d + 4, LD(s1 + 4));
} else if (op == 1) {
ST(d, _mm256_add_pd(LD(s1), LD(s2)));
ST(d + 4, _mm256_add_pd(LD(s1 + 4), LD(s2 + 4)));
} else {
ST(d, _mm256_sub_pd(LD(s1), LD(s2)));
ST(d + 4, _mm256_sub_pd(LD(s1 + 4), LD(s2 + 4)));
}
}
static void packA(int mc, int kc, const T *a1, const T *a2, int op, int lda, T *Ap) {
for (int p = 0; p < mc; p += MR) {
int mr = mc - p < MR ? mc - p : MR;
T *dst = Ap + (size_t)(p / MR) * MR * KCP;
for (int r = 0; r < mr; r++) {
const T *s1 = a1 + (size_t)(p + r) * lda;
/* E8PF: packA reads kc contiguous doubles per row but the ROWS are lda*8 = 16 KB
apart, so each row START is a demand miss the hw prefetcher cannot see. Start
the next row's stream early. (This is also the fix packB's k-outer swap gave.) */
if (r + 1 < mr) __builtin_prefetch(s1 + lda, 0, 1);
else if (p + MR + r < mc) __builtin_prefetch(s1 + (size_t)MR * lda, 0, 1);
const T *s2 = a2 ? a2 + (size_t)(p + r) * lda : (const T *)0;
/* lane m111: the E8PF row-start prefetch above has NEVER had its s2 twin. 5 of the
7 leaf products pass a second source (a2 != 0: op +1/-1) and read it with vop8 at
the SAME lda*8 = 16 KB row stride, so every one of its row starts is the same
hw-unpredictable demand miss that the s1 prefetch exists to hide. Mirror it. */
if (s2) {
if (r + 1 < mr) __builtin_prefetch(s2 + lda, 0, 1);
else if (p + MR + r < mc) __builtin_prefetch(s2 + (size_t)MR * lda, 0, 1);
}
T *d = dst + (size_t)r * KCP;
int k = 0;
if (op == 0) { for (; k + 8 <= kc; k += 8) vop8(d + k, s1 + k, 0, 0); }
else if (op == 1) { for (; k + 8 <= kc; k += 8) vop8(d + k, s1 + k, s2 + k, 1); }
else { for (; k + 8 <= kc; k += 8) vop8(d + k, s1 + k, s2 + k, 2); }
for (; k < kc; k++) d[k] = op == 0 ? s1[k] : (op == 1 ? s1[k] + s2[k] : s1[k] - s2[k]);
}
}
}
static void packB(int kc, int nc, const T *b1, const T *b2, int op, int lda, T *Bp) {
/* E8SWAP: k OUTER, g INNER -> the source read b1[k*lda + g] is CONTIGUOUS in g across the
whole row, so it becomes a hardware-prefetchable stream instead of kc lines at a
lda*8 = 16 KB stride (8 doubles = one full line per 16 KB jump = a demand miss each). */
#ifdef PB_PT
/* ONE dense contiguous sweep of the destination span, before the strided writes. No
counter, no timer, no other change: this file differs from the shipped artifact ONLY
by this block, so the board A/B isolates it exactly. */
{ static int pb_warmed = 0;
if (!pb_warmed) { pb_warmed = 1;
size_t _sp = (size_t)((nc + NR - 1) / NR) * KC * NR;
memset((void *)Bp, 0, _sp * sizeof(T)); } }
#endif
for (int k0 = 0; k0 < kc; k0 += 8)
for (int g = 0; g + NR <= nc; g += NR)
for (int k = k0; k < k0 + 8 && k < kc; k++)
vop8(Bp + (size_t)(g / NR) * KC * NR + (size_t)k * NR,
b1 + (size_t)k * lda + g, b2 ? b2 + (size_t)k * lda + g : (const T *)0, op);
for (int g = nc - (nc % NR); g < nc; g++) {
T *d = Bp + (size_t)(g / NR) * KC * NR;
for (int k = 0; k < kc; k++) {
T *dd = d + (size_t)k * NR;
for (int j = 0; j < NR; j++) dd[j] = (g + j < nc) ? (op == 0 ? b1[(size_t)k*lda+g+j] : (op == 1 ? b1[(size_t)k*lda+g+j]+b2[(size_t)k*lda+g+j] : b1[(size_t)k*lda+g+j]-b2[(size_t)k*lda+g+j])) : 0.0;
}
}
}
static void packB_orig(int kc, int nc, const T *b1, const T *b2, int op, int lda, T *Bp) {
for (int g = 0; g < nc; g += NR) {
int nr = nc - g < NR ? nc - g : NR;
T *dst = Bp + (size_t)(g / NR) * KC * NR;
for (int k = 0; k < kc; k++) {
T *d = dst + (size_t)k * NR;
const T *s1 = b1 + (size_t)k * lda + g;
const T *s2 = b2 ? b2 + (size_t)k * lda + g : (const T *)0;
if (nr == NR) { vop8(d, s1, s2, op); }
else {
for (int j = 0; j < nr; j++) d[j] = op == 0 ? s1[j] : (op == 1 ? s1[j] + s2[j] : s1[j] - s2[j]);
for (int j = nr; j < NR; j++) d[j] = 0.0;
}
}
}
}
#define KSTEP(U) \
do { \
__m256d b0 = LD(bp + (size_t)(k + (U)) * NR); \
__m256d b1 = LD(bp + (size_t)(k + (U)) * NR + 4); \
__m256d t; \
t = _mm256_broadcast_sd(ap + 0 * KCP + (k + (U))); \
c00 = _mm256_fmadd_pd(t, b0, c00); c01 = _mm256_fmadd_pd(t, b1, c01); \
t = _mm256_broadcast_sd(ap + 1 * KCP + (k + (U))); \
c10 = _mm256_fmadd_pd(t, b0, c10); c11 = _mm256_fmadd_pd(t, b1, c11); \
t = _mm256_broadcast_sd(ap + 2 * KCP + (k + (U))); \
c20 = _mm256_fmadd_pd(t, b0, c20); c21 = _mm256_fmadd_pd(t, b1, c21); \
t = _mm256_broadcast_sd(ap + 3 * KCP + (k + (U))); \
c30 = _mm256_fmadd_pd(t, b0, c30); c31 = _mm256_fmadd_pd(t, b1, c31); \
t = _mm256_broadcast_sd(ap + 4 * KCP + (k + (U))); \
c40 = _mm256_fmadd_pd(t, b0, c40); c41 = _mm256_fmadd_pd(t, b1, c41); \
t = _mm256_broadcast_sd(ap + 5 * KCP + (k + (U))); \
c50 = _mm256_fmadd_pd(t, b0, c50); c51 = _mm256_fmadd_pd(t, b1, c51); \
} while (0)
// ncols guard: B panel is zero-padded past ncols, so only the store needs masking.
__attribute__((always_inline)) static inline void micro(int kc, const T *__restrict ap, const T *__restrict bp,
T *__restrict C, int ldc, int accumulate, int mrows, int ncols,
int nt, T *__restrict C2, int ldc2, int acc2, int s2, int nt2) {
__m256d c00, c01, c10, c11, c20, c21, c30, c31, c40, c41, c50, c51;
c00 = c01 = c10 = c11 = c20 = c21 = c30 = c31 = c40 = c41 = c50 = c51 = _mm256_setzero_pd();
/* lane m106: port of lane m78's BOARD-PROVEN `csp_acc` (mmmd2k, -1.070 %): warm the six C
destination lines ONLY on the accumulating (RMW) path, at the TOP of the call so their
latency hides under the k-loop. With `nt` set the pure stores are vmovntpd -- they bypass
L1, so prefetching the destination there fetches a line that is overwritten without ever
being read. Per-row `mrows > r` guards keep the tail in bounds. */
if (accumulate) {
if (mrows > 0) __builtin_prefetch(C + 0 * (size_t)ldc, 1, 3);
if (mrows > 1) __builtin_prefetch(C + 1 * (size_t)ldc, 1, 3);
if (mrows > 2) __builtin_prefetch(C + 2 * (size_t)ldc, 1, 3);
if (mrows > 3) __builtin_prefetch(C + 3 * (size_t)ldc, 1, 3);
if (mrows > 4) __builtin_prefetch(C + 4 * (size_t)ldc, 1, 3);
if (mrows > 5) __builtin_prefetch(C + 5 * (size_t)ldc, 1, 3);
}
/* lane m106: the SAME guard applied to the SECOND destination -- this is m14's `h1_hoistc2`
INTENT (the six cold C2 lines demanded early so their latency hides under the k-loop)
done the way that cannot corrupt anything: a prefetch, not a register. `ST2`'s RMW reads
only fire for acc2 != 0; for acc2 == 0 and nt2 the store is vmovntpd and warming it would
fetch a line that is never read. `acc2 < 0` is also an RMW (subtract). */
if (C2 != 0 && acc2 != 0) {
if (mrows > 0) __builtin_prefetch(C2 + 0 * (size_t)ldc2, 1, 3);
if (mrows > 1) __builtin_prefetch(C2 + 1 * (size_t)ldc2, 1, 3);
if (mrows > 2) __builtin_prefetch(C2 + 2 * (size_t)ldc2, 1, 3);
if (mrows > 3) __builtin_prefetch(C2 + 3 * (size_t)ldc2, 1, 3);
if (mrows > 4) __builtin_prefetch(C2 + 4 * (size_t)ldc2, 1, 3);
if (mrows > 5) __builtin_prefetch(C2 + 5 * (size_t)ldc2, 1, 3);
}
/* E8SPF: the store block at the END touches 6 lines at an ldc*8 = 16 KB stride,
each a demand miss + RFO, measured at 96 cyc/mikro-call. Start them early. */
#define KSTEPQ(U) \
do { \
const char *bb = bpc + q*8 + (size_t)(U) * NR * 8; \
const char *aa = apc + q + (size_t)(U) * 8; \
__m256d b0 = LD((const T *)(bb)); \
__m256d b1 = LD((const T *)(bb + 32)); \
__m256d t; \
t = _mm256_broadcast_sd((const T *)(aa + 0 * KCP * 8)); \
c00 = _mm256_fmadd_pd(t, b0, c00); c01 = _mm256_fmadd_pd(t, b1, c01); \
t = _mm256_broadcast_sd((const T *)(aa + 1 * KCP * 8)); \
c10 = _mm256_fmadd_pd(t, b0, c10); c11 = _mm256_fmadd_pd(t, b1, c11); \
t = _mm256_broadcast_sd((const T *)(aa + 2 * KCP * 8)); \
c20 = _mm256_fmadd_pd(t, b0, c20); c21 = _mm256_fmadd_pd(t, b1, c21); \
t = _mm256_broadcast_sd((const T *)(aa + 3 * KCP * 8)); \
c30 = _mm256_fmadd_pd(t, b0, c30); c31 = _mm256_fmadd_pd(t, b1, c31); \
t = _mm256_broadcast_sd((const T *)(aa + 4 * KCP * 8)); \
c40 = _mm256_fmadd_pd(t, b0, c40); c41 = _mm256_fmadd_pd(t, b1, c41); \
t = _mm256_broadcast_sd((const T *)(aa + 5 * KCP * 8)); \
c50 = _mm256_fmadd_pd(t, b0, c50); c51 = _mm256_fmadd_pd(t, b1, c51); \
} while (0)
{
const char *apc = (const char *)(ap + kc);
const char *bpc = (const char *)(bp + (size_t)kc * NR);
/* lane m120 SI1 -- see rig/gen_si1.py. `rem` floors the panel run to a
multiple of 32 B so the asm loop's single `addq/jne` tail terminates
exactly; the C remainder loop below keeps the TRUE bases and `rem` steps. */
const int si_rem = kc & 3;
long qm = -(long)(kc - si_rem) * 8;
if (qm < 0) {
const char *apcR = (const char *)(ap + (kc - si_rem));
const char *bpcR = (const char *)(bp + (size_t)(kc - si_rem) * NR);
/* ---- lane m120 SI1: ONE induction register for BOTH panels ----
A advances 32 B/body, B advances 256 B/body (ratio exactly 8), so with the
index in %[q] the two address streams are
A: disp(%[pa],%[q],1) B: disp(%[pb],%[q],8)
and the loop tail is `addq $32,%[q] ; jne` -- 1 macro-fused uop instead of
the 3 uops (`2 addq + cmpq + jne`) that strength-reducing BOTH pointers
costs. %[pa]/%[pb] are the panels' bases shifted by -8*rem / -64*rem so
the index starts on a multiple of 32 and lands exactly on 0. */
{
__m256d t0, t1, t2, t3;
t0 = t1 = t2 = t3 = _mm256_setzero_pd(); /* keep all 16 ymm live so gcc must
give the 16 operands DISTINCT regs */
{
/* ---- SPLIT at 7/10: the six C2 prefetches move from the TOP of the
call to the top of the SECOND k-block, inside the L1 survival window. ---- */
const long L = (long)(kc - si_rem);
long kc1 = (long)(((long long)L * 7) / 10);
kc1 &= ~3L; if (kc1 > L) kc1 = L; if (kc1 < 0) kc1 = 0;
const long qA0 = -kc1 * 8;
const long qB0 = -(L - kc1) * 8;
const char *paA = (const char*)(ap + kc1);
const char *pbA = (const char*)(bp + (size_t)kc1 * NR);
const char *c2b0 = (C2 != 0 && acc2 != 0) ? (const char *)C2 : (const char *)C;
char *c2b = (char *)c2b0; const long c2s = (long)ldc2 * 8;
{ __m256d t0, t1, t2, t3; t0=t1=t2=t3=_mm256_setzero_pd();
__asm__ volatile (
"vxorpd %[c00], %[c00], %[c00]\n\t"
"vxorpd %[c01], %[c01], %[c01]\n\t"
"vxorpd %[c10], %[c10], %[c10]\n\t"
"vxorpd %[c11], %[c11], %[c11]\n\t"
"vxorpd %[c20], %[c20], %[c20]\n\t"
"vxorpd %[c21], %[c21], %[c21]\n\t"
"vxorpd %[c30], %[c30], %[c30]\n\t"
"vxorpd %[c31], %[c31], %[c31]\n\t"
"vxorpd %[c40], %[c40], %[c40]\n\t"
"vxorpd %[c41], %[c41], %[c41]\n\t"
"vxorpd %[c50], %[c50], %[c50]\n\t"
"vxorpd %[c51], %[c51], %[c51]\n\t"
"cmpq $0, %[qA0]\n\t"
"je .LsA%= \n\t"
"movq %[qA0], %[q]\n\t"
".LsiA%=:\n\t"
"prefetcht0 1024(%[pba],%[q],8)\n\t"
"prefetcht0 1088(%[pba],%[q],8)\n\t"
"prefetcht0 1152(%[pba],%[q],8)\n\t"
"prefetcht0 1216(%[pba],%[q],8)\n\t"
"/* ---- k+0 ---- */\n\t"
"vmovupd 0(%[pba],%[q],8), %[t2]\n\t"
"vmovupd 32(%[pba],%[q],8), %[t3]\n\t"
"vbroadcastsd 0(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c00]\n\t"
"vfmadd231pd %[t3], %[t0], %[c01]\n\t"
"vbroadcastsd 4192(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c10]\n\t"
"vfmadd231pd %[t3], %[t0], %[c11]\n\t"
"vbroadcastsd 8384(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c20]\n\t"
"vfmadd231pd %[t3], %[t0], %[c21]\n\t"
"vbroadcastsd 12576(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c30]\n\t"
"vfmadd231pd %[t3], %[t0], %[c31]\n\t"
"vbroadcastsd 16768(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c40]\n\t"
"vfmadd231pd %[t3], %[t0], %[c41]\n\t"
"vbroadcastsd 20960(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c50]\n\t"
"vfmadd231pd %[t3], %[t0], %[c51]\n\t"
"/* ---- k+1 ---- */\n\t"
"vmovupd 64(%[pba],%[q],8), %[t2]\n\t"
"vmovupd 96(%[pba],%[q],8), %[t3]\n\t"
"vbroadcastsd 8(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c00]\n\t"
"vfmadd231pd %[t3], %[t0], %[c01]\n\t"
"vbroadcastsd 4200(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c10]\n\t"
"vfmadd231pd %[t3], %[t0], %[c11]\n\t"
"vbroadcastsd 8392(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c20]\n\t"
"vfmadd231pd %[t3], %[t0], %[c21]\n\t"
"vbroadcastsd 12584(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c30]\n\t"
"vfmadd231pd %[t3], %[t0], %[c31]\n\t"
"vbroadcastsd 16776(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c40]\n\t"
"vfmadd231pd %[t3], %[t0], %[c41]\n\t"
"vbroadcastsd 20968(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c50]\n\t"
"vfmadd231pd %[t3], %[t0], %[c51]\n\t"
"/* ---- k+2 ---- */\n\t"
"vmovupd 128(%[pba],%[q],8), %[t2]\n\t"
"vmovupd 160(%[pba],%[q],8), %[t3]\n\t"
"vbroadcastsd 16(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c00]\n\t"
"vfmadd231pd %[t3], %[t0], %[c01]\n\t"
"vbroadcastsd 4208(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c10]\n\t"
"vfmadd231pd %[t3], %[t0], %[c11]\n\t"
"vbroadcastsd 8400(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c20]\n\t"
"vfmadd231pd %[t3], %[t0], %[c21]\n\t"
"vbroadcastsd 12592(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c30]\n\t"
"vfmadd231pd %[t3], %[t0], %[c31]\n\t"
"vbroadcastsd 16784(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c40]\n\t"
"vfmadd231pd %[t3], %[t0], %[c41]\n\t"
"vbroadcastsd 20976(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c50]\n\t"
"vfmadd231pd %[t3], %[t0], %[c51]\n\t"
"/* ---- k+3 ---- */\n\t"
"vmovupd 192(%[pba],%[q],8), %[t2]\n\t"
"vmovupd 224(%[pba],%[q],8), %[t3]\n\t"
"vbroadcastsd 24(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c00]\n\t"
"vfmadd231pd %[t3], %[t0], %[c01]\n\t"
"vbroadcastsd 4216(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c10]\n\t"
"vfmadd231pd %[t3], %[t0], %[c11]\n\t"
"vbroadcastsd 8408(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c20]\n\t"
"vfmadd231pd %[t3], %[t0], %[c21]\n\t"
"vbroadcastsd 12600(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c30]\n\t"
"vfmadd231pd %[t3], %[t0], %[c31]\n\t"
"vbroadcastsd 16792(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c40]\n\t"
"vfmadd231pd %[t3], %[t0], %[c41]\n\t"
"vbroadcastsd 20984(%[paa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c50]\n\t"
"vfmadd231pd %[t3], %[t0], %[c51]\n\t"
"addq $32, %[q]\n\t"
"jne .LsiA%=\n\t"
".LsA%=:\n\t"
"prefetcht0 (%[c2b])\n\t"
"addq %[c2s], %[c2b]\n\t"
"prefetcht0 (%[c2b])\n\t"
"addq %[c2s], %[c2b]\n\t"
"prefetcht0 (%[c2b])\n\t"
"addq %[c2s], %[c2b]\n\t"
"prefetcht0 (%[c2b])\n\t"
"addq %[c2s], %[c2b]\n\t"
"prefetcht0 (%[c2b])\n\t"
"addq %[c2s], %[c2b]\n\t"
"prefetcht0 (%[c2b])\n\t"
"addq %[c2s], %[c2b]\n\t"
"cmpq $0, %[qB0]\n\t"
"je .LsB%= \n\t"
"movq %[qB0], %[q]\n\t"
".LsiB%=:\n\t"
"prefetcht0 1024(%[pb],%[q],8)\n\t"
"prefetcht0 1088(%[pb],%[q],8)\n\t"
"prefetcht0 1152(%[pb],%[q],8)\n\t"
"prefetcht0 1216(%[pb],%[q],8)\n\t"
"/* ---- k+0 ---- */\n\t"
"vmovupd 0(%[pb],%[q],8), %[t2]\n\t"
"vmovupd 32(%[pb],%[q],8), %[t3]\n\t"
"vbroadcastsd 0(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c00]\n\t"
"vfmadd231pd %[t3], %[t0], %[c01]\n\t"
"vbroadcastsd 4192(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c10]\n\t"
"vfmadd231pd %[t3], %[t0], %[c11]\n\t"
"vbroadcastsd 8384(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c20]\n\t"
"vfmadd231pd %[t3], %[t0], %[c21]\n\t"
"vbroadcastsd 12576(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c30]\n\t"
"vfmadd231pd %[t3], %[t0], %[c31]\n\t"
"vbroadcastsd 16768(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c40]\n\t"
"vfmadd231pd %[t3], %[t0], %[c41]\n\t"
"vbroadcastsd 20960(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c50]\n\t"
"vfmadd231pd %[t3], %[t0], %[c51]\n\t"
"/* ---- k+1 ---- */\n\t"
"vmovupd 64(%[pb],%[q],8), %[t2]\n\t"
"vmovupd 96(%[pb],%[q],8), %[t3]\n\t"
"vbroadcastsd 8(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c00]\n\t"
"vfmadd231pd %[t3], %[t0], %[c01]\n\t"
"vbroadcastsd 4200(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c10]\n\t"
"vfmadd231pd %[t3], %[t0], %[c11]\n\t"
"vbroadcastsd 8392(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c20]\n\t"
"vfmadd231pd %[t3], %[t0], %[c21]\n\t"
"vbroadcastsd 12584(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c30]\n\t"
"vfmadd231pd %[t3], %[t0], %[c31]\n\t"
"vbroadcastsd 16776(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c40]\n\t"
"vfmadd231pd %[t3], %[t0], %[c41]\n\t"
"vbroadcastsd 20968(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c50]\n\t"
"vfmadd231pd %[t3], %[t0], %[c51]\n\t"
"/* ---- k+2 ---- */\n\t"
"vmovupd 128(%[pb],%[q],8), %[t2]\n\t"
"vmovupd 160(%[pb],%[q],8), %[t3]\n\t"
"vbroadcastsd 16(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c00]\n\t"
"vfmadd231pd %[t3], %[t0], %[c01]\n\t"
"vbroadcastsd 4208(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c10]\n\t"
"vfmadd231pd %[t3], %[t0], %[c11]\n\t"
"vbroadcastsd 8400(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c20]\n\t"
"vfmadd231pd %[t3], %[t0], %[c21]\n\t"
"vbroadcastsd 12592(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c30]\n\t"
"vfmadd231pd %[t3], %[t0], %[c31]\n\t"
"vbroadcastsd 16784(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c40]\n\t"
"vfmadd231pd %[t3], %[t0], %[c41]\n\t"
"vbroadcastsd 20976(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c50]\n\t"
"vfmadd231pd %[t3], %[t0], %[c51]\n\t"
"/* ---- k+3 ---- */\n\t"
"vmovupd 192(%[pb],%[q],8), %[t2]\n\t"
"vmovupd 224(%[pb],%[q],8), %[t3]\n\t"
"vbroadcastsd 24(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c00]\n\t"
"vfmadd231pd %[t3], %[t0], %[c01]\n\t"
"vbroadcastsd 4216(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c10]\n\t"
"vfmadd231pd %[t3], %[t0], %[c11]\n\t"
"vbroadcastsd 8408(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c20]\n\t"
"vfmadd231pd %[t3], %[t0], %[c21]\n\t"
"vbroadcastsd 12600(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c30]\n\t"
"vfmadd231pd %[t3], %[t0], %[c31]\n\t"
"vbroadcastsd 16792(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c40]\n\t"
"vfmadd231pd %[t3], %[t0], %[c41]\n\t"
"vbroadcastsd 20984(%[pa],%[q],1), %[t0]\n\t"
"vfmadd231pd %[t2], %[t0], %[c50]\n\t"
"vfmadd231pd %[t3], %[t0], %[c51]\n\t"
"addq $32, %[q]\n\t"
"jne .LsiB%=\n\t"
".LsB%=:\n\t"
: [c00]"=x"(c00), [c01]"=x"(c01), [c10]"=x"(c10), [c11]"=x"(c11),
[c20]"=x"(c20), [c21]"=x"(c21), [c30]"=x"(c30), [c31]"=x"(c31),
[c40]"=x"(c40), [c41]"=x"(c41), [c50]"=x"(c50), [c51]"=x"(c51),
[t0]"+x"(t0), [t2]"+x"(t2), [t3]"+x"(t3),
[q]"+&r"(qm), [c2b]"+&r"(c2b)
: [pa]"r"(apcR), [pb]"r"(bpcR), [paa]"r"(paA), [pba]"r"(pbA),
[qA0]"r"(qA0), [qB0]"r"(qB0), [c2s]"r"(c2s)
: "memory", "cc");
(void)t0;(void)t1;(void)t2;(void)t3; }
qm = qB0;
}
(void)t0;(void)t1;(void)t2;(void)t3;
}
}
long q = -(long)si_rem * 8;
for (; q < 0; q += 8) KSTEPQ(0);
}
#undef KSTEPQ
if (ncols == NR && mrows == MR) {
#define STO(r, a, b) \
{ \
T *p = C + (size_t)(r) * ldc; \
if (accumulate) { \
ST(p, _mm256_add_pd(LD(p), a)); \
ST(p + 4, _mm256_add_pd(LD(p + 4), b)); \
} else { if (nt) { STNT(p, a); STNT(p + 4, b); } \
else { ST(p, a); ST(p + 4, b); } } \
}
STO(0, c00, c01) STO(1, c10, c11) STO(2, c20, c21)
STO(3, c30, c31) STO(4, c40, c41) STO(5, c50, c51)
#undef STO
if (C2) {
#define ST2(r, a, b) \
{ \
T *q = C2 + (size_t)(r) * ldc2; \
if (acc2 > 0) { ST(q, _mm256_add_pd(LD(q), a)); ST(q + 4, _mm256_add_pd(LD(q + 4), b)); } \
else if (acc2 < 0) { ST(q, _mm256_sub_pd(LD(q), a)); ST(q + 4, _mm256_sub_pd(LD(q + 4), b)); } \
else if (nt2) { STNT(q, a); STNT(q + 4, b); } \
else { ST(q, a); ST(q + 4, b); } \
}
ST2(0, c00, c01) ST2(1, c10, c11) ST2(2, c20, c21)
ST2(3, c30, c31) ST2(4, c40, c41) ST2(5, c50, c51)
#undef ST2
}
} else {
__m256i ml = _mm256_set_epi64x(ncols > 3 ? -1LL : 0, ncols > 2 ? -1LL : 0,
ncols > 1 ? -1LL : 0, ncols > 0 ? -1LL : 0);
__m256i mh = _mm256_set_epi64x(ncols > 7 ? -1LL : 0, ncols > 6 ? -1LL : 0,
ncols > 5 ? -1LL : 0, ncols > 4 ? -1LL : 0);
#define STO2(r, a, b) \
if (r < mrows) { \
T *p = C + (size_t)(r) * ldc; \
if (accumulate) { \
_mm256_maskstore_pd(p, ml, _mm256_add_pd(_mm256_maskload_pd(p, ml), a)); \
_mm256_maskstore_pd(p + 4, mh, _mm256_add_pd(_mm256_maskload_pd(p + 4, mh), b)); \
} else { _mm256_maskstore_pd(p, ml, a); _mm256_maskstore_pd(p + 4, mh, b); } \
}
STO2(0, c00, c01) STO2(1, c10, c11) STO2(2, c20, c21)
STO2(3, c30, c31) STO2(4, c40, c41) STO2(5, c50, c51)
#undef STO2
if (C2) {
#define ST2M(r, a, b) \
if (r < mrows) { \
T *q = C2 + (size_t)(r) * ldc2; \
if (acc2 > 0) { _mm256_maskstore_pd(q, ml, _mm256_add_pd(_mm256_maskload_pd(q, ml), a)); \
_mm256_maskstore_pd(q + 4, mh, _mm256_add_pd(_mm256_maskload_pd(q + 4, mh), b)); } \
else if (acc2 < 0) { _mm256_maskstore_pd(q, ml, _mm256_sub_pd(_mm256_maskload_pd(q, ml), a)); \
_mm256_maskstore_pd(q + 4, mh, _mm256_sub_pd(_mm256_maskload_pd(q + 4, mh), b)); } \
else { _mm256_maskstore_pd(q, ml, a); _mm256_maskstore_pd(q + 4, mh, b); } \
}
ST2M(0, c00, c01) ST2M(1, c10, c11) ST2M(2, c20, c21)
ST2M(3, c30, c31) ST2M(4, c40, c41) ST2M(5, c50, c51)
#undef ST2M
}
}
}
/* E8FUSE: the two combine axpys of one Strassen product each read the WHOLE m x m scratch
M; doing them in ONE pass reads M once instead of twice (5 products x one 8 MB read saved).
Bit-identical: the same operations on the same values, only the visit order changes. */
static void blk_axpy2(int m, int n, T *dA, int ldA, T sA, int accA,
T *dB, int ldB, T sB, int accB,
const T *src, int lds) {
const int ntA = (!accA) && ((((unsigned long)dA) & 31) == 0) && ((ldA & 3) == 0);
const int ntB = (!accB) && ((((unsigned long)dB) & 31) == 0) && ((ldB & 3) == 0);
for (int i = 0; i < m; i++) {
T *a = dA + (size_t)i * ldA, *b = dB + (size_t)i * ldB;
const T *x = src + (size_t)i * lds;
/* E8RPF: next row starts are ldA*8 = ldc*8 = 16 KB apart -> demand misses */
if (i + 1 < m) {
if (i + 2 < m) {
}
}
int j = 0;
for (; j + 4 <= n; j += 4) {
__m256d v = LD(x + j);
if (accA) { if (sA > 0) ST(a + j, _mm256_add_pd(LD(a + j), v));
else ST(a + j, _mm256_sub_pd(LD(a + j), v)); }
else { if (!ntA) { if (sA > 0) ST(a + j, v);
else ST(a + j, _mm256_sub_pd(_mm256_setzero_pd(), v)); }
else { if (sA > 0) STNT(a + j, v);
else STNT(a + j, _mm256_sub_pd(_mm256_setzero_pd(), v)); } }
if (accB) { if (sB > 0) ST(b + j, _mm256_add_pd(LD(b + j), v));
else ST(b + j, _mm256_sub_pd(LD(b + j), v)); }
else { if (!ntB) { if (sB > 0) ST(b + j, v);
else ST(b + j, _mm256_sub_pd(_mm256_setzero_pd(), v)); }
else { if (sB > 0) STNT(b + j, v);
else STNT(b + j, _mm256_sub_pd(_mm256_setzero_pd(), v)); } }
}
for (; j < n; j++) { T v = x[j];
a[j] = accA ? (sA > 0 ? a[j] + v : a[j] - v) : (sA > 0 ? v : -v);
b[j] = accB ? (sB > 0 ? b[j] + v : b[j] - v) : (sB > 0 ? v : -v); }
}
if (ntA || ntB) __asm__ volatile("sfence" ::: "memory");
}
/* lane m31 TAIL2: 512 = 85*6 + 2, so the last ic block of EVERY (jr, product) calls the
mikro with mrows=2 while the mikro computes all SIX rows (mrows only masks the stores).
1.835 M FMAs per call = 0.78 % of the FMA count, thrown away. This kernel does exactly
the two real rows: 4 FMAs/k instead of 12, same B loads, same C/C2 semantics. */
__attribute__((always_inline)) static inline void micro2(int kc, const T *__restrict ap, const T *__restrict bp,
T *__restrict C, int ldc, int accumulate,
int nt, T *__restrict C2, int ldc2, int acc2, int nt2) {
__m256d c00, c01, c10, c11;
c00 = c01 = c10 = c11 = _mm256_setzero_pd();
#define KSTEP2(U) \
do { \
const char *bb = bpc + q*8 + (size_t)(U) * NR * 8; \
const char *aa = apc + q + (size_t)(U) * 8; \
__m256d b0 = LD((const T *)(bb)); \
__m256d b1 = LD((const T *)(bb + 32)); \
__m256d t; \
t = _mm256_broadcast_sd((const T *)(aa + 0 * KCP * 8)); \
c00 = _mm256_fmadd_pd(t, b0, c00); c01 = _mm256_fmadd_pd(t, b1, c01); \
t = _mm256_broadcast_sd((const T *)(aa + 1 * KCP * 8)); \
c10 = _mm256_fmadd_pd(t, b0, c10); c11 = _mm256_fmadd_pd(t, b1, c11); \
} while (0)
{
const char *apc = (const char *)(ap + kc);
const char *bpc = (const char *)(bp + (size_t)kc * NR);
long q = -(long)kc * 8;
for (; q + 32 <= 0; q += 32) {
__builtin_prefetch(bpc + q*8 + (size_t)(16 + 0) * NR * 8);
__builtin_prefetch(bpc + q*8 + (size_t)(16 + 1) * NR * 8);
__builtin_prefetch(bpc + q*8 + (size_t)(16 + 2) * NR * 8);
__builtin_prefetch(bpc + q*8 + (size_t)(16 + 3) * NR * 8);
KSTEP2(0); KSTEP2(1); KSTEP2(2); KSTEP2(3); }
for (; q < 0; q += 8) KSTEP2(0);
}
#undef KSTEP2
{
T *p = C;
if (accumulate) { ST(p, _mm256_add_pd(LD(p), c00)); ST(p + 4, _mm256_add_pd(LD(p + 4), c01)); }
else if (nt) { STNT(p, c00); STNT(p + 4, c01); }
else { ST(p, c00); ST(p + 4, c01); }
p = C + (size_t)ldc;
if (accumulate) { ST(p, _mm256_add_pd(LD(p), c10)); ST(p + 4, _mm256_add_pd(LD(p + 4), c11)); }
else if (nt) { STNT(p, c10); STNT(p + 4, c11); }
else { ST(p, c10); ST(p + 4, c11); }
}
if (C2) {
T *q2 = C2;
if (acc2 > 0) { ST(q2, _mm256_add_pd(LD(q2), c00)); ST(q2 + 4, _mm256_add_pd(LD(q2 + 4), c01)); }
else if (acc2 < 0) { ST(q2, _mm256_sub_pd(LD(q2), c00)); ST(q2 + 4, _mm256_sub_pd(LD(q2 + 4), c01)); }
else if (nt2) { STNT(q2, c00); STNT(q2 + 4, c01); }
else { ST(q2, c00); ST(q2 + 4, c01); }
q2 = C2 + (size_t)ldc2;
if (acc2 > 0) { ST(q2, _mm256_add_pd(LD(q2), c10)); ST(q2 + 4, _mm256_add_pd(LD(q2 + 4), c11)); }
else if (acc2 < 0) { ST(q2, _mm256_sub_pd(LD(q2), c10)); ST(q2 + 4, _mm256_sub_pd(LD(q2 + 4), c11)); }
else if (nt2) { STNT(q2, c10); STNT(q2 + 4, c11); }
else { ST(q2, c10); ST(q2 + 4, c11); }
}
}
/* +++ lane m175 `micro2w` -- the 2-row TAIL with TWO B-groups (16 columns) per call. +++
MEASURED on the judge (rig/bench.cpp, this lane, custom_test, judge rdtsc with TSC==core
proved by REF ticks/FMA = 0.5002):
mikro 24.589 ticks/body (p01 floor 48/2 = 24.0) => +2.5 %
micro2 16.325 ticks/body (p01 floor 8.0, issue 10.0) => +63 %
16.32 == 4 chained FMAs per accumulator x the 4-cycle FMA latency EXACTLY, so the shipped
2-row tail is DEPENDENCY-CHAIN bound, not port bound: each of its four accumulators takes one
FMA per k-step and a 4-k-step body cannot hide a 4-cycle chain.
THE CONTROL THAT NAMES THE BINDER: `mp2` (a software pipeline of micro2's broadcast AND B
load, built and benched) measured 16.3244 against the shipped 16.3254 -- a NULL. So it is
not load latency. Doubling the COLUMN tile to 2*NR gives 8 accumulators; per k-step 8 FMAs
against a 4-cycle chain => 16 cycles per 32 FMAs = 2.0 FMA/cycle, the port bound.
REGISTER NOTE: 8 accumulators + 4 B + 1 broadcast = 13 of 16 ymm -- it fits, and that is only
true because the tile is 2 rows. The 6x8 `micro` cannot do this (12 acc + 4 = 16, exactly). */
__attribute__((always_inline)) static inline void micro2w(int kc, const T *__restrict ap,
const T *__restrict bp0, const T *__restrict bp1,
T *__restrict C, int ldc, int accumulate,
int nt, T *__restrict C2, int ldc2, int acc2, int nt2) {
__m256d c00, c01, c02, c03, c10, c11, c12, c13;
c00 = c01 = c02 = c03 = c10 = c11 = c12 = c13 = _mm256_setzero_pd();
if (accumulate) { __builtin_prefetch(C + 0 * (size_t)ldc, 1, 3);
__builtin_prefetch(C + 1 * (size_t)ldc, 1, 3); }
if (C2 != 0 && acc2 != 0) { __builtin_prefetch(C2 + 0 * (size_t)ldc2, 1, 3);
__builtin_prefetch(C2 + 1 * (size_t)ldc2, 1, 3); }
#define K2W(U) \
do { \
const char *_a = aa + (size_t)(U) * 8; \
const char *_b0 = bb0 + (size_t)(U) * NR * 8; \
const char *_b1 = bb1 + (size_t)(U) * NR * 8; \
__m256d b0 = LD((const T *)(_b0)); \
__m256d b1 = LD((const T *)(_b0 + 32)); \
__m256d b2 = LD((const T *)(_b1)); \
__m256d b3 = LD((const T *)(_b1 + 32)); \
__m256d t; \
t = _mm256_broadcast_sd((const T *)(_a)); \
c00 = _mm256_fmadd_pd(t, b0, c00); c01 = _mm256_fmadd_pd(t, b1, c01); \
c02 = _mm256_fmadd_pd(t, b2, c02); c03 = _mm256_fmadd_pd(t, b3, c03); \
t = _mm256_broadcast_sd((const T *)(_a + KCP * 8)); \
c10 = _mm256_fmadd_pd(t, b0, c10); c11 = _mm256_fmadd_pd(t, b1, c11); \
c12 = _mm256_fmadd_pd(t, b2, c12); c13 = _mm256_fmadd_pd(t, b3, c13); \
} while (0)
{
const char *aa = (const char *)ap;
const char *bb0 = (const char *)bp0;
const char *bb1 = (const char *)bp1;
int nb = kc >> 2, rem = kc & 3;
for (int b = 0; b < nb; b++) {
/* lane m175 pass 2: the tail kernel reads TWO B groups where the mikro reads one, and
its prefetch optimum was never measured -- it was copied from the mikro's 4. Judge
bench rig/kbench.cpp, n=10/arm, spread 0.03 %: 4 -> 0.6746, 2 -> 0.6550, 0 -> 0.6690
cycles/FMA. TWO is this kernel's optimum. */
__builtin_prefetch(bb0 + 1024);
__builtin_prefetch(bb1 + 1024);
K2W(0); K2W(1); K2W(2); K2W(3);
aa += 32; bb0 += 4 * NR * 8; bb1 += 4 * NR * 8; /* A +4 dbl, B +4 lines */
}
for (int i = 0; i < rem; i++) { K2W(0); aa += 8; bb0 += NR * 8; bb1 += NR * 8; }
}
#undef K2W
{
__m256d a0 = c00, a1 = c01, a2 = c02, a3 = c03;
for (int r = 0; r < 2; r++) {
T *p = C + (size_t)r * ldc;
if (accumulate) { ST(p, _mm256_add_pd(LD(p), a0)); ST(p+4, _mm256_add_pd(LD(p+4), a1));
ST(p+8, _mm256_add_pd(LD(p+8), a2)); ST(p+12, _mm256_add_pd(LD(p+12), a3)); }
else if (nt) { STNT(p, a0); STNT(p+4, a1); STNT(p+8, a2); STNT(p+12, a3); }
else { ST(p, a0); ST(p+4, a1); ST(p+8, a2); ST(p+12, a3); }
if (C2) {
T *q2 = C2 + (size_t)r * ldc2;
if (acc2 > 0) { ST(q2, _mm256_add_pd(LD(q2), a0)); ST(q2+4, _mm256_add_pd(LD(q2+4), a1));
ST(q2+8, _mm256_add_pd(LD(q2+8), a2)); ST(q2+12, _mm256_add_pd(LD(q2+12), a3)); }
else if (acc2 < 0) { ST(q2, _mm256_sub_pd(LD(q2), a0)); ST(q2+4, _mm256_sub_pd(LD(q2+4), a1));
ST(q2+8, _mm256_sub_pd(LD(q2+8), a2)); ST(q2+12, _mm256_sub_pd(LD(q2+12), a3)); }
else if (nt2) { STNT(q2, a0); STNT(q2+4, a1); STNT(q2+8, a2); STNT(q2+12, a3); }
else { ST(q2, a0); ST(q2+4, a1); ST(q2+8, a2); ST(q2+12, a3); }
}
a0 = c10; a1 = c11; a2 = c12; a3 = c13;
}
}
}
static void gemm_core(
int m, int n, int k,
const T *a1, const T *a2, int oa, int lda,
const T *b1, const T *b2, int ob, int ldb,
T *C, int ldc, int accumulate,
T *C2, int ldc2, int acc2, int s2) {
if (m <= 0 || n <= 0 || k <= 0) return;
/* E8NT: 32B-aligned write-only destinations may use NT stores */
const int nt = ((((unsigned long)C) & 31) == 0) && ((ldc & 3) == 0);
const int nt2 = C2 && ((((unsigned long)C2) & 31) == 0) && ((ldc2 & 3) == 0);
for (int jc = 0; jc < n; jc += NC) {
int nc = n - jc < NC ? n - jc : NC;
for (int pc = 0; pc < k; pc += KC) {
int kc = k - pc < KC ? k - pc : KC;
packB(kc, nc, b1 + (size_t)pc * ldb + jc,
b2 ? b2 + (size_t)pc * ldb + jc : (const T *)0, ob, ldb, Bpack);
int acc0 = accumulate || pc > 0;
for (int ic = 0; ic < m; ic += MC) {
int mc = m - ic < MC ? m - ic : MC;
packA(mc, kc, a1 + (size_t)ic * lda + pc,
a2 ? a2 + (size_t)ic * lda + pc : (const T *)0, oa, lda, Apack);
T *Cp = C + (size_t)ic * ldc + jc;
int acc = acc0; /* (pc,ic) blocks cover DISJOINT rows: only pc>0 accumulates */
int ac2 = acc2 ? (acc2 > 0 ? 1 : -1) : 0; /* lane m8: sign-carrying accumulate flag */
if (pc > 0 && ac2 == 0) ac2 = 1;
if ((mc % MR) == 2 && (nc % (2 * NR)) == 0) {
/* lane m175: the 2-row row-tail meets a full-width column cover -> WIDE TAIL
(2 rows x 16 columns): 8 accumulators, 2.0 FMA/cycle, vs micro2's 0.98. */
const T *apt = Apack + (size_t)((mc - 2) / MR) * MR * KCP;
for (int jr = 0; jr < nc; jr += 2 * NR) {
const T *bpt = Bpack + (size_t)(jr / NR) * KC * NR;
micro2w(kc, apt, bpt, bpt + (size_t)KC * NR, Cp + (size_t)(mc - 2) * ldc + jr,
ldc, acc, nt,
C2 ? C2 + (size_t)(ic + mc - 2) * ldc2 + jr : (T *)0, ldc2, ac2, nt2);
}
continue;
}
for (int jr = 0; jr < nc; jr += NR) {
int ncc = nc - jr < NR ? nc - jr : NR;
for (int ir = 0; ir + MR <= mc; ir += MR) {
const T *ap = Apack + (size_t)(ir / MR) * MR * KCP;
const T *bp = Bpack + (size_t)(jr / NR) * KC * NR;
micro(kc, ap, bp, Cp + (size_t)ir * ldc + jr, ldc, acc, MR, ncc, nt,
C2 ? C2 + (size_t)(ic + ir) * ldc2 + jr : (T *)0, ldc2, ac2, s2, nt2);
}
int ir = mc - (mc % MR);
if (ir < mc) { /* row tail: phantom rows must not be stored */
const T *ap = Apack + (size_t)(ir / MR) * MR * KCP;
const T *bp = Bpack + (size_t)(jr / NR) * KC * NR;
if (mc - ir == 2 && ncc == NR)
micro2(kc, ap, bp, Cp + (size_t)ir * ldc + jr, ldc, acc, nt,
C2 ? C2 + (size_t)(ic + ir) * ldc2 + jr : (T *)0, ldc2, ac2, nt2);
else
micro(kc, ap, bp, Cp + (size_t)ir * ldc + jr, ldc, acc, mc - ir, ncc, nt,
C2 ? C2 + (size_t)(ic + ir) * ldc2 + jr : (T *)0, ldc2, ac2, s2, nt2);
}
}
}
}
if (nt) __asm__ volatile("sfence" ::: "memory");
}
}
static void blk_axpy(int m, int n, T *dst, int ldd, const T *src, int lds, T s, int accumulate) {
for (int i = 0; i < m; i++) {
T *d = dst + (size_t)i * ldd;
const T *x = src + (size_t)i * lds;
int j = 0;
if (accumulate) {
if (s > 0) {
for (; j + 4 <= n; j += 4) ST(d + j, _mm256_add_pd(LD(d + j), LD(x + j)));
for (; j < n; j++) d[j] += x[j];
} else {
for (; j + 4 <= n; j += 4) ST(d + j, _mm256_sub_pd(LD(d + j), LD(x + j)));
for (; j < n; j++) d[j] -= x[j];
}
} else {
if (s > 0) {
for (; j + 4 <= n; j += 4) ST(d + j, LD(x + j));
for (; j < n; j++) d[j] = x[j];
} else {
for (; j + 4 <= n; j += 4) ST(d + j, _mm256_sub_pd(_mm256_setzero_pd(), LD(x + j)));
for (; j < n; j++) d[j] = -x[j];
}
}
}
}
/* ============ LEAF-FUSED RECURSIVE STRASSEN (spliced driver) ============
* Same 7-product Strassen decomposition, but the LAST level does not
* materialise its operand sums: gemm_core already takes (s2, op) and folds the
* elementwise add/sub into its packing loop. That deletes the whole
* bottom-level sadd pass (2 reads + 1 write of m*m per operand) and the
* re-read of the materialised sum.
*/
#define NIL ((const T *)0)
#if GC_NT
#define GC(a1,a2,oa,lda,b1,b2,ob,ldb,c,ldc,acc) gemm_core(m,m,m,a1,a2,oa,lda,b1,b2,ob,ldb,c,ldc,acc,0)
#define GC0(sz,a1,a2,oa,lda,b1,b2,ob,ldb,c,ldc,acc) gemm_core(sz,sz,sz,a1,a2,oa,lda,b1,b2,ob,ldb,c,ldc,acc,0)
#else
#define GC(a1,a2,oa,lda,b1,b2,ob,ldb,c,ldc,acc) gemm_core(m,m,m,a1,a2,oa,lda,b1,b2,ob,ldb,c,ldc,acc,0,0,0,0)
#define GC0(sz,a1,a2,oa,lda,b1,b2,ob,ldb,c,ldc,acc) gemm_core(sz,sz,sz,a1,a2,oa,lda,b1,b2,ob,ldb,c,ldc,acc,0,0,0,0)
/* lane m8: two-destination product. c2==0 => identical to GC. */
#define GC2(a1,a2,oa,lda,b1,b2,ob,ldb,c,ldc,acc,c2,ldc2,acc2,s2) gemm_core(m,m,m,a1,a2,oa,lda,b1,b2,ob,ldb,c,ldc,acc,c2,ldc2,acc2,s2)
#endif
static void sadd(int m, int n, T *dst, int ldd, const T *s1, const T *s2, int sign, int ld1, int ld2) {
for (int i = 0; i < m; i++) {
T *d = dst + (size_t)i * ldd;
const T *p = s1 + (size_t)i * ld1;
const T *q = s2 ? s2 + (size_t)i * ld2 : (const T *)0;
int j = 0;
if (!q) { for (; j + 4 <= n; j += 4) ST(d + j, LD(p + j)); for (; j < n; j++) d[j] = p[j]; continue; }
if (sign > 0) { for (; j + 4 <= n; j += 4) ST(d + j, _mm256_add_pd(LD(p + j), LD(q + j)));
for (; j < n; j++) d[j] = p[j] + q[j]; }
else { for (; j + 4 <= n; j += 4) ST(d + j, _mm256_sub_pd(LD(p + j), LD(q + j)));
for (; j < n; j++) d[j] = p[j] - q[j]; }
}
}
static void dsgemm(int n, const T *A, const T *B, T *C, int lda, int ldb, int ldc, T *scratch) {
if (n <= SCUT) { GC0(n, A, 0, 0, lda, B, 0, 0, ldb, C, ldc, 0); return; }
int m = n >> 1;
int leaf = (m <= SCUT);
const T *A11 = A, *A12 = A + m, *A21 = A + (size_t)m * lda, *A22 = A + (size_t)m * lda + m;
const T *B11 = B, *B12 = B + m, *B21 = B + (size_t)m * ldb, *B22 = B + (size_t)m * ldb + m;
T *C11 = C, *C12 = C + m, *C21 = C + (size_t)m * ldc, *C22 = C + (size_t)m * ldc + m;
T *M = scratch;
T *TA = leaf ? (T *)0 : scratch + (size_t)m * m;
T *TB = leaf ? (T *)0 : scratch + 2 * (size_t)m * m;
T *sub = scratch + (leaf ? (size_t)m * m : 3 * (size_t)m * m);
/* one product: A-operand = Aa <op> Ab (op 0=copy,1=+,2=-), B likewise, into Dst */
#define PROD(oa, Aa, Ab, ob, Ba, Bb, Dst, lddst, acc) \
do { \
if (leaf) { \
GC(Aa, Ab, oa, lda, Ba, Bb, ob, ldb, Dst, lddst, acc); \
} else { \
const T *_a = Aa, *_b = Ba; int _la = lda, _lb = ldb; \
if (oa) { sadd(m, m, TA, m, Aa, Ab, oa, lda, lda); _a = TA; _la = m; } \
if (ob) { sadd(m, m, TB, m, Ba, Bb, ob, ldb, ldb); _b = TB; _lb = m; } \
dsgemm(m, _a, _b, (Dst), _la, _lb, (lddst), sub); \
} \
} while (0)
/* lane m8 x2: TWO-DESTINATION LEAF. `M` is never materialised: each Strassen product writes
its (one or two) C quadrants directly, with the accumulate flag and sign the identity implies.
C11 = P1 +P4 -P5 +P7 C12 = P3 +P5 C21 = P2 +P4 C22 = P1 -P2 +P3 +P6
P1 MUST be first: it OVERWRITES C11 and C22 (the shipped engine already overwrote C11 with
`blk_axpy2(...,accA=0)`), which is what makes the scheme legal at a CHILD level as well,
where `C` is the parent's recycled `M` plane. P6/P7 then accumulate, exactly as shipped.
acc2: +1 add, -1 subtract, 0 overwrite (and pc>0 turns an overwrite into an accumulate). */
if (leaf) {
GC2(A11, A22, +1, lda, B11, B22, +1, ldb, C11, ldc, 0, C22, ldc, 0, +1);
GC2(A21, A22, +1, lda, B11, NIL, 0, ldb, C21, ldc, 0, C22, ldc, -1, +1);
GC2(A11, NIL, 0, lda, B12, B22, -1, ldb, C22, ldc, 1, C12, ldc, 0, +1);
GC2(A22, NIL, 0, lda, B21, B11, -1, ldb, C11, ldc, 1, C21, ldc, 1, +1);
GC2(A11, A12, +1, lda, B22, NIL, 0, ldb, C12, ldc, 1, C11, ldc, -1, +1);
#if FUSE67
GC2(A21, A11, -1, lda, B11, B12, +1, ldb, C22, ldc, 1, (T *)0, 0, 0, +1);
GC2(A12, A22, -1, lda, B21, B22, +1, ldb, C11, ldc, 1, (T *)0, 0, 0, +1);
#else
PROD(-1, A21, A11, +1, B11, B12, M, m, 0); blk_axpy(m, m, C22, ldc, M, m, +1, 1);
PROD(-1, A12, A22, +1, B21, B22, M, m, 0); blk_axpy(m, m, C11, ldc, M, m, +1, 1);
#endif
return;
}
PROD(+1, A11, A22, +1, B11, B22, M, m, 0); blk_axpy2(m, m, C11, ldc, +1, 0, C22, ldc, +1, 0, M, m);
PROD(+1, A21, A22, 0, B11, NIL, M, m, 0); blk_axpy2(m, m, C21, ldc, +1, 0, C22, ldc, -1, 1, M, m);
PROD( 0, A11, NIL, -1, B12, B22, M, m, 0); blk_axpy2(m, m, C12, ldc, +1, 0, C22, ldc, +1, 1, M, m);
PROD( 0, A22, NIL, -1, B21, B11, M, m, 0); blk_axpy2(m, m, C11, ldc, +1, 1, C21, ldc, +1, 1, M, m);
PROD(+1, A11, A12, 0, B22, NIL, M, m, 0); blk_axpy2(m, m, C11, ldc, -1, 1, C12, ldc, +1, 1, M, m);
#if FUSE67
if (!leaf) {
PROD(-1, A21, A11, +1, B11, B12, M, m, 0); blk_axpy(m, m, C22, ldc, M, m, +1, 1);
PROD(-1, A12, A22, +1, B21, B22, M, m, 0); blk_axpy(m, m, C11, ldc, M, m, +1, 1);
}
#else
PROD(-1, A21, A11, +1, B11, B12, M, m, 0); blk_axpy(m, m, C22, ldc, M, m, +1, 1);
PROD(-1, A12, A22, +1, B21, B22, M, m, 0); blk_axpy(m, m, C11, ldc, M, m, +1, 1);
#endif
#undef PROD
}
static T *g_pool = 0;
static size_t g_pool_n = 0;
static void strassen1(int n, const T *A, const T *B, T *C, int ldc) {
size_t need = 2 * (size_t)n * n + 1024; /* lane m8 x2: headroom for the non-leaf planes */
if (g_pool_n < need) {
free(g_pool);
if (posix_memalign((void **)&g_pool, 64, need * sizeof(T))) g_pool = 0;
g_pool_n = g_pool ? need : 0;
}
if (!g_pool) { GC0(n, A, 0, 0, n, B, 0, 0, n, C, n, 0); return; }
dsgemm(n, A, B, C, n, n, ldc, g_pool);
}
#pragma GCC pop_options
#if defined(LOCAL_DRIVER) || defined(LOCAL_ONLY)
#else
void matrix_multiply(int n, const double *A, const double *B, double *C) {
if (n > SCUT) strassen1(n, A, B, C, n);
else gemm_core(n, n, n, A, 0, 0, n, B, 0, 0, n, C, n, 0, 0, 0, 0, 0);
}
#endif
/*
*/
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 38.113 ms | 10 MB + 52 KB | Accepted | Score: 100 | 显示更多 |