#define FUSE67 1
#define GC_NT 0
#define MC 6
#define NC 1024
#define KC 512
#define KCP (KC+8)
#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;
static inline unsigned long long MIKRD(void) {
unsigned lo, hi; __asm__ volatile("rdtsc" : "=a"(lo), "=d"(hi));
return ((unsigned long long)hi << 32) | lo;
}
static unsigned long long g_mik, g_kloop, g_last, g_ph;
static volatile unsigned long long g_sink;
#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 Apack[(size_t)((MC + MR - 1) / MR) * MR * KCP + 64] __attribute__((aligned(64)));
static T Bpack[(size_t)((NC + NR - 1) / NR) * NR * KC + 64] __attribute__((aligned(64)));
// 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;
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). */
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;
}
}
}
}
/* ================= lane m206 -- RELOCATED DEAD SCRATCH (CSCR) =================
MEASURED, on the board, this session: 32 768 FRESH 4 KiB pages added to the footprint cost
+8.089 ms (sid 118007: 287.727 vs 279.638) = 0.2469 us/page = 889 ticks/page. Our pool is
~37 700 touched pages = ~9.3 ms of first-touch faults; the ten sum frames are 20 157 of those.
The caller's C buffer is n*n doubles (32 MiB) that the HARNESS zeroes (problem page: an
output-only array starts at 0), so its pages are ALREADY MAPPED AND WARM. fused_comb writes
every element of C at the end and leaves nothing stale, so frames that are DEAD before the
combine may live in C's pages: same traffic, same layout discipline (+128-double skew kept),
and ~26.3 MiB of first-touch faults that are then never paid.
Relocated: the LEAF level's M1..M5 (5 x 266 368 doubles) + two top-level sum frames.
No arithmetic changes at all -- the same frames at different addresses. */
static T *g_rl = 0; /* relocation base = C; 0 disables the whole trick */
#define RL_LEAF 32 /* leaf scratch at C+32 doubles (256 B: 64-B aligned, not 4K) */
/* lane m207 -- THE LEAF NOW SEATS 7 FRAMES, NOT 5. Under the fix below the leaf's two accumulating
products P6/P7 materialise into M6/M7 (the leaf's frames 6 and 7) instead of accumulating into the
destination, so the leaf's carve-out is 7 x 266 368 = 1 864 576 doubles and the two relocated sums
start after it. n=2048: RL_NEED = 3 978 464 <= n*n = 4 194 304. */
/* lane m209 ILV: the leaf's seven M planes are now INTERLEAVED row-wise, so the leaf's
frame block is 512 * (7*520) + 128 = 1 863 808 doubles, not 7 * 266 368. */
#define RL_SUM 1863840 /* = RL_LEAF + (512*3640 + 128) */
#define RL_NEED (RL_SUM + 2 * 1056896 + 64)
#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) {
__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();
/* 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. */
/* lane m78 csp_acc: warm C ONLY on the accumulating (RMW) path. With `nt` set the pure
stores are vmovntpd -- they bypass L1, so a prefetcht0 of the destination fetches a line
that is then overwritten without ever being read: 6 x 64 B of wasted DRAM per call. */
if (accumulate) {
if (mrows > 0) __builtin_prefetch(C + 0 * ldc, 1, 3);
if (mrows > 1) __builtin_prefetch(C + 1 * ldc, 1, 3);
if (mrows > 2) __builtin_prefetch(C + 2 * ldc, 1, 3);
if (mrows > 3) __builtin_prefetch(C + 3 * ldc, 1, 3);
if (mrows > 4) __builtin_prefetch(C + 4 * ldc, 1, 3);
if (mrows > 5) __builtin_prefetch(C + 5 * ldc, 1, 3);
}
int k = 0;
for (; k + 5 <= kc; k += 5) {
__builtin_prefetch(bp + (size_t)(k + 16 + 0) * NR);
__builtin_prefetch(bp + (size_t)(k + 16 + 1) * NR);
__builtin_prefetch(bp + (size_t)(k + 16 + 2) * NR);
__builtin_prefetch(bp + (size_t)(k + 16 + 3) * NR);
__builtin_prefetch(bp + (size_t)(k + 16 + 4) * NR);
KSTEP(0); KSTEP(1); KSTEP(2); KSTEP(3); KSTEP(4); }
for (; k + 4 <= kc; k += 4) {
__builtin_prefetch(bp + (size_t)(k + 16 + 0) * NR);
__builtin_prefetch(bp + (size_t)(k + 16 + 1) * NR);
__builtin_prefetch(bp + (size_t)(k + 16 + 2) * NR);
__builtin_prefetch(bp + (size_t)(k + 16 + 3) * NR);
KSTEP(0); KSTEP(1); KSTEP(2); KSTEP(3); }
for (; k < kc; k++) KSTEP(0);
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
} 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
}
}
/* lane m118 -- TAIL KERNEL. The mikro's k-loop is hardcoded to MR=6 rows and the `mrows`
argument only guards the STORE block (this is mmml4k's §N237 shape, re-derived here from
src: gemm_core's ic loop is `ic += MC(6)` and the leaf is m=512 = 85*6 + 2, so exactly ONE
group per leaf pays 6 rows of k-loop for 2 -- 4/6 of one group, 6.42 M ticks/row).
microR<2> is the same instruction stream restricted to the two rows that exist. */
template <int R>
__attribute__((always_inline)) static inline void microR(int kc, const T *__restrict ap, const T *__restrict bp,
T *__restrict C, int ldc, int accumulate, int ncols, int nt) {
__m256d c[6][2];
for (int r = 0; r < 6; r++) { c[r][0] = _mm256_setzero_pd(); c[r][1] = _mm256_setzero_pd(); }
int k = 0;
for (; k + 4 <= kc; k += 4) {
__builtin_prefetch(bp + (size_t)(k + 16 + 0) * NR);
__builtin_prefetch(bp + (size_t)(k + 16 + 1) * NR);
__builtin_prefetch(bp + (size_t)(k + 16 + 2) * NR);
__builtin_prefetch(bp + (size_t)(k + 16 + 3) * NR);
for (int u = 0; u < 4; u++) {
__m256d b0 = LD(bp + (size_t)(k + u) * NR), b1 = LD(bp + (size_t)(k + u) * NR + 4);
for (int r = 0; r < R; r++) {
__m256d t = _mm256_broadcast_sd(ap + (size_t)r * KCP + (k + u));
c[r][0] = _mm256_fmadd_pd(t, b0, c[r][0]);
c[r][1] = _mm256_fmadd_pd(t, b1, c[r][1]);
}
}
}
for (; k < kc; k++) {
__m256d b0 = LD(bp + (size_t)k * NR), b1 = LD(bp + (size_t)k * NR + 4);
for (int r = 0; r < R; r++) {
__m256d t = _mm256_broadcast_sd(ap + (size_t)r * KCP + k);
c[r][0] = _mm256_fmadd_pd(t, b0, c[r][0]);
c[r][1] = _mm256_fmadd_pd(t, b1, c[r][1]);
}
}
if (ncols == NR) {
for (int r = 0; r < R; r++) {
T *p = C + (size_t)r * ldc;
if (accumulate) { ST(p, _mm256_add_pd(LD(p), c[r][0])); ST(p + 4, _mm256_add_pd(LD(p + 4), c[r][1])); }
else { if (nt) { STNT(p, c[r][0]); STNT(p + 4, c[r][1]); } else { ST(p, c[r][0]); ST(p + 4, c[r][1]); } }
}
} 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);
for (int r = 0; r < R; r++) {
T *p = C + (size_t)r * ldc;
if (accumulate) {
_mm256_maskstore_pd(p, ml, _mm256_add_pd(_mm256_maskload_pd(p, ml), c[r][0]));
_mm256_maskstore_pd(p + 4, mh, _mm256_add_pd(_mm256_maskload_pd(p + 4, mh), c[r][1]));
} else { _mm256_maskstore_pd(p, ml, c[r][0]); _mm256_maskstore_pd(p + 4, mh, c[r][1]); }
}
}
}
/* 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");
}
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) {
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);
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;
{ unsigned long long _t=MIKRD(); packB(kc, nc, b1 + (size_t)pc * ldb + jc,
b2 ? b2 + (size_t)pc * ldb + jc : (const T *)0, ob, ldb, Bpack); g_ph += MIKRD()-_t; }
int acc0 = accumulate || pc > 0;
for (int ic = 0; ic < m; ic += MC) {
int mc = m - ic < MC ? m - ic : MC;
{ unsigned long long _t=MIKRD(); packA(mc, kc, a1 + (size_t)ic * lda + pc,
a2 ? a2 + (size_t)ic * lda + pc : (const T *)0, oa, lda, Apack); g_ph += MIKRD()-_t; }
T *Cp = C + (size_t)ic * ldc + jc;
int acc = acc0; /* (pc,ic) blocks cover DISJOINT rows: only pc>0 accumulates */
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 * KC;
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);
}
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 * KC;
const T *bp = Bpack + (size_t)(jr / NR) * KC * NR;
T *cp = Cp + (size_t)ir * ldc + jr;
if (mc - ir == 2) microR<2>(kc, ap, bp, cp, ldc, acc, ncc, nt);
else micro(kc, ap, bp, cp, ldc, acc, mc - ir, ncc, nt);
}
}
}
}
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)
#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)
#endif
/* ================= lane m125 -- SADD5 ================= */
static void sadd5(int m, T *S1, T *S2, T *S5, T *S6, T *S7, int lds,
const T *q11, const T *q12, const T *q21, const T *q22, int ldq) {
for (int i = 0; i < m; i++) {
const T *a11 = q11 + (size_t)i * ldq, *a12 = q12 + (size_t)i * ldq,
*a21 = q21 + (size_t)i * ldq, *a22 = q22 + (size_t)i * ldq;
T *d1 = S1 + (size_t)i * lds, *d2 = S2 + (size_t)i * lds, *d5 = S5 + (size_t)i * lds,
*d6 = S6 + (size_t)i * lds, *d7 = S7 + (size_t)i * lds;
if (i + 1 < m) { /* every stream's next row starts ldq*8 bytes away */
__builtin_prefetch(q11 + (size_t)(i + 1) * ldq, 0, 1);
__builtin_prefetch(q12 + (size_t)(i + 1) * ldq, 0, 1);
__builtin_prefetch(q21 + (size_t)(i + 1) * ldq, 0, 1);
__builtin_prefetch(q22 + (size_t)(i + 1) * ldq, 0, 1);
__builtin_prefetch(S1 + (size_t)(i + 1) * lds, 1, 3);
__builtin_prefetch(S2 + (size_t)(i + 1) * lds, 1, 3);
__builtin_prefetch(S5 + (size_t)(i + 1) * lds, 1, 3);
__builtin_prefetch(S6 + (size_t)(i + 1) * lds, 1, 3);
__builtin_prefetch(S7 + (size_t)(i + 1) * lds, 1, 3);
}
int j = 0;
/* FIVE NT destination streams interleaved at 32 B per iteration thrash the WC
buffers: the phase measured 4.4 B/tick instead of the 6.0 a 3-stream pass gets
(lane m125, judge, this session). m118's fix for the same signature -- put each
stream's two 32-B halves back to back -- i.e. unroll j by 8. */
for (; j + 8 <= m; j += 8) {
__m256d v11 = LD(a11 + j), v11b = LD(a11 + j + 4);
__m256d v12 = LD(a12 + j), v12b = LD(a12 + j + 4);
__m256d v21 = LD(a21 + j), v21b = LD(a21 + j + 4);
__m256d v22 = LD(a22 + j), v22b = LD(a22 + j + 4);
STNT(d1 + j, _mm256_add_pd(v11, v22));
STNT(d1 + j + 4, _mm256_add_pd(v11b, v22b));
STNT(d2 + j, _mm256_add_pd(v21, v22));
STNT(d2 + j + 4, _mm256_add_pd(v21b, v22b));
STNT(d5 + j, _mm256_add_pd(v11, v12));
STNT(d5 + j + 4, _mm256_add_pd(v11b, v12b));
STNT(d6 + j, _mm256_sub_pd(v21, v11));
STNT(d6 + j + 4, _mm256_sub_pd(v21b, v11b));
STNT(d7 + j, _mm256_sub_pd(v12, v22));
STNT(d7 + j + 4, _mm256_sub_pd(v12b, v22b));
}
for (; j + 4 <= m; j += 4) {
__m256d v11 = LD(a11 + j), v12 = LD(a12 + j), v21 = LD(a21 + j), v22 = LD(a22 + j);
STNT(d1 + j, _mm256_add_pd(v11, v22));
STNT(d2 + j, _mm256_add_pd(v21, v22));
STNT(d5 + j, _mm256_add_pd(v11, v12));
STNT(d6 + j, _mm256_sub_pd(v21, v11));
STNT(d7 + j, _mm256_sub_pd(v12, v22));
}
for (; j < m; j++) {
T x11 = a11[j], x12 = a12[j], x21 = a21[j], x22 = a22[j];
d1[j] = x11 + x22; d2[j] = x21 + x22; d5[j] = x11 + x12;
d6[j] = x21 - x11; d7[j] = x12 - x22;
}
}
__asm__ volatile("sfence" ::: "memory");
}
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;
if (i + 1 < m) { __builtin_prefetch(s1 + (size_t)(i + 1) * ld1, 0, 3);
if (s2) __builtin_prefetch(s2 + (size_t)(i + 1) * ld2, 0, 3);
__builtin_prefetch(dst + (size_t)(i + 1) * ldd, 1, 3); }
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) STNT(d + j, LD(p + j)); for (; j < n; j++) d[j] = p[j]; continue; }
if (sign > 0) { for (; j + 4 <= n; j += 4) STNT(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) STNT(d + j, _mm256_sub_pd(LD(p + j), LD(q + j)));
for (; j < n; j++) d[j] = p[j] - q[j]; }
} __asm__ volatile("sfence" ::: "memory");
}
/* lane m76 -- the SINGLE-PASS FUSED COMBINE. One pass over the m x m plane that reads the five
product planes M1..M5 and writes each of the four C quadrants exactly ONCE, instead of reading a
single reused M five times and touching the quadrants sixteen times. The expressions reproduce
blk_axpy2's accumulation order exactly, so this must be BIT-IDENTICAL to the base. */
#if 1
#define FCOMB_ST(p, v) STNT(p, v)
#else
#define FCOMB_ST(p, v) ST(p, v)
#endif
static void fused_comb(int n, T *C11, T *C12, T *C21, T *C22, int ldc,
const T *M1, const T *M2, const T *M3, const T *M4, const T *M5,
const T *M6, const T *M7, int Mp) {
for (int i = 0; i < n; i++) {
const T *m1 = M1 + (size_t)i * Mp, *m2 = M2 + (size_t)i * Mp, *m3 = M3 + (size_t)i * Mp,
*m4 = M4 + (size_t)i * Mp, *m5 = M5 + (size_t)i * Mp;
T *p11 = C11 + (size_t)i * ldc, *p12 = C12 + (size_t)i * ldc,
*p21 = C21 + (size_t)i * ldc, *p22 = C22 + (size_t)i * ldc;
int j = 0;
if (M6) { /* NON-LEAF: P6/P7 are materialised, so fold them in -- ONE add each, same order */
for (; j + 8 <= n; j += 8) {
__m256d a1 = LD(m1 + j), a2 = LD(m2 + j), a3 = LD(m3 + j), a4 = LD(m4 + j), a5 = LD(m5 + j);
__m256d b1 = LD(m1 + j + 4), b2 = LD(m2 + j + 4), b3 = LD(m3 + j + 4), b4 = LD(m4 + j + 4), b5 = LD(m5 + j + 4);
__m256d a7 = LD(M7 + (size_t)i * Mp + j), b7 = LD(M7 + (size_t)i * Mp + j + 4);
__m256d a6 = LD(M6 + (size_t)i * Mp + j), b6 = LD(M6 + (size_t)i * Mp + j + 4);
FCOMB_ST(p11 + j, _mm256_add_pd(_mm256_sub_pd(_mm256_add_pd(a1, a4), a5), a7));
FCOMB_ST(p11 + j + 4, _mm256_add_pd(_mm256_sub_pd(_mm256_add_pd(b1, b4), b5), b7));
FCOMB_ST(p12 + j, _mm256_add_pd(a3, a5));
FCOMB_ST(p12 + j + 4, _mm256_add_pd(b3, b5));
FCOMB_ST(p21 + j, _mm256_add_pd(a2, a4));
FCOMB_ST(p21 + j + 4, _mm256_add_pd(b2, b4));
FCOMB_ST(p22 + j, _mm256_add_pd(_mm256_add_pd(_mm256_sub_pd(a1, a2), a3), a6));
FCOMB_ST(p22 + j + 4, _mm256_add_pd(_mm256_add_pd(_mm256_sub_pd(b1, b2), b3), b6));
}
for (; j + 4 <= n; j += 4) {
__m256d a1 = LD(m1 + j), a2 = LD(m2 + j), a3 = LD(m3 + j), a4 = LD(m4 + j), a5 = LD(m5 + j);
FCOMB_ST(p11 + j, _mm256_add_pd(_mm256_sub_pd(_mm256_add_pd(a1, a4), a5), LD(M7 + (size_t)i * Mp + j)));
FCOMB_ST(p22 + j, _mm256_add_pd(_mm256_add_pd(_mm256_sub_pd(a1, a2), a3), LD(M6 + (size_t)i * Mp + j)));
FCOMB_ST(p12 + j, _mm256_add_pd(a3, a5));
FCOMB_ST(p21 + j, _mm256_add_pd(a2, a4));
}
for (; j < n; j++) { const T *z6 = M6 + (size_t)i * Mp, *z7 = M7 + (size_t)i * Mp;
p11[j] = ((m1[j] + m4[j]) - m5[j]) + z7[j];
p22[j] = ((m1[j] - m2[j]) + m3[j]) + z6[j];
p12[j] = m3[j] + m5[j]; p21[j] = m2[j] + m4[j]; }
continue;
}
for (; j + 8 <= n; j += 8) {
__m256d a1 = LD(m1 + j), a2 = LD(m2 + j), a3 = LD(m3 + j), a4 = LD(m4 + j), a5 = LD(m5 + j);
__m256d b1 = LD(m1 + j + 4), b2 = LD(m2 + j + 4), b3 = LD(m3 + j + 4), b4 = LD(m4 + j + 4), b5 = LD(m5 + j + 4);
FCOMB_ST(p11 + j, _mm256_sub_pd(_mm256_add_pd(a1, a4), a5));
FCOMB_ST(p11 + j + 4, _mm256_sub_pd(_mm256_add_pd(b1, b4), b5));
FCOMB_ST(p12 + j, _mm256_add_pd(a3, a5));
FCOMB_ST(p12 + j + 4, _mm256_add_pd(b3, b5));
FCOMB_ST(p21 + j, _mm256_add_pd(a2, a4));
FCOMB_ST(p21 + j + 4, _mm256_add_pd(b2, b4));
FCOMB_ST(p22 + j, _mm256_add_pd(_mm256_sub_pd(a1, a2), a3));
FCOMB_ST(p22 + j + 4, _mm256_add_pd(_mm256_sub_pd(b1, b2), b3));
}
for (; j + 4 <= n; j += 4) {
__m256d a1 = LD(m1 + j), a2 = LD(m2 + j), a3 = LD(m3 + j), a4 = LD(m4 + j), a5 = LD(m5 + j);
FCOMB_ST(p11 + j, _mm256_sub_pd(_mm256_add_pd(a1, a4), a5));
FCOMB_ST(p22 + j, _mm256_add_pd(_mm256_sub_pd(a1, a2), a3));
FCOMB_ST(p12 + j, _mm256_add_pd(a3, a5));
FCOMB_ST(p21 + j, _mm256_add_pd(a2, a4));
}
for (; j < n; j++) {
p11[j] = (m1[j] + m4[j]) - m5[j];
p22[j] = (m1[j] - m2[j]) + m3[j];
p12[j] = m3[j] + m5[j];
p21[j] = m2[j] + m4[j];
}
}
#if 1
__asm__ volatile("sfence" ::: "memory");
#endif
}
#undef FCOMB_ST
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;
/* m43: PAD THE M STRIDE. M's rows are ldc = m doubles apart, i.e. 8192 B at m=1024 and
4096 B at m=512 -- BOTH exact multiples of 4096, so the 6 rows the mikro stores per call
all land in the SAME L1 set and the same 4K page offset. +8 doubles per row rotates the
set by 1. Bit-identical by construction: the pad columns are never written (nc == m) and
never read (all loops run j < n). */
const int Mp = m + 8;
/* lane m209 -- ILV (interleaved M planes) AT THE LEAF.
fused_comb reads the seven M planes IN PARALLEL at the same in-frame offset. With the
shipped layout they are MPL = m*Mp+128 doubles apart: 2 130 944 B at the leaf, i.e. only
TWO distinct 4 KiB page offsets and FOUR distinct L1 set-groups for SEVEN concurrent read
streams. ILV puts plane k at base + k*Mp with row stride 7*Mp, so the seven planes of a
row are ONE contiguous 3640-double window and the seven streams land on SEVEN distinct page
offsets and SEVEN distinct L1 sets. Addresses only: same values, same order, BIT-IDENTICAL.
(mmmd4k boarded the same construction for 5 planes: `ILV`, sid 118126/118133, -9 ms.) */
/* m76: the five planes are m*Mp doubles apart = 8519680 B, an EXACT multiple of 4096, so all
five would be 4K-congruent -- the same aliasing class the boarded pad fixed. Each plane
therefore carries a +64-double (512 B) skew. */
const size_t MPL = (size_t)m * Mp + 128;
const int MpR = leaf ? 7 * Mp : Mp;
const size_t MST = leaf ? (size_t)Mp : MPL;
T *M1 = scratch, *M2 = M1 + MST, *M3 = M2 + MST, *M4 = M3 + MST, *M5 = M4 + MST,
*M6 = M5 + MST, *M7 = M6 + MST;
/* m125 SADD5: five A-sums and five B-sums, each in its own frame. */
T *S1 = M7 + MST, *S2 = S1 + MPL, *S5 = S2 + MPL, *S6 = S5 + MPL, *S7 = S6 + MPL;
T *T1 = S7 + MPL, *T3 = T1 + MPL, *T4 = T3 + MPL, *T6 = T4 + MPL, *T7 = T6 + MPL;
T *sub = scratch + 17 * MPL;
if (g_rl) { /* m206 CSCR: dead-until-the-combine frames live in C's pre-touched pages */
if (!leaf) { T *r = g_rl + RL_SUM; S6 = r; S7 = r + MPL; }
sub = g_rl + RL_LEAF; /* the leaf's M1..M5 are carved from here */
}
/* ================= lane m206 -- M-PLANE OVERLAY =================
The seven M planes and the ten operand sums are 17 distinct frames, and the PAGE BILL is
charged on every one of them (measured: 0.082 us/page by elimination, sid 118019). But each
M plane is written by its OWN product and stays alive only until the combine, while each sum
is dead the moment its LAST consumer has run. So an M plane may REUSE the frame of a sum that
is already dead -- the page is then touched ONCE instead of twice, with ZERO extra traffic,
no reordering (the shipped PROD order already satisfies every assignment) and bit-identical
output. Last consumers in the shipped order: S1->P1 S2->P2 S5->P5 S6->P6 S7->P7 ; T1->P1,P2
T3->P3 T4->P4 T6->P6 T7->P7. M1 is first, so it takes a fresh frame; M2..M7 overlay. */
if (!leaf) { M1 = T1; M2 = S1; M3 = S2; M4 = T3; M5 = T4; M6 = S5; M7 = T6; }
if (!leaf) {
{ unsigned long long _t=MIKRD(); sadd5(m, S1, S2, S5, S6, S7, MpR, A11, A12, A21, A22, lda); g_ph += MIKRD()-_t; }
/* The B-side sums are a PERMUTATION of the same five forms: T1=B11+B22 (F1),
T3=B12-B22 (F5), T4=B21-B11 (F4), T6=B11+B12 (F3), T7=B21+B22 (F2). sadd5 emits
F1..F5 in that order, so the B destinations are passed permuted. */
{ unsigned long long _t=MIKRD(); sadd5(m, T1, T7, T6, T4, T3, MpR, B11, B12, B21, B22, ldb); g_ph += MIKRD()-_t; }
}
/* one product: A-operand = Aa <op> Ab (op 0=copy,1=+,2=-), B likewise, into Dst */
#define PROD(oa, Aa, Ab, As, ob, Ba, Bb, Bs, Dst, lddst, acc) \
do { \
if (leaf) { \
GC(Aa, Ab, oa, lda, Ba, Bb, ob, ldb, Dst, lddst, acc); \
} else { \
const T *_a = (oa) ? (As) : (Aa), *_b = (ob) ? (Bs) : (Ba); \
int _la = (oa) ? MpR : lda, _lb = (ob) ? MpR : ldb; \
dsgemm(m, _a, _b, (Dst), _la, _lb, (lddst), sub); \
} \
} while (0)
PROD(+1, A11, A22, S1, +1, B11, B22, T1, M1, MpR, 0);
PROD(+1, A21, A22, S2, 0, B11, NIL, T1, M2, MpR, 0);
PROD( 0, A11, NIL, S1, -1, B12, B22, T3, M3, MpR, 0);
PROD( 0, A22, NIL, S1, -1, B21, B11, T4, M4, MpR, 0);
PROD(+1, A11, A12, S5, 0, B22, NIL, T1, M5, MpR, 0);
#if FUSE67
{ /* lane m207 FIX -- MATERIALISE P6/P7 INTO M6/M7 AT THE LEAF TOO (was: accumulate into C
after the combine). WHY: the leaf's own C11/C22 ARE the destination, and the destination
may now BE one of the leaf's own operands (the M-plane overlay's self-aliased frame), so
the base's order makes P6/P7 read a quadrant the combine has just overwritten. Measured
signature of that failure (rig/mkpat.py, n=2048): exactly the 4 diagonal 512-blocks,
262 144 entries each, worst 9.9e4. BIT-IDENTICAL: the base computes
((M1+M4)-M5) then += P7, which is ((M1+M4)-M5) + M7 in the same add order. */
PROD(-1, A21, A11, S6, +1, B11, B12, T6, M6, MpR, 0);
PROD(-1, A12, A22, S7, +1, B21, B22, T7, M7, MpR, 0);
}
#else
{ PROD(-1, A21, A11, +1, B11, B12, M6, MpR, 0);
PROD(-1, A12, A22, +1, B21, B22, M7, MpR, 0); }
#endif
{ unsigned long long _t=MIKRD(); fused_comb(m, C11, C12, C21, C22, ldc, M1, M2, M3, M4, M5, M6, M7, MpR); g_ph += MIKRD()-_t; }
/* LEAF and NON-LEAF are now the same shape: P6/P7 materialised, the combine folds them. */
#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 _m0 = (size_t)(n >> 1);
/* m125 SADD5: 17 frames at the top level (7 M + 10 sums) and the leaf's own
17-frame carve-out. The leaf never runs sadd5 (its PROD takes the GC path),
so it actually seats 7; the carve-out is sized alike for simplicity. */
size_t _P1 = _m0 * (_m0 + 8) + 128, _P2 = (_m0 / 2) * ((_m0 / 2) + 8) + 128;
size_t need = 17 * _P1 + 17 * _P2 + 64 * (size_t)n + 1024;
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; }
/* m206 CSCR: engage only when the caller's C is big enough to host the relocated frames */
g_rl = ((size_t)n * (size_t)n >= (size_t)RL_NEED) ? C : (T *)0;
dsgemm(n, A, B, C, n, n, ldc, g_pool);
}
#pragma GCC pop_options
void matrix_multiply(int n, const double *A, const double *B, double *C) {
unsigned long long t0 = MIKRD();
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);
unsigned long long t1 = MIKRD(), x = 0;
unsigned long long q = ((t1 - t0) - g_ph); g_ph = 0; g_last = q;
{ unsigned long long target = t1 + (q >> 2);
do { __asm__ volatile("addq $1, %0" : "+r"(x)); } while (MIKRD() < target); }
g_sink = x;
(void)t0; (void)t1;
}
unsigned long long mw_aux(void) { return g_last; }
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 334.924 ms | 98 MB + 568 KB | Accepted | Score: 100 | 显示更多 |