/* m211 arm: A-panel ROW-STRIDE pad, KCPAD=8 (stride KC+8 = 2112 B = 33 lines) */
#define KCPAD 8
#ifndef FUSE67
#define FUSE67 1
#endif
#ifndef FUSE_NL
#define FUSE_NL 1
#endif
#if ILV
#endif
#ifndef ILV
#define ILV 1
#endif
#ifndef GC_NT
#define GC_NT 1
#endif
#ifndef MC_G
#define MC_G 6
#endif
#ifndef NC_G
#define NC_G 1024
#endif
#define KC_G 512
#ifndef SCUT_G
#ifndef SCUT_OVR
#define SCUT_G 512
#else
#define SCUT_G SCUT_OVR
#endif
#endif
// 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>
#define M5_MODE 3
typedef double T;
#ifndef MR
#define MR 6
#endif
#ifndef NR
#define NR 8
#endif
#ifndef KC
#define KC KC_G
#endif
#ifndef MC
#define MC MC_G
#endif
#ifndef NC
#define NC NC_G
#endif
#ifndef SCUT
#define SCUT SCUT_G
#endif
static unsigned long long g_nlm_nodes = 0;
static int g_top = 0;
#pragma GCC optimize("O3","unroll-loops")
#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));
}
// 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)
#ifndef KCPAD
#define KCPAD 0
#endif
#define KCP (KC + KCPAD) /* m211: A-panel ROW STRIDE. At KCPAD=0 the stride is KC = 2048 B,
exactly HALF the L1 set period (64 sets x 64 B = 4096 B), so all
six micro-rows land on two set-groups and the 8-way L1 is ~88 %
occupied before the C stores. KCPAD=8 -> 2112 B = 33 lines. */
/* =================== lane m4k_skel arm knobs =================== */
#define g_nt 1
#define g_pf 0
#define g_pfa 0
#define g_pfb 0
#define g_arm 0
static inline void stn(int nt, T *p, __m256d v){ if(nt) _mm256_stream_pd(p,v); else ST(p,v); }
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 nt = g_nt && !accA && !accB;
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;
if(g_pf){__builtin_prefetch(x+g_pf,0,1);__builtin_prefetch(a+g_pf,0,1);__builtin_prefetch(b+g_pf,0,1);}
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(sA>0) stn(nt,a+j,v); else stn(nt,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(sB>0) stn(nt,b+j,v); else stn(nt,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); }
}
}
static void blk_axpy(int m,int n,T*dst,int ldd,const T*src,int lds,T s,int accumulate){
const int nt = g_nt && !accumulate;
for(int i=0;i<m;i++){
T*d=dst+(size_t)i*ldd; const T*x=src+(size_t)i*lds;
if(g_pf){__builtin_prefetch(x+g_pf,0,1);}
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)stn(nt,d+j,LD(x+j));for(;j<n;j++)d[j]=x[j];}
else {for(;j+4<=n;j+=4)stn(nt,d+j,_mm256_sub_pd(_mm256_setzero_pd(),LD(x+j)));for(;j<n;j++)d[j]=-x[j];}
}
}
}
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;
if(g_pf){__builtin_prefetch(p+g_pf,0,1); if(q)__builtin_prefetch(q+g_pf,0,1);}
int j=0;
if(!q){for(;j+4<=n;j+=4)stn(g_nt,d+j,LD(p+j));for(;j<n;j++)d[j]=p[j];continue;}
if(sign>0){for(;j+4<=n;j+=4)stn(g_nt,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)stn(g_nt,d+j,_mm256_sub_pd(LD(p+j),LD(q+j)));for(;j<n;j++)d[j]=p[j]-q[j];}
}
}
static void comb5(int m,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,int ldm){
const int nt=g_nt;
for(int i=0;i<m;i++){
T*a=C11+(size_t)i*ldc,*b=C12+(size_t)i*ldc,*c=C21+(size_t)i*ldc,*d=C22+(size_t)i*ldc;
const T*p1=m1+(size_t)i*ldm,*p2=m2+(size_t)i*ldm,*p3=m3+(size_t)i*ldm,*p4=m4+(size_t)i*ldm,*p5=m5+(size_t)i*ldm;
if(g_pf){__builtin_prefetch(p1+g_pf,0,1);__builtin_prefetch(p2+g_pf,0,1);__builtin_prefetch(p3+g_pf,0,1);__builtin_prefetch(p4+g_pf,0,1);__builtin_prefetch(p5+g_pf,0,1);}
for(int j=0;j+4<=m;j+=4){
__m256d x1=LD(p1+j),x2=LD(p2+j),x3=LD(p3+j),x4=LD(p4+j),x5=LD(p5+j);
stn(nt,a+j,_mm256_add_pd(_mm256_sub_pd(x1,x5),x4));
stn(nt,b+j,_mm256_add_pd(x3,x5));
stn(nt,c+j,_mm256_add_pd(x2,x4));
stn(nt,d+j,_mm256_sub_pd(_mm256_add_pd(x1,x3),x2));
}
for(int j=(m&~3);j<m;j++){
double x1=p1[j],x2=p2[j],x3=p3[j],x4=p4[j],x5=p5[j];
a[j]=x1-x5+x4; b[j]=x3+x5; c[j]=x2+x4; d[j]=x1-x2+x3;
}
}
}
/* =============================================================== */
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;
if (r + 1 < mr) __builtin_prefetch(s1 + lda, 0, 1); /* E8PF */
else if (p + MR + r < mc) __builtin_prefetch(s1 + (size_t)MR * lda, 0, 1);
for (int _q = 1; _q <= g_pfa; _q++) {
int _r = r + _q;
if (_r < mr) __builtin_prefetch(s1 + (size_t)_q * lda, 0, 1);
else if (p + MR + _r < mc) __builtin_prefetch(s1 + (size_t)(MR + _q) * 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 (ported from mmmd2k, where it is board-priced at -2.21 %): the shipped loop is
g OUTER / k INNER, so it reads 8 doubles (ONE FULL 64-B LINE) per lda*8 = 32 KB jump and
kc of them in a row. THE HW PREFETCHER CANNOT FOLLOW A 32 KB STRIDE, so every line was a
DEMAND MISS. k OUTER / g INNER makes the source read nc*8 bytes CONTIGUOUS per row. */
for (int k0 = 0; k0 < kc; k0 += 8)
for (int g = 0; g + NR <= nc; g += NR) {
if (g_pfb) __builtin_prefetch(b1 + (size_t)(k0 + 8*g_pfb) * lda + g, 0, 1);
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;
}
}
}
#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();
/* lane m164c SPF: the C destination is written in NR-double row segments ldc*8 bytes apart
(16 KB at ldc=2048) and an ACCUMULATING store is a read-modify-write, so its first touch is
a cold line fetch paid at the very end of the call, where the OoO window is closing.
One prefetch per row, issued before the k-loop, is hidden by ~kc*6.5 cycles of mikro work.
Output-preserving by construction (a prefetch cannot change a value). */
if (accumulate) {
for (int r = 0; r < mrows; r++)
_mm_prefetch((const char *)(C + (size_t)r * ldc), _MM_HINT_T0);
}
int k = 0;
for (; k + 4 <= kc; k += 4) {
__builtin_prefetch(bp + (size_t)(k + 16) * NR);
__builtin_prefetch(bp + (size_t)(k + 17) * NR);
__builtin_prefetch(bp + (size_t)(k + 18) * NR);
__builtin_prefetch(bp + (size_t)(k + 19) * NR);
KSTEP(0); KSTEP(1); KSTEP(2); KSTEP(3);
}
for (; k < kc; k++) KSTEP(0);
if (ncols == NR) {
#define STO(r, a, b) \
if (r < mrows) { \
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) { \
_mm256_stream_pd(p, a); _mm256_stream_pd(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
}
}
/* E8FUSE (ported from mmmd2k, board-priced at -0.6 % there): the two combine axpys of one
Strassen product each read the WHOLE m x m scratch M; one pass reads M once, not twice. */
/* moved */
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, int nt_allowed) {
if (m <= 0 || n <= 0 || k <= 0) return;
/* non-temporal C stores need every written 64 B line to be complete and aligned */
int ntok = nt_allowed && ((ldc & 7) == 0) && (((size_t)C & 63) == 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 */
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, ntok && !acc);
}
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;
micro(kc, ap, bp, Cp + (size_t)ir * ldc + jr, ldc, acc, mc - ir, ncc, 0);
}
}
}
}
}
}
/* moved */
/* ============ 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,1)
#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
/* moved */
/* ---- DA1: fused single-pass Strassen combine (lane_da1_mmmd24k) ----
The shipped leaf combine is 5 sequential blk_axpy2 passes: each C quadrant is read+written
2-4 times, 9 quadrant touches where 4 would do. Holding M1..M5 LIVE lets one pass produce all
four quadrants. MEASURED on the real judge, whole row, interleaved base/fuse/fuse/base at
n=2048: 1 071 333 406 / 1 071 273 606 -> 1 062 728 456 / 1 062 600 628 ticks = -0.807 %.
CORRECTNESS: comb5 OVERWRITES all four quadrants and FUSE67's two direct products are applied
AFTERWARDS. That ordering is forced: at a CHILD level the destination `C` is the parent's
recycled `M` plane, so ACCUMULATING into C11/C22 would pick up the previous product -- which is
exactly how this arm returned a WRONG ANSWER (max|diff| 2.36e+02) before the fix.
Max |diff| vs the shipped engine at n=2048 = 2.84e-14 against the problem's 3n^2*eps
= 2.79e-09 => ratio 1.0e-5. */
/* moved */
/* ==================== lane m4k_nl: NLM helpers ==================== */
/* ONE pass over the four operand quadrants produces the five Strassen sums.
base: 5 sadd calls x (2 reads + 1 write) = 15 m^2 of traffic per side.
NLM : 4 reads + 5 writes = 9 m^2 per side. */
static void sum5A(int m,const T*p11,const T*p12,const T*p21,const T*p22,int ld,
T*d1,T*d2,T*d5,T*d6,T*d7,int ldo,int lq,int nt){
nt = nt && ((((size_t)d1|(size_t)d2|(size_t)d5|(size_t)d6|(size_t)d7) & 31)==0);
for(int i=0;i<m;i++){
const T*a11=p11+(size_t)i*ld,*a12=p12+(size_t)i*ld,*a21=p21+(size_t)i*ld,*a22=p22+(size_t)i*ld;
T*r1=d1+(size_t)i*ldo,*r2=d2+(size_t)i*ldo,*r5=d5+(size_t)i*ldo,*r6=d6+(size_t)i*lq,*r7=d7+(size_t)i*lq;
int j=0;
for(;j+4<=m;j+=4){
__m256d x11=LD(a11+j),x12=LD(a12+j),x21=LD(a21+j),x22=LD(a22+j);
stn(nt,r1+j,_mm256_add_pd(x11,x22));
stn(nt,r2+j,_mm256_add_pd(x21,x22));
stn(nt,r5+j,_mm256_add_pd(x11,x12));
stn(nt,r6+j,_mm256_sub_pd(x21,x11));
stn(nt,r7+j,_mm256_sub_pd(x12,x22));
}
for(;j<m;j++){ T x11=a11[j],x12=a12[j],x21=a21[j],x22=a22[j];
r1[j]=x11+x22; r2[j]=x21+x22; r5[j]=x11+x12; r6[j]=x21-x11; r7[j]=x12-x22; }
}
}
static void sum5B(int m,const T*p11,const T*p12,const T*p21,const T*p22,int ld,
T*d1,T*d3,T*d4,T*d6,T*d7,int ldo,int lq,int nt){
nt = nt && ((((size_t)d1|(size_t)d3|(size_t)d4|(size_t)d6|(size_t)d7) & 31)==0);
for(int i=0;i<m;i++){
const T*b11=p11+(size_t)i*ld,*b12=p12+(size_t)i*ld,*b21=p21+(size_t)i*ld,*b22=p22+(size_t)i*ld;
T*r1=d1+(size_t)i*ldo,*r3=d3+(size_t)i*ldo,*r4=d4+(size_t)i*ldo,*r6=d6+(size_t)i*lq,*r7=d7+(size_t)i*lq;
int j=0;
for(;j+4<=m;j+=4){
__m256d x11=LD(b11+j),x12=LD(b12+j),x21=LD(b21+j),x22=LD(b22+j);
stn(nt,r1+j,_mm256_add_pd(x11,x22));
stn(nt,r3+j,_mm256_sub_pd(x12,x22));
stn(nt,r4+j,_mm256_sub_pd(x21,x11));
stn(nt,r6+j,_mm256_add_pd(x11,x12));
stn(nt,r7+j,_mm256_add_pd(x21,x22));
}
for(;j<m;j++){ T x11=b11[j],x12=b12[j],x21=b21[j],x22=b22[j];
r1[j]=x11+x22; r3[j]=x12-x22; r4[j]=x21-x11; r6[j]=x11+x12; r7[j]=x21+x22; }
}
}
/* SEVEN live M planes, ONE pass, all four quadrants out (base: comb5 = 9 m^2 and TWO
materialised FUSE67 products added by blk_axpy = 6 m^2 => 15 m^2). */
static void comb7(int m,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 ldm){
const int nt=g_nt;
for(int i=0;i<m;i++){
T*a=C11+(size_t)i*ldc,*b=C12+(size_t)i*ldc,*c=C21+(size_t)i*ldc,*d=C22+(size_t)i*ldc;
const T*q1=m1+(size_t)i*ldm,*q2=m2+(size_t)i*ldm,*q3=m3+(size_t)i*ldm,*q4=m4+(size_t)i*ldm,
*q5=m5+(size_t)i*ldm,*q6=m6+(size_t)i*ldm,*q7=m7+(size_t)i*ldm;
int j=0;
for(;j+4<=m;j+=4){
__m256d x1=LD(q1+j),x2=LD(q2+j),x3=LD(q3+j),x4=LD(q4+j),x5=LD(q5+j),x6=LD(q6+j),x7=LD(q7+j);
stn(nt,a+j,_mm256_add_pd(_mm256_add_pd(_mm256_sub_pd(x1,x5),x4),x7));
stn(nt,b+j,_mm256_add_pd(x3,x5));
stn(nt,c+j,_mm256_add_pd(x2,x4));
stn(nt,d+j,_mm256_add_pd(_mm256_add_pd(_mm256_sub_pd(x1,x2),x3),x6));
}
for(;j<m;j++){
double x1=q1[j],x2=q2[j],x3=q3[j],x4=q4[j],x5=q5[j],x6=q6[j],x7=q7[j];
a[j]=x1-x5+x4+x7; b[j]=x3+x5; c[j]=x2+x4; d[j]=x1-x2+x3+x6;
}
}
}
/* ================= lane m4k5_rv : MERGED LEAF, MIKRO ORDER PRESERVED =================
*
* WHAT THE BASE DOES AT A LEAF (m = 512, KC = 512, NC = 1024 -> nc = m):
* seven gemm_core(m,m,m) calls, each with its own packA and packB.
* The seven products read their A operand from 12 (a1,a2) quadrant strips and their
* B operand from 12 (b1,b2) quadrant panels -- while only FOUR quadrants per side
* exist. That is the 12-uses-of-4-quadrants multiplicity.
*
* WHAT THIS ARM DOES:
* packB7 reads the four B quadrants ONCE and emits SEVEN packed B panels.
* packA7 reads the four A quadrant strips ONCE and emits SEVEN packed A strips.
* comb7 consumes the seven live F planes in one pass (the base needs comb5 + 2
* materialised FUSE67 products through blk_axpy = 17 m^2 of touches against 11).
*
* ** THE ORDERING, WHICH IS THE WHOLE POINT AND THE THING m4k_mg's +13.04 % DID NOT DO:
* m4k_mg emitted `for (jr) { for (g) { micro } }` -- seven A panels and seven B slices
* live at once against a 32 KB L1, and the mikro paid +64 %. Here the product loop is
* OUTSIDE the jr loop: `for (g) { for (jr) { micro } }`. During product g's sweep the
* only packed A strip being read is Ap[g] (24 KB, L1-resident across all 64 jr blocks)
* and the only packed B panel being streamed is Bp[g] (2 MB, L3-resident) -- i.e. the
* mikro's access pattern is the base's. The mikro is at 100 % of its FMA port floor
* (5.947 ticks/k-step x 1.0089 = 6.00 cycles), so it must not be perturbed.
*
* SOURCE-READ LEDGER per leaf per side (m^2 = 262144 doubles = 2.097 MB):
* base 12 m^2 = 25.17 MB of quadrant reads arm 4 m^2 = 8.39 MB
* PRICE PAID: seven packed B panels (7 * KC * m = 14.7 MB) live instead of one 2 MB
* panel re-used seven times. That is MEASURED here, not argued.
* ============================================================================ */
#define M5_MC 6
#define M5_NPL 7 /* the seven Strassen products */
#ifndef M5_NC
#define M5_NC 128 /* the B column block (m/K): the residency knob */
#endif
#ifndef M5_MODE
#define M5_MODE 0 /* 0 = merge BOTH sides | 1 = merge B ONLY, product-outer */
#endif
#define M5_BPSTRIDE (((512 + NR - 1) / NR) * NR * KC)
/* The seven leaf products, in the order the arm emits them (comb7's x1..x7):
* g0 = M1 = (A11+A22)(B11+B22) g4 = M5 = (A11+A12) B22
* g1 = M2 = (A21+A22) B11 g5 = M6 = (A21-A11)(B11+B12)
* g2 = M3 = A11 (B12-B22) g6 = M7 = (A12-A22)(B21+B22)
* g3 = M4 = A22 (B21-B11)
* comb7 forms C11 = x1-x5+x4+x7, C12 = x3+x5, C21 = x2+x4, C22 = x1-x2+x3+x6
* which is EXACTLY comb5 + the two FUSE67 products the base runs afterwards.
*/
static T M5_Ap[7 * 6 * KCP + 64] __attribute__((aligned(64)));
static T M5_Bp[7 * M5_BPSTRIDE + 512] __attribute__((aligned(64)));
static unsigned long long g_l7_calls = 0, g_l7_pb = 0, g_l7_pa = 0;
static int g_ph = 0; /* 1 skip packB7 | 2 skip packA7 | 4 skip mikro | 8 skip comb7 | 16 force base leaf */
static const int g_nc = M5_NC;
/* ONE pass over the four A quadrant strips -> SEVEN packed A strips.
A is the leaf's 2m x 2m operand (stride lda); the strip is rows ic..ic+mc-1 of
each quadrant, kc = m columns. Exactly the 12 -> 4 quadrant-read reduction. */
static void packA7(int mc, int kc, const T *A, int lda, int m, T **Ap) {
for (int r = 0; r < mc; r++) {
const T *s11 = A + (size_t)r * lda;
const T *s12 = s11 + m;
const T *s21 = s11 + (size_t)m * lda;
const T *s22 = s21 + m;
T *o0 = Ap[0] + (size_t)r * KCP, *o1 = Ap[1] + (size_t)r * KCP,
*o2 = Ap[2] + (size_t)r * KCP, *o3 = Ap[3] + (size_t)r * KCP,
*o4 = Ap[4] + (size_t)r * KCP, *o5 = Ap[5] + (size_t)r * KCP,
*o6 = Ap[6] + (size_t)r * KCP;
int k = 0;
for (; k + 8 <= kc; k += 8) {
__m256d l11 = LD(s11 + k), h11 = LD(s11 + k + 4);
__m256d l12 = LD(s12 + k), h12 = LD(s12 + k + 4);
__m256d l21 = LD(s21 + k), h21 = LD(s21 + k + 4);
__m256d l22 = LD(s22 + k), h22 = LD(s22 + k + 4);
ST(o0 + k, _mm256_add_pd(l11, l22)); ST(o0 + k + 4, _mm256_add_pd(h11, h22));
ST(o1 + k, _mm256_add_pd(l21, l22)); ST(o1 + k + 4, _mm256_add_pd(h21, h22));
ST(o2 + k, l11); ST(o2 + k + 4, h11);
ST(o3 + k, l22); ST(o3 + k + 4, h22);
ST(o4 + k, _mm256_add_pd(l11, l12)); ST(o4 + k + 4, _mm256_add_pd(h11, h12));
ST(o5 + k, _mm256_sub_pd(l21, l11)); ST(o5 + k + 4, _mm256_sub_pd(h21, h11));
ST(o6 + k, _mm256_sub_pd(l12, l22)); ST(o6 + k + 4, _mm256_sub_pd(h12, h22));
}
for (; k < kc; k++) {
T a11 = s11[k], a12 = s12[k], a21 = s21[k], a22 = s22[k];
o0[k] = a11 + a22; o1[k] = a21 + a22; o2[k] = a11; o3[k] = a22;
o4[k] = a11 + a12; o5[k] = a21 - a11; o6[k] = a12 - a22;
}
}
}
/* ONE pass over the four B quadrants -> SEVEN packed B panels.
B is the leaf's 2m x 2m operand (stride ldb); kc = m rows, nc = m columns.
Source reads: 4 * kc * nc doubles. Base: 12 * kc * nc (7 products x their 1 or 2
quadrant sources). The panel k index and the quadrant row index coincide. */
static void packB7(int kc, int nc, const T *B, int ldb, int m, T **Bp) {
/* E8SWAP ORDER (the base's own packB carries the comment for this; it is board-priced at
-2.21 % on the sibling row mmmd2k). The k-inner form reads 64 B per ldb*8 = 8 KB jump
and the HW prefetcher cannot follow it -- every line is a demand miss. k0 OUTER +
column-group INNER makes consecutive reads walk a ROW CONTIGUOUSLY, four streams. */
for (int k0 = 0; k0 < kc; k0 += 8) {
int kmax = kc - k0 < 8 ? kc - k0 : 8;
for (int j0 = 0; j0 < nc; j0 += NR) {
T *o0 = Bp[0] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
T *o1 = Bp[1] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
T *o2 = Bp[2] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
T *o3 = Bp[3] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
T *o4 = Bp[4] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
T *o5 = Bp[5] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
T *o6 = Bp[6] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
for (int kk = 0; kk < kmax; kk++) {
const int k = k0 + kk;
const T *b11 = B + (size_t)k*ldb + j0, *b12 = b11 + m;
const T *b21 = b11 + (size_t)m*ldb, *b22 = b21 + m;
T *r0 = o0 + (size_t)kk*NR, *r1 = o1 + (size_t)kk*NR, *r2 = o2 + (size_t)kk*NR,
*r3 = o3 + (size_t)kk*NR, *r4 = o4 + (size_t)kk*NR, *r5 = o5 + (size_t)kk*NR,
*r6 = o6 + (size_t)kk*NR;
if (j0 + NR <= nc) {
__m256d l11 = LD(b11), h11 = LD(b11+4), l12 = LD(b12), h12 = LD(b12+4);
__m256d l21 = LD(b21), h21 = LD(b21+4), l22 = LD(b22), h22 = LD(b22+4);
ST(r0, _mm256_add_pd(l11, l22)); ST(r0+4, _mm256_add_pd(h11, h22));
ST(r1, l11); ST(r1+4, h11);
ST(r2, _mm256_sub_pd(l12, l22)); ST(r2+4, _mm256_sub_pd(h12, h22));
ST(r3, _mm256_sub_pd(l21, l11)); ST(r3+4, _mm256_sub_pd(h21, h11));
ST(r4, l22); ST(r4+4, h22);
ST(r5, _mm256_add_pd(l11, l12)); ST(r5+4, _mm256_add_pd(h11, h12));
ST(r6, _mm256_add_pd(l21, l22)); ST(r6+4, _mm256_add_pd(h21, h22));
} else {
for (int j = 0; j < NR; j++) {
T x11 = b11[j], x12 = b12[j], x21 = b21[j], x22 = b22[j];
r0[j] = x11+x22; r1[j] = x11; r2[j] = x12-x22; r3[j] = x21-x11;
r4[j] = x22; r5[j] = x11+x12; r6[j] = x21+x22;
}
}
}
}
}
}
static void leaf7b(int m, const T *A, int lda, const T *B, int ldb, T *C, int ldc, T *F0) {
const int ld7 = 7 * m;
T *F[7];
for (int g = 0; g < M5_NPL; g++) F[g] = F0 + (size_t)g * m;
T *Bp[7];
for (int g = 0; g < M5_NPL; g++) Bp[g] = M5_Bp + (size_t)g * M5_BPSTRIDE;
const int nt = ((ldc & 7) == 0) && (((size_t)C & 63) == 0) &&
((ld7 & 7) == 0) && (((size_t)F0 & 63) == 0);
const T *A11 = A, *A12 = A + m, *A21 = A + (size_t)m * lda, *A22 = A21 + m;
T *C11 = C, *C12 = C + m, *C21 = C + (size_t)m * ldc, *C22 = C + (size_t)m * ldc + m;
const T *pa1[7] = {A11, A21, A11, A22, A11, A21, A12};
const T *pa2[7] = {A22, A22, NIL, NIL, A12, A11, A22};
const int poa[7] = { 1, 1, 0, 0, 1, -1, -1};
g_l7_calls++;
g_l7_pb++;
if (!(g_ph & 1)) packB7(m, m, B, ldb, m, Bp); /* ONE pass: 4 B quadrants -> 7 panels */
for (int g = 0; g < M5_NPL; g++) {
for (int ic = 0; ic < m; ic += M5_MC) {
int mc = m - ic < M5_MC ? m - ic : M5_MC;
g_l7_pa++;
if (!(g_ph & 2))
packA(mc, m, pa1[g] + (size_t)ic * lda,
pa2[g] ? pa2[g] + (size_t)ic * lda : (const T *)0, poa[g], lda, M5_Ap);
if (g_ph & 4) continue;
for (int jr = 0; jr < m; jr += NR) {
int ncc = m - jr < NR ? m - jr : NR;
for (int ir = 0; ir + MR <= mc; ir += MR)
micro(m, M5_Ap + (size_t)(ir / MR) * MR * KCP, Bp[g] + (size_t)(jr / NR) * KC * NR,
F[g] + (size_t)(ic + ir) * ld7 + jr, ld7, 0, MR, ncc, nt);
int ir = mc - (mc % MR);
if (ir < mc)
micro(m, M5_Ap + (size_t)(ir / MR) * MR * KCP, Bp[g] + (size_t)(jr / NR) * KC * NR,
F[g] + (size_t)(ic + ir) * ld7 + jr, ld7, 0, mc - ir, ncc, 0);
}
}
}
if (!(g_ph & 8)) comb7(m, C11, C12, C21, C22, ldc, F[0], F[1], F[2], F[3], F[4], F[5], F[6], ld7);
}
static void leaf7(int m, const T *A, int lda, const T *B, int ldb, T *C, int ldc, T *F0) {
const int ld7 = 7 * m;
T *F[7];
for (int g = 0; g < M5_NPL; g++) F[g] = F0 + (size_t)g * m;
T *Ap[7], *Bp[7];
for (int g = 0; g < M5_NPL; g++) {
Ap[g] = M5_Ap + (size_t)g * 6 * KCP;
Bp[g] = M5_Bp + (size_t)g * M5_BPSTRIDE;
}
const int nt = ((ldc & 7) == 0) && (((size_t)C & 63) == 0) &&
((ld7 & 7) == 0) && (((size_t)F0 & 63) == 0);
T *C11 = C, *C12 = C + m, *C21 = C + (size_t)m * ldc, *C22 = C + (size_t)m * ldc + m;
g_l7_calls++;
/* ---------------------------------------------------------------------------
★ THE COLUMN BLOCK IS THE RESIDENCY KNOB, AND IT IS NOT OPTIONAL.
The mikro re-reads the whole packed B set once per `ic` (86 times per leaf).
With NC7 = m that set is 7 * KC * m = 14.7 MB against a 6 MB L3: the re-read
goes to DRAM and the mikro -- which has ZERO issue slack (12 vfmadd231pd on
2 ports = 6 cycles = exactly its 12 memory uops on 2 load ports) -- pays it
directly. Measured locally: 7 x 2 MB panels = +74 % on the leaf.
NC7 = m/K makes the set 7*KC*m/K, and the B SOURCE read volume is INVARIANT
at 4 m^2 per side per leaf for every K (each jc pass reads only its own
quadrant columns). The A source read is 4*K m^2: the price of K.
--------------------------------------------------------------------------- */
for (int jc = 0; jc < m; jc += M5_NC) {
int nc = m - jc < M5_NC ? m - jc : M5_NC;
g_l7_pb++;
if (!(g_ph & 1)) packB7(m, nc, B + jc, ldb, m, Bp); /* * 4 B quadrant reads, 7 panels */
for (int ic = 0; ic < m; ic += M5_MC) {
int mc = m - ic < M5_MC ? m - ic : M5_MC;
g_l7_pa++;
if (!(g_ph & 2)) packA7(mc, m, A + (size_t)ic * lda, lda, m, Ap); /* * 4 A strip reads, 7 strips */
if (g_ph & 4) continue;
/* PRODUCT OUTERMOST -- the mikro sees one A strip and one B panel, as in the base */
for (int g = 0; g < M5_NPL; g++) {
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) {
micro(m, Ap[g] + (size_t)(ir / MR) * MR * KCP, Bp[g] + (size_t)(jr / NR) * KC * NR,
F[g] + (size_t)(ic + ir) * ld7 + jc + jr, ld7, 0, MR, ncc, nt);
}
int ir = mc - (mc % MR);
if (ir < mc) {
micro(m, Ap[g] + (size_t)(ir / MR) * MR * KCP, Bp[g] + (size_t)(jr / NR) * KC * NR,
F[g] + (size_t)(ic + ir) * ld7 + jc + jr, ld7, 0, mc - ir, ncc, 0);
}
}
}
}
}
if (!(g_ph & 8)) comb7(m, C11, C12, C21, C22, ldc, F[0], F[1], F[2], F[3], F[4], F[5], F[6], ld7);
}
#if M5_MODE == 0
#define M5_LEAF(m,A,lda,B,ldb,C,ldc,sc) leaf7 (m,A,lda,B,ldb,C,ldc,sc)
#elif M5_MODE == 1
#define M5_LEAF(m,A,lda,B,ldb,C,ldc,sc) leaf7b(m,A,lda,B,ldb,C,ldc,sc)
#elif M5_MODE == 2
#define M5_LEAF(m,A,lda,B,ldb,C,ldc,sc) leaf7g(m,A,lda,B,ldb,C,ldc,sc)
#elif M5_MODE == 3
#define M5_LEAF(m,A,lda,B,ldb,C,ldc,sc) leaf7p(m,A,lda,B,ldb,C,ldc,sc)
#endif
/* ============================================================================
* ★★★ THE GROUP MERGE (mode 2 / mode 3) -- THE FORM WITH **NO** MATERIALISATION PENALTY.
*
* leaf7 (all seven panels live) measured +66/+49/+13/+10 % across nc = 512/256/128/64.
* leaf7b (product-outer, 7 panels) measured +3.5 %.
* BOTH LOSE FOR THE SAME REASON, AND IT IS NOT THE MIKRO'S L1 FOOTPRINT:
* a merged pass produces SEVEN panels at once, so 14.7 MB (nc = 512) of packed
* output must be written and then re-read. The base writes the same bytes into
* ONE reused 2 MB buffer that stays L3-resident; the merged form writes them to
* seven buffers that do not, so 16.8 MB of saved source reads are paid back as
* ~29 MB of DRAM write + DRAM read. Measured exchange: mode 1 = +1.5 ms/leaf.
*
* THE FIX IS TO MERGE IN GROUPS SMALL ENOUGH THAT THE GROUP'S PANEL SET FITS L3.
* OPTIMAL 3-GROUP PARTITION (union sizes, quadrant index 0..3 = 11/12/21/22):
* G0 = {M1, M3, M6}: A needs {0,2,3} (3 of 5 base reads) B needs {0,1,3} (3 of 6)
* G1 = {M2, M4, M7}: A needs {1,2,3} (3 of 5) B needs {0,2,3} (3 of 5)
* G2 = {M5} : A needs {0,1} (2 of 2) B needs {3} (1 of 1)
* => SOURCE READS A 12 -> 8 m^2 , B 12 -> 7 m^2 per leaf per side
* => PANEL SET PER GROUP <= 3 * KC * m = 6.0 MB (vs 14.7 MB), i.e. L3-sized.
* ========================================================================== */
#define M5_PBINIT \
const T *b11 = B + (size_t)k * ldb + j0, *b12 = b11 + m; \
const T *b21 = b11 + (size_t)m * ldb, *b22 = b21 + m; \
T *r0 = o0 + (size_t)k * NR, *r1 = o1 + (size_t)k * NR, *r2 = o2 + (size_t)k * NR; \
(void)r0; (void)r1; (void)r2
#define M5_PBTGT \
T *o0 = Bp[0] + (size_t)(j0 / NR) * KC * NR; \
T *o1 = Bp[1] + (size_t)(j0 / NR) * KC * NR; \
T *o2 = Bp[2] + (size_t)(j0 / NR) * KC * NR
/* ---- G0 = {M1, M3, M6} : B needs q0, q1, q3 ---- */
static void packBg0(int kc, int nc, const T *B, int ldb, int m, T **Bp) {
for (int k0 = 0; k0 < kc; k0 += 8) {
int kmax = kc - k0 < 8 ? kc - k0 : 8;
for (int j0 = 0; j0 < nc; j0 += NR) {
T *o0 = Bp[0] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
T *o1 = Bp[1] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
T *o2 = Bp[2] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
for (int kk = 0; kk < kmax; kk++) {
const int k = k0 + kk;
const T *b11 = B + (size_t)k*ldb + j0, *b12 = b11 + m;
const T *b21 = b11 + (size_t)m*ldb, *b22 = b21 + m;
T *r0 = o0 + (size_t)kk*NR, *r1 = o1 + (size_t)kk*NR, *r2 = o2 + (size_t)kk*NR;
(void)r1; (void)r2;
__m256d x0 = LD(b11), x0h = LD(b11+4), x1 = LD(b12), x1h = LD(b12+4);
__m256d x3 = LD(b22), x3h = LD(b22+4);
ST(r0, _mm256_add_pd(x0, x3)); ST(r0+4, _mm256_add_pd(x0h, x3h)); /* M1 = B11+B22 */
ST(r1, _mm256_sub_pd(x1, x3)); ST(r1+4, _mm256_sub_pd(x1h, x3h)); /* M3 = B12-B22 */
ST(r2, _mm256_add_pd(x0, x1)); ST(r2+4, _mm256_add_pd(x0h, x1h)); /* M6 = B11+B12 */
}
}
}
}
static void packBg1(int kc, int nc, const T *B, int ldb, int m, T **Bp) {
for (int k0 = 0; k0 < kc; k0 += 8) {
int kmax = kc - k0 < 8 ? kc - k0 : 8;
for (int j0 = 0; j0 < nc; j0 += NR) {
T *o0 = Bp[0] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
T *o1 = Bp[1] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
T *o2 = Bp[2] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
for (int kk = 0; kk < kmax; kk++) {
const int k = k0 + kk;
const T *b11 = B + (size_t)k*ldb + j0, *b12 = b11 + m;
const T *b21 = b11 + (size_t)m*ldb, *b22 = b21 + m;
T *r0 = o0 + (size_t)kk*NR, *r1 = o1 + (size_t)kk*NR, *r2 = o2 + (size_t)kk*NR;
(void)r1; (void)r2;
__m256d x0 = LD(b11), x0h = LD(b11+4), x2 = LD(b21), x2h = LD(b21+4);
__m256d x3 = LD(b22), x3h = LD(b22+4);
ST(r0, x0); ST(r0+4, x0h); /* M2 = B11 */
ST(r1, _mm256_sub_pd(x2, x0)); ST(r1+4, _mm256_sub_pd(x2h, x0h)); /* M4 = B21-B11 */
ST(r2, _mm256_add_pd(x2, x3)); ST(r2+4, _mm256_add_pd(x2h, x3h)); /* M7 = B21+B22 */
}
}
}
}
static void packBg2(int kc, int nc, const T *B, int ldb, int m, T **Bp) {
for (int k0 = 0; k0 < kc; k0 += 8) {
int kmax = kc - k0 < 8 ? kc - k0 : 8;
for (int j0 = 0; j0 < nc; j0 += NR) {
T *o0 = Bp[0] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
T *o1 = Bp[1] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
T *o2 = Bp[2] + (size_t)(j0/NR)*KC*NR + (size_t)k0*NR;
for (int kk = 0; kk < kmax; kk++) {
const int k = k0 + kk;
const T *b11 = B + (size_t)k*ldb + j0, *b12 = b11 + m;
const T *b21 = b11 + (size_t)m*ldb, *b22 = b21 + m;
T *r0 = o0 + (size_t)kk*NR, *r1 = o1 + (size_t)kk*NR, *r2 = o2 + (size_t)kk*NR;
(void)r1; (void)r2;
__m256d x3 = LD(b22), x3h = LD(b22+4);
ST(r0, x3); ST(r0+4, x3h); /* M5 = B22 */
}
}
}
}
/* ---- A groups: same shape; the strip is rows ic..ic+mc-1 of each 2m x 2m operand ---- */
#define M5_PAINIT \
const T *s11 = A + (size_t)r * lda; \
const T *s12 = s11 + m; \
const T *s21 = s11 + (size_t)m * lda; \
const T *s22 = s21 + m; \
T *o0 = Ap[0] + (size_t)r * KCP, *o1 = Ap[1] + (size_t)r * KCP, \
*o2 = Ap[2] + (size_t)r * KCP; \
int k = 0; \
for (; k + 8 <= kc; k += 8)
#define M5_PATAIL \
for (; k < kc; k++)
static void packAg0(int mc, int kc, const T *A, int lda, int m, T **Ap) {
for (int r = 0; r < mc; r++) { M5_PAINIT {
__m256d x0 = LD(s11 + k), x0h = LD(s11 + k + 4);
__m256d x2 = LD(s21 + k), x2h = LD(s21 + k + 4);
__m256d x3 = LD(s22 + k), x3h = LD(s22 + k + 4);
ST(o0 + k, _mm256_add_pd(x0, x3)); ST(o0 + k + 4, _mm256_add_pd(x0h, x3h)); /* M1 = A11+A22 */
ST(o1 + k, x0); ST(o1 + k + 4, x0h); /* M3 = A11 */
ST(o2 + k, _mm256_sub_pd(x2, x0)); ST(o2 + k + 4, _mm256_sub_pd(x2h, x0h)); /* M6 = A21-A11 */
}
M5_PATAIL { T x0 = s11[k], x2 = s21[k], x3 = s22[k];
o0[k] = x0 + x3; o1[k] = x0; o2[k] = x2 - x0; } }
}
/* ---- G1 A = {M2, M4, M7} : needs q1, q2, q3 ---- */
static void packAg1(int mc, int kc, const T *A, int lda, int m, T **Ap) {
for (int r = 0; r < mc; r++) { M5_PAINIT {
__m256d x1 = LD(s12 + k), x1h = LD(s12 + k + 4);
__m256d x2 = LD(s21 + k), x2h = LD(s21 + k + 4);
__m256d x3 = LD(s22 + k), x3h = LD(s22 + k + 4);
ST(o0 + k, _mm256_add_pd(x2, x3)); ST(o0 + k + 4, _mm256_add_pd(x2h, x3h)); /* M2 = A21+A22 */
ST(o1 + k, x3); ST(o1 + k + 4, x3h); /* M4 = A22 */
ST(o2 + k, _mm256_sub_pd(x1, x3)); ST(o2 + k + 4, _mm256_sub_pd(x1h, x3h)); /* M7 = A12-A22 */
}
M5_PATAIL { T x1 = s12[k], x2 = s21[k], x3 = s22[k];
o0[k] = x2 + x3; o1[k] = x3; o2[k] = x1 - x3; } }
}
/* ---- G2 A = {M5} : needs q0, q1 ---- */
static void packAg2(int mc, int kc, const T *A, int lda, int m, T **Ap) {
for (int r = 0; r < mc; r++) { M5_PAINIT {
__m256d x0 = LD(s11 + k), x0h = LD(s11 + k + 4);
__m256d x1 = LD(s12 + k), x1h = LD(s12 + k + 4);
ST(o0 + k, _mm256_add_pd(x0, x1)); ST(o0 + k + 4, _mm256_add_pd(x0h, x1h)); /* M5 = A11+A12 */
}
M5_PATAIL { T x0 = s11[k], x1 = s12[k]; o2[k] = x0 + x1; } }
}
/* THE GROUP MERGE: one merged B pass and (per ic) one merged A pass PER GROUP,
panels <= 3 * KC * m = 6 MB, mikro order product-outer inside the group. */
static void leaf7g(int m, const T *A, int lda, const T *B, int ldb, T *C, int ldc, T *F0) {
const int ld7 = 7 * m;
T *F[7];
for (int g = 0; g < M5_NPL; g++) F[g] = F0 + (size_t)g * m;
const int nt = ((ldc & 7) == 0) && (((size_t)C & 63) == 0) &&
((ld7 & 7) == 0) && (((size_t)F0 & 63) == 0);
T *C11 = C, *C12 = C + m, *C21 = C + (size_t)m * ldc, *C22 = C + (size_t)m * ldc + m;
/* the three groups, by GLOBAL product id */
static const int gid[3][3] = {{0, 2, 5}, {1, 3, 6}, {4, -1, -1}};
const int gsz[3] = {3, 3, 1};
g_l7_calls++;
for (int G = 0; G < 3; G++) {
T *Bp[3], *Ap[3];
for (int l = 0; l < 3; l++) {
Bp[l] = M5_Bp + (size_t)l * M5_BPSTRIDE;
Ap[l] = M5_Ap + (size_t)l * 6 * KCP;
}
g_l7_pb++;
if (!(g_ph & 1)) {
if (G == 0) packBg0(m, m, B, ldb, m, Bp);
else if (G == 1) packBg1(m, m, B, ldb, m, Bp);
else packBg2(m, m, B, ldb, m, Bp);
}
for (int ic = 0; ic < m; ic += M5_MC) {
int mc = m - ic < M5_MC ? m - ic : M5_MC;
g_l7_pa++;
if (!(g_ph & 2)) {
if (G == 0) packAg0(mc, m, A + (size_t)ic * lda, lda, m, Ap);
else if (G == 1) packAg1(mc, m, A + (size_t)ic * lda, lda, m, Ap);
else packAg2(mc, m, A + (size_t)ic * lda, lda, m, Ap);
}
if (g_ph & 4) continue;
for (int l = 0; l < gsz[G]; l++) {
const int P = gid[G][l];
for (int jr = 0; jr < m; jr += NR) {
int ncc = m - jr < NR ? m - jr : NR;
for (int ir = 0; ir + MR <= mc; ir += MR)
micro(m, Ap[l] + (size_t)(ir / MR) * MR * KCP, Bp[l] + (size_t)(jr / NR) * KC * NR,
F[P] + (size_t)(ic + ir) * ld7 + jr, ld7, 0, MR, ncc, nt);
int ir = mc - (mc % MR);
if (ir < mc)
micro(m, Ap[l] + (size_t)(ir / MR) * MR * KCP, Bp[l] + (size_t)(jr / NR) * KC * NR,
F[P] + (size_t)(ic + ir) * ld7 + jr, ld7, 0, mc - ir, ncc, 0);
}
}
}
}
if (!(g_ph & 8)) comb7(m, C11, C12, C21, C22, ldc, F[0], F[1], F[2], F[3], F[4], F[5], F[6], ld7);
}
/* ============================================================================
* leaf7p -- GROUP-MERGED **B**, PRODUCT-OUTER. mode 3.
*
* THIS IS THE ONLY FORM THAT KEEPS BOTH PROPERTIES AT ONCE:
* (i) the B source is read ONCE PER GROUP (union of the group's quadrants:
* 3+3+1 = 7 m^2 instead of 12), and
* (ii) only ONE packed panel is hot at a time -- product-outer -- so the mikro's
* B re-read working set is 2 MB, exactly the base's, and the group's 6 MB of
* panels never has to be co-resident.
* The A side is NOT merged (per-product packA into the base's own 24 KB strip
* buffer): the A merge requires the ic-outer order, and ic-outer is what broke
* leaf7 (nc sweep) and leaf7g (+45 %), both times because the mikro then cycles
* the whole panel set once per ic.
* ========================================================================== */
static void leaf7p(int m, const T *A, int lda, const T *B, int ldb, T *C, int ldc, T *F0) {
const int ld7 = 7 * m;
T *F[7];
for (int g = 0; g < M5_NPL; g++) F[g] = F0 + (size_t)g * m;
T *Bp[3];
for (int l = 0; l < 3; l++) Bp[l] = M5_Bp + (size_t)l * M5_BPSTRIDE;
const int nt = ((ldc & 7) == 0) && (((size_t)C & 63) == 0) &&
((ld7 & 7) == 0) && (((size_t)F0 & 63) == 0);
const T *A11 = A, *A12 = A + m, *A21 = A + (size_t)m * lda, *A22 = A21 + m;
T *C11 = C, *C12 = C + m, *C21 = C + (size_t)m * ldc, *C22 = C + (size_t)m * ldc + m;
const T *pa1[7] = {A11, A21, A11, A22, A11, A21, A12};
const T *pa2[7] = {A22, A22, NIL, NIL, A12, A11, A22};
const int poa[7] = { 1, 1, 0, 0, 1, -1, -1};
static const int gid[3][3] = {{0, 2, 5}, {1, 3, 6}, {4, -1, -1}};
const int gsz[3] = {3, 3, 1};
g_l7_calls++;
for (int G = 0; G < 3; G++) {
g_l7_pb++;
if (!(g_ph & 1)) {
if (G == 0) packBg0(m, m, B, ldb, m, Bp);
else if (G == 1) packBg1(m, m, B, ldb, m, Bp);
else packBg2(m, m, B, ldb, m, Bp);
}
for (int l = 0; l < gsz[G]; l++) {
const int P = gid[G][l];
for (int ic = 0; ic < m; ic += M5_MC) {
int mc = m - ic < M5_MC ? m - ic : M5_MC;
g_l7_pa++;
if (!(g_ph & 2))
packA(mc, m, pa1[P] + (size_t)ic * lda,
pa2[P] ? pa2[P] + (size_t)ic * lda : (const T *)0, poa[P], lda, M5_Ap);
if (g_ph & 4) continue;
for (int jr = 0; jr < m; jr += NR) {
int ncc = m - jr < NR ? m - jr : NR;
for (int ir = 0; ir + MR <= mc; ir += MR)
micro(m, M5_Ap + (size_t)(ir / MR) * MR * KCP, Bp[l] + (size_t)(jr / NR) * KC * NR,
F[P] + (size_t)(ic + ir) * ld7 + jr, ld7, 0, MR, ncc, nt);
int ir = mc - (mc % MR);
if (ir < mc)
micro(m, M5_Ap + (size_t)(ir / MR) * MR * KCP, Bp[l] + (size_t)(jr / NR) * KC * NR,
F[P] + (size_t)(ic + ir) * ld7 + jr, ld7, 0, mc - ir, ncc, 0);
}
}
}
}
if (!(g_ph & 8)) comb7(m, C11, C12, C21, C22, ldc, F[0], F[1], F[2], F[3], F[4], F[5], F[6], ld7);
}
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;
#if FUSE_NL
T *TA = leaf ? (T *)0 : scratch + 5 * (size_t)m * m;
T *TB = leaf ? (T *)0 : scratch + 6 * (size_t)m * m;
T *sub = scratch + (leaf ? (size_t)m * m : 7 * (size_t)m * m);
#else
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);
#endif
/* 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)
if (leaf) {
if (m == 512 && KC == 512 && NR == 8 && MR == 6) {
if (!(g_ph & 16)) { M5_LEAF(m, A, lda, B, ldb, C, ldc, scratch); return; }
}
/* DA1 (FALLBACK for m != 512, or when the driver forces the base leaf:
this path is LIVE, not dead code): */
#if ILV
/* m211 ILV: the five live M planes INTERLEAVED by row (plane i at base+i*m, row stride
5*m). comb5's five read streams become ONE contiguous 5*m-double stream per output
row -- one prefetch stream, one page set -- while the mikro's per-row NT stores only
change their stride. BIT-IDENTICAL (same values, same order, addresses only). */
const int ld5 = 5 * m;
T *M1 = scratch, *M2 = M1 + m, *M3 = M2 + m, *M4 = M3 + m, *M5 = M4 + m;
GC(A11, A12, 1, lda, B22, NIL, 0, ldb, M5, ld5, 0);
GC(A11, NIL, 0, lda, B12, B22, 2, ldb, M3, ld5, 0);
GC(A11, A22, 1, lda, B11, B22, 1, ldb, M1, ld5, 0);
GC(A21, A22, 1, lda, B11, NIL, 0, ldb, M2, ld5, 0);
GC(A22, NIL, 0, lda, B21, B11, 2, ldb, M4, ld5, 0);
comb5(m, C11, C12, C21, C22, ldc, M1, M2, M3, M4, M5, ld5);
#else
T *M1 = scratch, *M2 = M1 + (size_t)m * m, *M3 = M2 + (size_t)m * m,
*M4 = M3 + (size_t)m * m, *M5 = M4 + (size_t)m * m;
GC(A11, A22, 1, lda, B11, B22, 1, ldb, M1, m, 0);
GC(A21, A22, 1, lda, B11, NIL, 0, ldb, M2, m, 0);
GC(A11, NIL, 0, lda, B12, B22, 2, ldb, M3, m, 0);
GC(A22, NIL, 0, lda, B21, B11, 2, ldb, M4, m, 0);
GC(A11, A12, 1, lda, B22, NIL, 0, ldb, M5, m, 0);
comb5(m, C11, C12, C21, C22, ldc, M1, M2, M3, M4, M5, m);
#endif
PROD(-1, A21, A11, +1, B11, B12, C22, ldc, 1);
PROD(-1, A12, A22, +1, B21, B22, C11, ldc, 1);
return;
}
{
/* ---------- lane m4k_nl MS5: THE C-DEAD WINDOW AT EVERY NON-LEAF NODE ----------
A node's own C is dead from its entry until its own combine (z4 SS1), and with a
7-input combine there are NO post-combine products, so nothing pool-resident needs
that window. C is n^2 = 4*m^2 doubles = EXACTLY four m^2 slots. The four slots
used are s3/s4/s8/s9 -- the sums consumed by M6/M7, whose OUTPUTS land in pool
slots -- so comb7 never reads a C-resident plane while it writes C.
=> 7 slots from the pool per non-leaf level, i.e. THE BASE'S OWN pool size. */
int istop = g_top; g_top = 0;
if (istop) {
/* ---------- lane m4k_rst NLT: the top node's ORIGINAL pre-FUSE_NL branch ---------- */
const int ld5 = 5 * m;
T *F1 = scratch, *F2 = F1 + m, *F3 = F2 + m, *F4 = F3 + m, *F5 = F4 + m;
PROD(+1, A11, A22, +1, B11, B22, F1, ld5, 0);
PROD(+1, A21, A22, 0, B11, NIL, F2, ld5, 0);
PROD( 0, A11, NIL, -1, B12, B22, F3, ld5, 0);
PROD( 0, A22, NIL, -1, B21, B11, F4, ld5, 0);
PROD(+1, A11, A12, 0, B22, NIL, F5, ld5, 0);
comb5(m, C11, C12, C21, C22, ldc, F1, F2, F3, F4, F5, ld5);
PROD(-1, A21, A11, +1, B11, B12, F1, m, 0); blk_axpy(m, m, C22, ldc, F1, m, +1, 1);
PROD(-1, A12, A22, +1, B21, B22, F2, m, 0); blk_axpy(m, m, C11, ldc, F2, m, +1, 1);
return;
}
g_nlm_nodes++;
const int ldo = m;
const size_t SP = (size_t)m*m + 8;
T *s3, *s4, *s8, *s9, *subn;
T *s0 = scratch + 0*SP, *s1 = scratch + 1*SP, *s2 = scratch + 2*SP;
T *s5 = scratch + 3*SP, *s6 = scratch + 4*SP, *s7 = scratch + 5*SP, *sX = scratch + 6*SP;
/* m4k_rst: the four C slots ARE the four quadrant buffers, at the C row stride ldc.
They exactly tile the node's own C region for ANY parent plane layout, and they
are the addresses comb7 writes -- so s3/s4/s8/s9 are dead exactly when it overwrites
them. Identical to the shipped placement when ldc == 2*m (an MS5-form parent). */
const int lq = ldc;
s3 = C; s4 = C + (size_t)m; s8 = C + (size_t)m*ldc; s9 = C + (size_t)m*ldc + m;
subn = scratch + 7*SP;
const int nts = g_nt;
sum5A(m, A11, A12, A21, A22, lda, s0, s1, s2, s3, s4, ldo, lq, nts);
sum5B(m, B11, B12, B21, B22, ldb, s5, s6, s7, s8, s9, ldo, lq, nts);
dsgemm(m, s0, s5, sX, ldo, ldo, m, subn); /* M1 = (A11+A22)(B11+B22) -> sX */
dsgemm(m, s1, B11, s0, ldo, ldb, m, subn); /* M2 = (A21+A22) B11 -> s0 */
dsgemm(m, A11, s6, s5, lda, ldo, m, subn); /* M3 = A11 (B12-B22) -> s5 */
dsgemm(m, A22, s7, s1, lda, ldo, m, subn); /* M4 = A22 (B21-B11) -> s1 */
dsgemm(m, s2, B22, s6, ldo, ldb, m, subn); /* M5 = (A11+A12) B22 -> s6 */
dsgemm(m, s3, s8, s7, lq, lq, m, subn); /* M6 = (A21-A11)(B11+B12) -> s7 */
dsgemm(m, s4, s9, s2, lq, lq, m, subn); /* M7 = (A12-A22)(B21+B22) -> s2 */
comb7(m, C11, C12, C21, C22, ldc, sX, s0, s5, s1, s6, s7, s2, ldo);
return;
}
#undef PROD
}
static T *g_pool = 0;
static size_t g_pool_n = 0;
/* m211: the pool must hold 5 live M planes + TA + TB at EVERY non-leaf level
(FUSE_NL), or 5 at a leaf level; the levels below share the tail sequentially. */
static size_t pool_need(int n) {
size_t need = 0;
int m = n >> 1;
for (;;) {
if (m <= SCUT) { need += 7 * (size_t)m * m; break; }
need += 7 * ((size_t)m * m + 8);
m >>= 1;
}
return need + 1024;
}
static void strassen1(int n, const T *A, const T *B, T *C, int ldc) {
#if FUSE_NL
size_t need = pool_need(n);
#else
size_t need = 2 * (size_t)n * n + 1024; /* DA1: 5 live m^2 planes at a leaf level */
#endif
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; }
g_top = 1;
dsgemm(n, A, B, C, n, n, ldc, g_pool); g_top = 0;
}
#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, 1);
}
#endif
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 2.082 s | 428 MB + 48 KB | Accepted | Score: 100 | 显示更多 |