// lane m164 arm cpf_t1: combine, next-row prefetch of all 7 M streams, T1 (L2)
#define MK5O 1
#pragma GCC optimize("O3")
#include <immintrin.h>
#include <stddef.h>
#include <stdint.h>
static unsigned long long PHK[5], PHKC[5], _mk5t;
static unsigned long long PHN[5], PHNC[5];
static inline unsigned long long RDN(void){ }
static inline void MK5B(void){ }
static int MK5IDX(int h){ return h==2048?0:(h==1024?1:(h==512?2:(h==256?3:4))); }
static inline void MK5E(int hh){ }
static int MMNIDX(int n){ return n>=4096?0:(n>=2048?1:(n>=1024?2:3)); }
static inline void MMNA(int n, unsigned long long t0){ }
#include <stdlib.h>
#pragma GCC push_options
#pragma GCC target("avx2")
#ifndef ORD
#define ORD 1
#endif
typedef short TE;
#define PAD 32
/* lane m4k_h1024_streams: the padded copies of the four RAW operand blocks. */
static TE *g_cc[4];
static int g_ch = 0, g_cok = 0;
/* ARM KNOBS, and two measured traps avoided: (i) a PLAIN `static int` is
constant-propagated by gcc-9 and each arm gets its own .text (measured: 522419 /
523443 / 523571 B) -- the CSQ lane's documented trap; (ii) two separate `volatile int`s
still differ in .text, because gcc reorders file-scope statics per arm, so every
RIP-relative displacement to them moves (measured: 13 differing bytes each). A single
struct pins the layout and the arms then differ ONLY in .data. */
static TE *g_rc[4];
static int g_rh = 0;
struct Knobs { volatile int copy; volatile int use; volatile int alloc; volatile int wmode;
volatile int copy2; volatile int use2;
volatile int wA11; volatile int wA22; volatile int wB11; volatile int wB22;
volatile int lv; volatile int bp; volatile int abl; };
static struct Knobs g_k = { 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0 };
#define g_kbp (g_k.bp)
#define g_abl (g_k.abl)
#define g_copy (g_k.copy)
#define g_use (g_k.use)
#define g_alloc (g_k.alloc)
#define g_wmode (g_k.wmode)
#define g_copy2 (g_k.copy2)
#define g_use2 (g_k.use2)
#define g_wA11 (g_k.wA11)
#define g_wA22 (g_k.wA22)
#define g_wB11 (g_k.wB11)
#define g_wB22 (g_k.wB22)
#define g_lv (g_k.lv)
#ifndef CUT
#define CUT 256
#endif
#define MR 6
#define NR 16
#define KC 256
#define NC 512
#define UNROLL 4 /* in k-PAIRS */
#define TRS (((KC) + 15) & ~15) /* tr[] row stride: keep every row 32B-aligned */
// Every 256-bit address here is provably 32-byte aligned: every lead dimension is
// n+PAD (PAD=16 elements = 32 B) or h+PAD, every column offset is a multiple of 16
// elements, and every pool offset is a multiple of 64 elements. Using the *aligned*
// intrinsics stops GCC-9 (the judge's compiler) from splitting each access into
// 128-bit halves (vmovdqu xmm + vinserti128 / vextracti128 + vmovups xmm).
#define STV(p, v) _mm256_store_si256((__m256i *)(p), (v))
#define LDV(p) _mm256_load_si256((const __m256i *)(p))
#define LDU(p) _mm256_loadu_si256((const __m256i *)(p))
/* lane m4k_layout: Bpanel's OWN base moved by exactly 4096 B. Bpanel is the hottest
read structure on this row (196 608 B = 6 x 32 KB, read by every leaf's mikro) and
its base is the last address here that is ours. qbp8 moves it a FULL 32-KB cycle and
is the probe's MATCHED NULL: a page-phase effect must read it like the control. */
#define BPADMAX 128 /* FIXED allocation: g_bp's size is BPAD-independent */
#define BPAD 32 /* the swept knob: block-row stride, in ELEMENTS */
struct __attribute__((aligned(64))) { char _pad[4096]; TE _b[(size_t)NC * KC + (size_t)(NC / NR) * BPADMAX + 256]; char _x[16384]; } g_bp;
#define Bpanel g_bp._b
// K=2 panel entry: 32 bytes holding [a_k, a_{k+1}] repeated 8 times.
// vpbroadcastd of the 32-bit little-endian pair does exactly this in ONE uop.
static inline __m256i b32(const TE *p) {
__m256i r; __asm__("vpbroadcastd %1, %0" : "=x"(r) : "m"(*(const int *)p)); return r;
}
// B operand for a k-pair: [B[k][j], B[k+1][j], ...] interleaved (8 j's per 32 B).
static inline __m256i bilv(const TE *lo, const TE *hi) {
return _mm256_unpacklo_epi16(_mm256_castsi128_si256(_mm_loadu_si128((const __m128i *)lo)),
_mm256_castsi128_si256(_mm_loadu_si128((const __m128i *)hi)));
}
// --- k-PAIR mikro: A operand is a fused memory operand of vpmaddwd ---
// panel entry (row r, k-pair kk) = 32 B = [a_{2kk}, a_{2kk+1}] x 8 (vpbroadcastd)
// B operand (kk, j-half h) = 32 B = [B[k][j], B[k+1][j], ...] x 8 j's (unpck)
// vpmaddwd lane t = a_k*B[k][j_t] + a_{k+1}*B[k+1][j_t] -> 16 products / uop
// vpaddd accumulates the 8 i32 lanes; only the low 16 bits are ever read.
// COMPACT A operand: ONE vpbroadcastd (4-B load) per (row,k-pair) produces the whole
// [a_k,a_{k+1}]x8 lane pattern, and the register is used by BOTH the L and H muls.
#define KST(off, C0, C1) \
"vpbroadcastd " #off "(%[a]), %%ymm12\n\t" \
"vpmaddwd %%ymm12, %%ymm14, %%ymm13\n\t" \
"vpaddd %%ymm13, %" C0 ", %" C0 "\n\t" \
"vpmaddwd %%ymm12, %%ymm15, %%ymm13\n\t" \
"vpaddd %%ymm13, %" C1 ", %" C1 "\n\t"
#define KST5D(a5) \
"vpbroadcastd " #a5 "(%[a]), %%ymm12\n\t" \
"vpmaddwd %%ymm12, %%ymm14, %%ymm13\n\t" \
"vpmaddwd %%ymm12, %%ymm15, %%ymm12\n\t"
#define KSTF(boff, off, C0, C1) \
"vpbroadcastd " #off "(%[a]), %%ymm12\n\t" \
"vpmaddwd " #boff "(%[b]), %%ymm12, %%ymm13\n\t" \
"vpaddd %%ymm13, %" C0 ", %" C0 "\n\t" \
"vpmaddwd %%ymm12, %%ymm15, %%ymm13\n\t" \
"vpaddd %%ymm13, %" C1 ", %" C1 "\n\t"
#define KST5DF(boff, a5) \
"vpbroadcastd " #a5 "(%[a]), %%ymm12\n\t" \
"vpmaddwd " #boff "(%[b]), %%ymm12, %%ymm13\n\t" \
"vpmaddwd %%ymm12, %%ymm15, %%ymm12\n\t"
#define KGRP_P2(a0,a1,a2,a3,a4,a5,b0,b1) \
"vmovdqa " #b1 "(%[b]), %%ymm15\n\t" \
"vpaddd %%ymm13, %10, %10\n\t" \
"vpaddd %%ymm12, %11, %11\n\t" \
KSTF(b0,a0,"0","1") KSTF(b0,a1,"2","3") KSTF(b0,a2,"4","5") \
KSTF(b0,a3,"6","7") KSTF(b0,a4,"8","9") KST5DF(b0,a5)
/* ---- PROGRAM-ORDER PIPELINE OF THE A BROADCAST (lane mmml_lane, mmms4k) -------------
The measured signature on this row is a ~1.85-cycle-per-iteration STALL on a body that
is otherwise front-end bound (32.75 cycles theoretical, 34.6 measured), and the mikro
is 81.9 % of the row. The shortest consumption distance in KGRP_P2 belongs to the A
operand: `vpbroadcastd a(%[a]), %ymm12` was consumed by the very next instruction, six
times per group, and that broadcast is a load whose latency is 5-7 cycles.
THIS IS A PURE REORDER. ymm14 is dead in KGRP_P2 (the L half of B is a folded memory
operand of vpmaddwd), so the six broadcasts alternate between ymm12 and ymm14 and each
one is issued ~5 instructions ahead of its own multiplies instead of one behind. Same
six broadcasts, same thirty multiply/accumulate instructions, same operands, same
addresses: only the ISSUE ORDER changes, so it is byte-exact by construction.
ymm12 still carries the deferred row-5 ring, and that is why the odd rows keep it: the
a0 broadcast is hoisted into ymm14 so the group-top ring drain still reads ymm12. */
#define KST_A(boff, off, C0, C1, AR) \
"vpmaddwd " #boff "(%[b]), " AR ", %%ymm13\n\t" \
"vpmaddwd " AR ", %%ymm15, " AR "\n\t" \
"vpaddd %%ymm13, %" C0 ", %" C0 "\n\t" \
"vpaddd " AR ", %" C1 ", %" C1 "\n\t"
#define KST5D_A(boff, AR) \
"vpmaddwd " #boff "(%[b]), " AR ", %%ymm13\n\t" \
"vpmaddwd " AR ", %%ymm15, %%ymm12\n\t"
#define KGRP_P2A(a0,a1,a2,a3,a4,a5,b0,b1) \
"vpbroadcastd " #a0 "(%[a]), %%ymm14\n\t" \
"vmovdqa " #b1 "(%[b]), %%ymm15\n\t" \
"vpaddd %%ymm13, %10, %10\n\t" \
"vpaddd %%ymm12, %11, %11\n\t" \
"vpbroadcastd " #a1 "(%[a]), %%ymm12\n\t" \
KST_A(b0,a0,"0","1","%%ymm14") \
"vpbroadcastd " #a2 "(%[a]), %%ymm14\n\t" \
KST_A(b0,a1,"2","3","%%ymm12") \
"vpbroadcastd " #a3 "(%[a]), %%ymm12\n\t" \
KST_A(b0,a2,"4","5","%%ymm14") \
"vpbroadcastd " #a4 "(%[a]), %%ymm14\n\t" \
KST_A(b0,a3,"6","7","%%ymm12") \
"vpbroadcastd " #a5 "(%[a]), %%ymm12\n\t" \
KST_A(b0,a4,"8","9","%%ymm14") \
KST5D_A(b0,"%%ymm12")
#define KFULL_A_DR PF(1024) KGRP_P2A(0,2112,4224,6336,8448,10560, 0,32) \
PF(1088) KGRP_P2A(4,2116,4228,6340,8452,10564, 64,96) \
PF(1152) KGRP_P2A(8,2120,4232,6344,8456,10568, 128,160) \
PF(1216) KGRP_P2A(12,2124,4236,6348,8460,10572, 192,224)
#define KFULL_A PF(1024) KGRP_P2A(0,512,1024,1536,2048,2560, 0,32) \
PF(1088) KGRP_P2A(4,516,1028,1540,2052,2564, 64,96) \
PF(1152) KGRP_P2A(8,520,1032,1544,2056,2568, 128,160) \
PF(1216) KGRP_P2A(12,524,1036,1548,2060,2572, 192,224)
/* lane MS1 arm a4_ordDR -- a4_ord AND the same issue-order edit on the SECOND mikro body
(micro_dr / KFULL_A_DR, taken when lda == 1056), whose A offsets are 2112-strided.
SAME instructions, SAME operands, SAME accumulator updates as KFULL_A_DR -- only the
ISSUE ORDER of the six vpbroadcastd changes. Byte-exact by construction. */
#define KO_BC(o) "vpbroadcastd " #o "(%[a]), %%ymm12\n\t"
#define KO_MAD(bo) "vpmaddwd " #bo "(%[b]), %%ymm12, %%ymm13\n\t" "vpmaddwd %%ymm12, %%ymm15, %%ymm14\n\t"
#define KO_ADD(C0,C1) "vpaddd %%ymm13, %" C0 ", %" C0 "\n\t" "vpaddd %%ymm14, %" C1 ", %" C1 "\n\t"
#define KO_MAD5D(bo) "vpmaddwd " #bo "(%[b]), %%ymm12, %%ymm13\n\t" "vpmaddwd %%ymm12, %%ymm15, %%ymm12\n\t"
#define KGRP_P2O(a0,a1,a2,a3,a4,a5,b0,b1) \
"vmovdqa " #b1 "(%[b]), %%ymm15\n\t" \
"vpaddd %%ymm13, %10, %10\n\t" \
"vpaddd %%ymm12, %11, %11\n\t" \
KO_BC(a0) KO_MAD(b0) KO_BC(a1) KO_ADD("0","1") KO_MAD(b0) KO_BC(a2) KO_ADD("2","3") \
KO_MAD(b0) KO_BC(a3) KO_ADD("4","5") KO_MAD(b0) KO_BC(a4) KO_ADD("6","7") \
KO_MAD(b0) KO_BC(a5) KO_ADD("8","9") KO_MAD5D(b0)
#define KFULL_O PF(1024) KGRP_P2O(0,512,1024,1536,2048,2560, 0,32) \
PF(1088) KGRP_P2O(4,516,1028,1540,2052,2564, 64,96) \
PF(1152) KGRP_P2O(8,520,1032,1544,2056,2568, 128,160) \
PF(1216) KGRP_P2O(12,524,1036,1548,2060,2572, 192,224)
#define KFULL_O_DR PF(1024) KGRP_P2O(0,2112,4224,6336,8448,10560, 0,32) \
PF(1088) KGRP_P2O(4,2116,4228,6340,8452,10564, 64,96) \
PF(1152) KGRP_P2O(8,2120,4232,6344,8456,10568, 128,160) \
PF(1216) KGRP_P2O(12,2124,4236,6348,8460,10572, 192,224)
#define KGRP(a0,a1,a2,a3,a4,a5,b0,b1) \
"vmovdqa " #b0 "(%[b]), %%ymm14\n\t" \
"vmovdqa " #b1 "(%[b]), %%ymm15\n\t" \
KST(a0,"0","1") KST(a1,"2","3") KST(a2,"4","5") KST(a3,"6","7") KST(a4,"8","9") KST(a5,"10","11")
#define TRB (TRS * 2) /* bytes per tr row */
#if TRS != 256
#error "compact A operand hard-codes the tr row stride in the asm immediates"
#endif
#define KGRP1 KGRP(0,512,1024,1536,2048,2560, 0,32)
#define PF(o) "prefetcht0 " #o "(%[b])\n\t"
#define KFULL PF(512) KGRP_P2(0,512,1024,1536,2048,2560, 0,32) \
PF(576) KGRP_P2(4,516,1028,1540,2052,2564, 64,96) \
PF(640) KGRP_P2(8,520,1032,1544,2056,2568, 128,160) \
PF(704) KGRP_P2(12,524,1036,1548,2060,2572, 192,224)
#define MICRO_CLOB "ymm12", "ymm13", "ymm14", "ymm15", "cc", "memory"
// 8 i32 -> 8 i16 packed into the low 16 bytes
#define FOLD16(v) _mm256_permute4x64_epi64(_mm256_shuffle_epi8((v), FOLDM), 0x08)
static const __m256i FOLDM = { (long long)0x0D0C090805040100LL, (long long)0x8080808080808080LL,
(long long)0x0D0C090805040100LL, (long long)0x8080808080808080LL };
// L/H panel lane order: b0 = L = [j0..j3 | j8..j11], b1 = H = [j4..j7 | j12..j15]
// (lane index t of the acc = column perm(t)); the store must re-interleave the two
// halves, which is exactly one vpunpcklqdq of the two folded vectors.
#define FH(v) _mm256_shuffle_epi8((v), FOLDM)
#define ST_ROW(r, v0, v1) \
if (rows > (r)) { \
TE *p0 = C + (size_t)(r) * ldc; \
__m256i zz = _mm256_unpacklo_epi64(FH(v0), FH(v1)); \
if (OP == 0) { STV(p0, zz); } \
else { STV(p0, _mm256_add_epi16(LDV(p0), zz)); } \
}
#define MICRO_DEF(NAME, BODY, ADV_A, ADV_B) \
template <int OP> \
static inline void NAME(long k, const TE *ap, const TE *bp, TE *C, int ldc, int rows) { \
__m256i c00 = _mm256_setzero_si256(), c01 = c00, c10 = c00, c11 = c00, c20 = c00, c21 = c00; \
__m256i c30 = c00, c31 = c00, c40 = c00, c41 = c00, c50 = c00, c51 = c00; \
__asm__ volatile( \
"vpxor %%ymm12, %%ymm12, %%ymm12\n\t" \
"vpxor %%ymm13, %%ymm13, %%ymm13\n\t" \
"1:\n\t" BODY \
"addq $" #ADV_A ", %[a]\n\t" \
"addq $" #ADV_B ", %[b]\n\t" \
"decq %[k]\n\t" \
"jnz 1b\n\t" \
"vpaddd %%ymm13, %10, %10\n\t" \
"vpaddd %%ymm12, %11, %11\n\t" \
: "+x"(c00), "+x"(c01), "+x"(c10), "+x"(c11), "+x"(c20), "+x"(c21), \
"+x"(c30), "+x"(c31), "+x"(c40), "+x"(c41), "+x"(c50), "+x"(c51), \
[a] "+r"(ap), [b] "+r"(bp), [k] "+r"(k) \
: \
: MICRO_CLOB); \
ST_ROW(0, c00, c01) ST_ROW(1, c10, c11) ST_ROW(2, c20, c21) \
ST_ROW(3, c30, c31) ST_ROW(4, c40, c41) ST_ROW(5, c50, c51) \
}
MICRO_DEF(micro, KFULL_O, 16, 256) /* four k-pairs per iteration */
MICRO_DEF(micro_dr, KFULL_O_DR, 16, 256) /* four k-pairs per iteration */
MICRO_DEF(micro1, KGRP1, 4, 64) /* one k-pair */
// tail helper (nc not a multiple of 2*NR): same L/H lane order as the fast path.
// entry lo = [ j0..j3 | j8..j11 ] entry hi = [ j4..j7 | j12..j15 ]
static inline __m256i bp2x(const TE *p0, const TE *p1, const TE *q0, const TE *q1, int op, int hi) {
__m128i a0 = _mm_loadu_si128((const __m128i *)p0);
__m128i b0 = _mm_loadu_si128((const __m128i *)p1);
__m128i a1 = _mm_loadu_si128((const __m128i *)(p0 + 8));
__m128i b1 = _mm_loadu_si128((const __m128i *)(p1 + 8));
if (op == 1) { a0 = _mm_add_epi16(a0, _mm_loadu_si128((const __m128i *)q0));
b0 = _mm_add_epi16(b0, _mm_loadu_si128((const __m128i *)q1));
a1 = _mm_add_epi16(a1, _mm_loadu_si128((const __m128i *)(q0 + 8)));
b1 = _mm_add_epi16(b1, _mm_loadu_si128((const __m128i *)(q1 + 8))); }
else if (op == 2) { a0 = _mm_sub_epi16(a0, _mm_loadu_si128((const __m128i *)q0));
b0 = _mm_sub_epi16(b0, _mm_loadu_si128((const __m128i *)q1));
a1 = _mm_sub_epi16(a1, _mm_loadu_si128((const __m128i *)(q0 + 8)));
b1 = _mm_sub_epi16(b1, _mm_loadu_si128((const __m128i *)(q1 + 8))); }
if (!hi) return _mm256_inserti128_si256(_mm256_castsi128_si256(_mm_unpacklo_epi16(a0, b0)),
_mm_unpacklo_epi16(a1, b1), 1);
return _mm256_inserti128_si256(_mm256_castsi128_si256(_mm_unpackhi_epi16(a0, b0)),
_mm_unpackhi_epi16(a1, b1), 1);
}
#define OP_RAW 0
#define OP_ADD 1
#define OP_SUB 2
// =====================================================================================
// THE FUSED tr[] BUILD -- lane m4k_4cut
// operand = (P1 aop P2) aop2 (P3 aop P4)
// OP1 = aop (within-pair op for BOTH pairs) ; OP2 = aop2 (OP_RAW => pair 2 absent)
// P1,P2 use stride `lda`; P3,P4 use stride `lda2`.
// With OP1 == OP_RAW this reduces to the pre-existing single-term copy; with OP2 ==
// OP_RAW it is the pre-existing 2-term build, instruction-for-instruction.
// NEVER a second pass over tr[]: the join happens in registers.
// =====================================================================================
// generic scalar tail (only reached when kc is not a multiple of 32 -- never at
// n = 256 with KC = 256, kept for correctness in general)
static void trb_tail(TE tr[MR][TRS], int rows, int kc, int k0,
const TE *w1, const TE *w2, const TE *w3, const TE *w4, int lda, int lda2,
int op1, int op2) {
for (int k = k0; k < kc; k++) {
for (int r = 0; r < rows; r++) {
int v = (int)w1[(size_t)r * lda + k];
if (op1 != OP_RAW) {
int y = (int)w2[(size_t)r * lda2 + k];
v += (op1 == OP_ADD) ? y : -y;
}
if (op2 != OP_RAW) {
int y = (int)w3[(size_t)r * lda2 + k];
if (op1 != OP_RAW) {
int z = (int)w4[(size_t)r * lda2 + k];
y += (op1 == OP_ADD) ? z : -z;
}
v += (op2 == OP_ADD) ? y : -y;
}
tr[r][k] = (TE)v;
}
}
}
template <int OP1, int OP2>
static inline void trb(TE tr[MR][TRS], int rows, int kc,
const TE *w1, const TE *w2, const TE *w3, const TE *w4,
int lda, int lda2) {
int k = 0;
for (; k + 32 <= kc; k += 32, w1 += 32, w2 += 32, w3 += 32, w4 += 32)
for (int r = 0; r < rows; r++) {
const TE *s1 = w1 + (size_t)r * lda;
__m256i x0 = LDV(s1), x1 = LDV(s1 + 16);
if (OP1 != OP_RAW) {
const TE *s2 = w2 + (size_t)r * lda2;
if (OP1 == OP_ADD) { x0 = _mm256_add_epi16(x0, LDV(s2)); x1 = _mm256_add_epi16(x1, LDV(s2 + 16)); }
else { x0 = _mm256_sub_epi16(x0, LDV(s2)); x1 = _mm256_sub_epi16(x1, LDV(s2 + 16)); }
}
if (OP2 == OP_RAW) {
STV(tr[r] + k, x0); STV(tr[r] + k + 16, x1);
} else {
const TE *s3 = w3 + (size_t)r * lda2;
__m256i y0 = LDV(s3), y1 = LDV(s3 + 16);
if (OP1 != OP_RAW) {
const TE *s4 = w4 + (size_t)r * lda2;
if (OP1 == OP_ADD) { y0 = _mm256_add_epi16(y0, LDV(s4)); y1 = _mm256_add_epi16(y1, LDV(s4 + 16)); }
else { y0 = _mm256_sub_epi16(y0, LDV(s4)); y1 = _mm256_sub_epi16(y1, LDV(s4 + 16)); }
}
if (OP2 == OP_ADD) { STV(tr[r] + k, _mm256_add_epi16(x0, y0)); STV(tr[r] + k + 16, _mm256_add_epi16(x1, y1)); }
else { STV(tr[r] + k, _mm256_sub_epi16(x0, y0)); STV(tr[r] + k + 16, _mm256_sub_epi16(x1, y1)); }
}
}
if (k < kc) trb_tail(tr, rows, kc, k, w1, w2, w3, w4, lda, lda2, OP1, OP2);
}
// =====================================================================================
// THE FUSED B-PANEL PACK -- lane m4k_4cut
// packed value = ((B[x] bop B2[x]) bop2 (B3[x] bop B4[x])), x = row k and row k+1
// B, B2 stride `ldb`; B3, B4 stride `ldb2`. OP2 == OP_RAW => pair 2 absent.
// Same lane order, same addresses, same loop nests as the shipped packer.
// =====================================================================================
static void packb_tail(const TE *B, int ldb, const TE *B2, const TE *B3, const TE *B4, int ldb2,
int kc, int nc, int pc, int jc, int jr, int bop, int bop2) {
const int kn = kc >> 1;
for (; jr < nc; jr += NR) {
TE *q = Bpanel + (size_t)(jr / NR) * (kc * NR + BPAD);
for (int kp = 0; kp < kn; kp++) {
const int k = kp << 1;
/* mirror the 256-bit fast path element for element:
vec 0 -> q+0 (row k, cols jr..jr+7 interleaved with row k+1)
vec 1 -> q+16 (row k, cols jr+8..jr+15)
vec 2 -> q2+0 (row k, cols jr+16..jr+23)
vec 3 -> q2+16(row k, cols jr+24..jr+31) */
for (int vec = 0; vec < 4; vec++) {
TE *dst = q + (vec & 1) * 16 + (size_t)(vec >> 1) * (kc * NR + BPAD);
const int colbase = jc + jr + (vec >> 1) * 16 + ((vec & 1) ? 8 : 0);
for (int t = 0; t < 16; t++) {
const int col = colbase + (t >> 1);
const int rw = pc + k + (t & 1);
int a = (int)B[(size_t)rw * ldb + col];
if (bop != OP_RAW) {
int x = (int)B2[(size_t)rw * ldb2 + col];
a += (bop == OP_ADD) ? x : -x;
}
if (bop2 != OP_RAW) {
int x = (int)B3[(size_t)rw * ldb2 + col];
if (bop != OP_RAW) {
int x2 = (int)B4[(size_t)rw * ldb2 + col];
x += (bop == OP_ADD) ? x2 : -x2;
}
a += (bop2 == OP_ADD) ? x : -x;
}
dst[t] = (TE)a;
}
}
q += 2 * NR;
}
}
}
template <int OP1, int OP2>
static void packb(const TE *B, int ldb, const TE *B2, const TE *B3, const TE *B4, int ldb2,
int kc, int nc, int pc, int jc) {
const int kn = kc >> 1;
int jr = 0;
const TE *b2 = (OP1 == OP_RAW) ? B : B2;
const TE *b3 = (OP2 == OP_RAW) ? B : B3;
const TE *b4 = (OP1 == OP_RAW) ? B3 : B4;
for (; jr + 2 * NR <= nc; jr += 2 * NR) {
TE *q = Bpanel + (size_t)(jr / NR) * (kc * NR + BPAD);
TE *q2 = q + (size_t)(kc * NR + BPAD);
for (int kp = 0; kp < kn; kp++) {
const int k = kp << 1;
const TE *s0 = B + (size_t)(pc + k) * ldb + jc + jr;
const TE *s1 = B + (size_t)(pc + k + 1) * ldb + jc + jr;
__m256i A0 = LDV(s0), A1 = LDV(s0 + 16), B0 = LDV(s1), B1 = LDV(s1 + 16);
if (OP1 != OP_RAW) {
const TE *t0 = b2 + (size_t)(pc + k) * ldb2 + jc + jr;
const TE *t1 = b2 + (size_t)(pc + k + 1) * ldb2 + jc + jr;
if (OP1 == OP_ADD) { A0 = _mm256_add_epi16(A0, LDV(t0)); A1 = _mm256_add_epi16(A1, LDV(t0 + 16));
B0 = _mm256_add_epi16(B0, LDV(t1)); B1 = _mm256_add_epi16(B1, LDV(t1 + 16)); }
else { A0 = _mm256_sub_epi16(A0, LDV(t0)); A1 = _mm256_sub_epi16(A1, LDV(t0 + 16));
B0 = _mm256_sub_epi16(B0, LDV(t1)); B1 = _mm256_sub_epi16(B1, LDV(t1 + 16)); }
}
if (OP2 != OP_RAW) {
const TE *u0 = b3 + (size_t)(pc + k) * ldb2 + jc + jr;
const TE *u1 = b3 + (size_t)(pc + k + 1) * ldb2 + jc + jr;
__m256i C0 = LDV(u0), C1 = LDV(u0 + 16), D0 = LDV(u1), D1 = LDV(u1 + 16);
if (OP1 != OP_RAW) {
const TE *v0 = b4 + (size_t)(pc + k) * ldb2 + jc + jr;
const TE *v1 = b4 + (size_t)(pc + k + 1) * ldb2 + jc + jr;
if (OP1 == OP_ADD) { C0 = _mm256_add_epi16(C0, LDV(v0)); C1 = _mm256_add_epi16(C1, LDV(v0 + 16));
D0 = _mm256_add_epi16(D0, LDV(v1)); D1 = _mm256_add_epi16(D1, LDV(v1 + 16)); }
else { C0 = _mm256_sub_epi16(C0, LDV(v0)); C1 = _mm256_sub_epi16(C1, LDV(v0 + 16));
D0 = _mm256_sub_epi16(D0, LDV(v1)); D1 = _mm256_sub_epi16(D1, LDV(v1 + 16)); }
}
if (OP2 == OP_ADD) { A0 = _mm256_add_epi16(A0, C0); A1 = _mm256_add_epi16(A1, C1);
B0 = _mm256_add_epi16(B0, D0); B1 = _mm256_add_epi16(B1, D1); }
else { A0 = _mm256_sub_epi16(A0, C0); A1 = _mm256_sub_epi16(A1, C1);
B0 = _mm256_sub_epi16(B0, D0); B1 = _mm256_sub_epi16(B1, D1); }
}
STV(q, _mm256_unpacklo_epi16(A0, B0));
STV(q + 16, _mm256_unpackhi_epi16(A0, B0));
STV(q2, _mm256_unpacklo_epi16(A1, B1));
STV(q2 + 16, _mm256_unpackhi_epi16(A1, B1));
q += 2 * NR; q2 += 2 * NR;
}
}
if (jr < nc) packb_tail(B, ldb, B2, B3, B4, ldb2, kc, nc, pc, jc, jr, OP1, OP2);
}
// =====================================================================================
// leaf(): C = (A op A2) * (B op B2) -- with the m4k_4cut second operand pair.
// A-side operand = (A aop A2) aop2 (A3 aop A4) [aop2 == OP_RAW => pair 2 absent]
// B-side operand = (B bop B2) bop2 (B3 bop B4)
// Both aop2 and bop2 apply the leaf's OWN within-pair op to the second pair; that is
// exactly the algebra the n == 4*CUT fusion produces (see the header).
// =====================================================================================
#pragma GCC push_options
#pragma GCC target("avx2")
/* ---------------------------------------------------------------------------------------
lane m4k_fe_b -- THE RAGGED-TAIL MIKRO.
`leaf` steps its panel loop by MR = 6, and 256 % 6 == 4, so the LAST ic block of every
leaf(256) has rows = 4 while the mikro body unconditionally computes 6 rows: 2 of every
258 row-slots, 0.78 % of every mikro call, are wasted. This adds a FOUR-ROW mikro that
is used for that block only, so the waste is not paid -- without changing MR (changing
MR to 4 costs more in operand reuse: 2 B loads per 4 rows instead of per 6 rows, measured
+1.752 % on the judge as v/mr4b.cpp).
The four-row body is the row's own mechanism with two fewer rows:
2 vmovdqa + 4 vpbroadcastd + 8 vpmaddwd + 8 vpaddd = 22 insns per KGRP (128 MACs),
the same 2 ALU uops per 16 16-bit MACs and the same |S| = 2 deferred ring.
------------------------------------------------------------------------------------- */
#define KGRP4_B(a0,a1,a2,a3,b0,b1) \
"vmovdqa " #b0 "(%[b]), %%ymm14\n\t" \
"vmovdqa " #b1 "(%[b]), %%ymm15\n\t" \
"vpaddd %%ymm13, %6, %6\n\t" \
"vpaddd %%ymm12, %7, %7\n\t" \
KST(a0,"0","1") KST(a1,"2","3") KST(a2,"4","5") KST5D(a3)
#define KGRP4 KGRP4_B(0,512,1024,1536, 0,32) \
KGRP4_B(4,516,1028,1540, 64,96) \
KGRP4_B(8,520,1032,1544, 128,160) \
KGRP4_B(12,524,1036,1548, 192,224)
template <int OP>
static inline void micro4(long k, const TE *ap, const TE *bp, TE *C, int ldc, int rows) {
__m256i c00 = _mm256_setzero_si256(), c01 = c00, c10 = c00, c11 = c00;
__m256i c20 = c00, c21 = c00, c30 = c00, c31 = c00;
__asm__ volatile(
"vpxor %%ymm12, %%ymm12, %%ymm12\n\t"
"vpxor %%ymm13, %%ymm13, %%ymm13\n\t"
"1:\n\t" KGRP4
"addq $16, %[a]\n\t"
"addq $256, %[b]\n\t"
"decq %[k]\n\t"
"jnz 1b\n\t"
"vpaddd %%ymm13, %6, %6\n\t"
"vpaddd %%ymm12, %7, %7\n\t"
: "+x"(c00), "+x"(c01), "+x"(c10), "+x"(c11), "+x"(c20), "+x"(c21),
"+x"(c30), "+x"(c31), [a] "+r"(ap), [b] "+r"(bp), [k] "+r"(k)
:
: MICRO_CLOB);
ST_ROW(0, c00, c01) ST_ROW(1, c10, c11) ST_ROW(2, c20, c21) ST_ROW(3, c30, c31)
}
#pragma GCC pop_options
static void leaf(TE *C, int ldc, const TE *A, int lda, const TE *A2, int lda2,
const TE *A3, const TE *A4, int aop, int aop2,
const TE *B, int ldb, const TE *B2, int ldb2,
const TE *B3, const TE *B4, int bop, int bop2, int n) {
if (g_abl) { lda = lda2 = ldb = ldb2 = g_abl; }
const int nbc = n - (n % NR);
for (int jc = 0; jc < nbc; jc += NC) {
const int nc = (nbc - jc < NC) ? (nbc - jc) : NC;
for (int pc = 0; pc < n; pc += KC) {
const int kc = (n - pc < KC) ? (n - pc) : KC;
{
if (bop2 == OP_RAW) {
if (bop == OP_RAW) packb<OP_RAW, OP_RAW>(B, ldb, 0, 0, 0, ldb2, kc, nc, pc, jc);
else if (bop == OP_ADD) packb<OP_ADD, OP_RAW>(B, ldb, B2, 0, 0, ldb2, kc, nc, pc, jc);
else packb<OP_SUB, OP_RAW>(B, ldb, B2, 0, 0, ldb2, kc, nc, pc, jc);
} else if (bop == OP_RAW) {
if (bop2 == OP_ADD) packb<OP_RAW, OP_ADD>(B, ldb, 0, B3, B3, ldb2, kc, nc, pc, jc);
else packb<OP_RAW, OP_SUB>(B, ldb, 0, B3, B3, ldb2, kc, nc, pc, jc);
} else if (bop == OP_ADD) {
if (bop2 == OP_ADD) packb<OP_ADD, OP_ADD>(B, ldb, B2, B3, B4, ldb2, kc, nc, pc, jc);
else packb<OP_ADD, OP_SUB>(B, ldb, B2, B3, B4, ldb2, kc, nc, pc, jc);
} else {
if (bop2 == OP_ADD) packb<OP_SUB, OP_ADD>(B, ldb, B2, B3, B4, ldb2, kc, nc, pc, jc);
else packb<OP_SUB, OP_SUB>(B, ldb, B2, B3, B4, ldb2, kc, nc, pc, jc);
}
}
for (int ic = 0; ic < n; ic += MR) {
const int rows = (n - ic < MR) ? (n - ic) : MR;
// COMPACT A operand: tr[r][k] = ((A aop A2) aop2 (A3 aop A4))[ic+r][pc+k],
// row stride TRS. No 32-B replicated panel: the mikro broadcasts each 4-B
// k-pair itself.
struct __attribute__((aligned(64))) { char pad[160u]; TE t[MR][TRS]; char tail[864u]; } g_trb;
if (n < 0) { g_trb.pad[0] = 0; g_trb.tail[0] = 0; }
#define tr (g_trb.t)
{ /* lane m4k_re: k-major, 32-element chunks, hoisted row pointers.
Board-verified -0.44 % on the compact KC=128 design (677.504 -> 674.509,
sid 91304); composed onto KC=256 here unchanged. */
const TE *w1 = A + (size_t)ic * lda + pc;
const TE *w2 = A2 ? A2 + (size_t)ic * lda2 + pc : 0;
const TE *w3 = A3 ? A3 + (size_t)ic * lda2 + pc : 0;
const TE *w4 = A4 ? A4 + (size_t)ic * lda2 + pc : 0;
if (aop2 == OP_RAW) {
if (aop == OP_RAW && aop2 == OP_RAW && lda == 1056 && rows == MR) {
/* trb would be a byte-for-byte copy of A rows at a 2112 B stride: read A directly.
GUARD: only when the leaf row stride is exactly 1056 elements. */
} else if (aop == OP_RAW) trb<OP_RAW, OP_RAW>(tr, rows, kc, w1, w1, w1, w1, lda, lda2);
else if (aop == OP_ADD) trb<OP_ADD, OP_RAW>(tr, rows, kc, w1, w2, w1, w1, lda, lda2);
else trb<OP_SUB, OP_RAW>(tr, rows, kc, w1, w2, w1, w1, lda, lda2);
} else if (aop == OP_RAW) {
if (aop2 == OP_ADD) trb<OP_RAW, OP_ADD>(tr, rows, kc, w1, w1, w3, w3, lda, lda2);
else trb<OP_RAW, OP_SUB>(tr, rows, kc, w1, w1, w3, w3, lda, lda2);
} else if (aop == OP_ADD) {
if (aop2 == OP_ADD) trb<OP_ADD, OP_ADD>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
else trb<OP_ADD, OP_SUB>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
} else {
if (aop2 == OP_ADD) trb<OP_SUB, OP_ADD>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
else trb<OP_SUB, OP_SUB>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
}
}
for (int r = rows; r < MR; r++) // deterministic, never stored
for (int k = 0; k < kc; k++) tr[r][k] = 0;
const TE *trp = (aop == OP_RAW && aop2 == OP_RAW && lda == 1056 && rows == MR) ? (A + (size_t)ic * lda + pc) : (const TE *)tr;
const TE *b0 = Bpanel;
TE *cp = C + (size_t)ic * ldc + jc;
for (int jr = 0; jr < nc; jr += NR) {
const int kn = kc >> 1; /* k-pairs */
const int kf = kn & ~(UNROLL - 1);
if (rows == MR && aop == OP_RAW && aop2 == OP_RAW && lda == 1056) micro_dr<0>(kf / UNROLL, trp, b0, cp + jr, ldc, rows);
else if (rows == MR) micro<0>(kf / UNROLL, trp, b0, cp + jr, ldc, rows);
else micro4<0>(kf / UNROLL, trp, b0, cp + jr, ldc, rows);
b0 += (size_t)(kc * NR + BPAD);
}
}
}
}
if (nbc < n) {
for (int i = 0; i < n; i++)
for (int j = nbc; j < n; j++) {
unsigned acc = 0;
for (int k = 0; k < n; k++)
acc += (unsigned)(unsigned short)A[(size_t)i * lda + k] * (unsigned)(unsigned short)B[(size_t)k * ldb + j];
C[(size_t)i * ldc + j] = (TE)acc;
}
}
}
// ---------------- Strassen driver ----------------
static TE *g_pool;
static size_t g_off;
/* lane m4k_page2: the DEEP arena. Below the top node every pool region is dead
long before the top-level combine writes C, and C itself is dead until then, so
the deep regions are allocated INSIDE C's own buffer (g_scr) instead of the pool.
g_topn is the top node's n; every node with n < g_topn draws from g_scr, and the
top node -- and only it -- draws from g_pool. g_topn = -1 disables the split
(used when the deep stack does not fit in C). */
static TE *g_scr;
static size_t g_soff;
static int g_topn = -1;
// form the 5 A-sums and 5 B-sums of an h-split block
static void mk5(TE *AS, TE *BS, const TE *A, int lda, const TE *B, int ldb, int h,
int bipA, int bipB) {
const size_t q = (size_t)(h + PAD) * h;
const size_t qs = q + 2048;
/* LANE m4k_copyres: THE BLOCK-PACKED SUM LAYOUT. `bpw` (h == 2048 = the row's top
node) writes each sum as FOUR (h/2) x (h/2 + PAD) blocks instead of one h x (h+PAD)
array, i.e. in the layout the n = h consumer reads its quadrants in. qsx = the slot
stride = 4 blocks; hph = the block row stride = the CONSUMER's h + PAD.
`bipA`/`bipB` say the operand handed IN is already in that layout (so its own four
quadrants live at A + k*qs and their row stride is this node's h + PAD). */
const size_t qb = (size_t)((h >> 1) + PAD) * (h >> 1) + 2048;
const int bpw = (g_kbp >= 1 && h == 2048) || (g_kbp >= 2 && h == 1024);
const int hb = (h == 2048) ? 10 : 9; /* log2(h/2) */
const int hm = (h == 2048) ? 1023 : 511;
/* qsx is the SLOT stride and does NOT depend on the knob: both arms allocate and
address the same bytes, so the only difference is the byte ORDER inside a slot. */
const size_t qsx = (h == 2048 || h == 1024) ? 4 * qb : qs;
const int hph = (h >> 1) + PAD;
const TE *a11, *a12, *a21, *a22, *b11, *b12, *b21, *b22;
if (bipA) { a11 = A; a12 = A + qs; a21 = A + 2 * qs; a22 = A + 3 * qs; lda = h + PAD; }
else { a11 = A; a12 = A + h; a21 = A + (size_t)h * lda; a22 = a21 + h; }
if (bipB) { b11 = B; b12 = B + qs; b21 = B + 2 * qs; b22 = B + 3 * qs; ldb = h + PAD; }
else { b11 = B; b12 = B + h; b21 = B + (size_t)h * ldb; b22 = b21 + h; }
TE *const *const cd = (h == g_ch) ? g_cc : g_rc;
const int carm = (g_copy && g_cok && h == g_ch && g_cc[0]) ||
(g_copy2 && g_cok && g_rh && h == g_rh && g_rc[0]);
/* a block-packed side needs NO copy: its quadrant IS already a contiguous block. */
const int copyA = carm && !bipA, copyB = carm && !bipB;
if (copyA || copyB) {
for (int i = 0; i < h; i++) {
const TE *r11 = a11 + (size_t)i * lda, *r12 = a12 + (size_t)i * lda;
const TE *r21 = a21 + (size_t)i * lda, *r22 = a22 + (size_t)i * lda;
TE *o0 = bpw ? (AS + ((size_t)((i >> hb) << 1)) * qb + (size_t)(i & hm) * (size_t)hph)
: (AS + (size_t)i * (h + PAD));
TE *o1 = o0 + qsx, *o2 = o1 + qsx, *o3 = o2 + qsx, *o4 = o3 + qsx;
TE *c11r = cd[0] + (size_t)i * (h + PAD), *c22r = cd[1] + (size_t)i * (h + PAD);
int j = 0;
#define MKAC(J) { const int jj = (J); \
__m256i q11 = LDV(r11 + jj), q12 = LDV(r12 + jj), q21 = LDV(r21 + jj), q22 = LDV(r22 + jj); \
__m256i w11 = LDV(r11 + jj + 16), w12 = LDV(r12 + jj + 16), w21 = LDV(r21 + jj + 16), w22 = LDV(r22 + jj + 16); \
_mm256_stream_si256((__m256i *)(o0 + jj), _mm256_add_epi16(q11, q22)); \
_mm256_stream_si256((__m256i *)(o0 + jj + 16), _mm256_add_epi16(w11, w22)); \
_mm256_stream_si256((__m256i *)(o1 + jj), _mm256_add_epi16(q21, q22)); \
_mm256_stream_si256((__m256i *)(o1 + jj + 16), _mm256_add_epi16(w21, w22)); \
_mm256_stream_si256((__m256i *)(o2 + jj), _mm256_add_epi16(q11, q12)); \
_mm256_stream_si256((__m256i *)(o2 + jj + 16), _mm256_add_epi16(w11, w12)); \
_mm256_stream_si256((__m256i *)(o3 + jj), _mm256_sub_epi16(q21, q11)); \
_mm256_stream_si256((__m256i *)(o3 + jj + 16), _mm256_sub_epi16(w21, w11)); \
_mm256_stream_si256((__m256i *)(o4 + jj), _mm256_sub_epi16(q12, q22)); \
_mm256_stream_si256((__m256i *)(o4 + jj + 16), _mm256_sub_epi16(w12, w22)); \
_mm256_stream_si256((__m256i *)(c11r + jj), q11); \
_mm256_stream_si256((__m256i *)(c11r + jj + 16), w11); \
_mm256_stream_si256((__m256i *)(c22r + jj), q22); \
_mm256_stream_si256((__m256i *)(c22r + jj + 16), w22); }
#define MKACS(J) { const int jj = (J); \
__m256i q11 = LDV(r11 + jj), q12 = LDV(r12 + jj), q21 = LDV(r21 + jj), q22 = LDV(r22 + jj); \
__m256i w11 = LDV(r11 + jj + 16), w12 = LDV(r12 + jj + 16), w21 = LDV(r21 + jj + 16), w22 = LDV(r22 + jj + 16); \
_mm256_stream_si256((__m256i *)(o0 + jj), _mm256_add_epi16(q11, q22)); \
_mm256_stream_si256((__m256i *)(o0 + jj + 16), _mm256_add_epi16(w11, w22)); \
_mm256_stream_si256((__m256i *)(o1 + jj), _mm256_add_epi16(q21, q22)); \
_mm256_stream_si256((__m256i *)(o1 + jj + 16), _mm256_add_epi16(w21, w22)); \
_mm256_stream_si256((__m256i *)(o2 + jj), _mm256_add_epi16(q11, q12)); \
_mm256_stream_si256((__m256i *)(o2 + jj + 16), _mm256_add_epi16(w11, w12)); \
_mm256_stream_si256((__m256i *)(o3 + jj), _mm256_sub_epi16(q21, q11)); \
_mm256_stream_si256((__m256i *)(o3 + jj + 16), _mm256_sub_epi16(w21, w11)); \
_mm256_stream_si256((__m256i *)(o4 + jj), _mm256_sub_epi16(q12, q22)); \
_mm256_stream_si256((__m256i *)(o4 + jj + 16), _mm256_sub_epi16(w12, w22)); }
if (copyA) {
if (bpw) { const int h1 = h >> 1;
for (; j + 32 <= h1; j += 32) { MKAC(j) }
o0 += qb - h1; o1 += qb - h1; o2 += qb - h1; o3 += qb - h1; o4 += qb - h1;
}
for (; j + 32 <= h; j += 32) { MKAC(j) }
} else {
if (bpw) { const int h1 = h >> 1;
for (; j + 32 <= h1; j += 32) { MKACS(j) }
o0 += qb - h1; o1 += qb - h1; o2 += qb - h1; o3 += qb - h1; o4 += qb - h1;
}
for (; j + 32 <= h; j += 32) { MKACS(j) }
}
#undef MKACS
#undef MKAC
}
for (int i = 0; i < h; i++) {
const TE *s11 = b11 + (size_t)i * ldb, *s12 = b12 + (size_t)i * ldb;
const TE *s21 = b21 + (size_t)i * ldb, *s22 = b22 + (size_t)i * ldb;
TE *p0 = bpw ? (BS + ((size_t)((i >> hb) << 1)) * qb + (size_t)(i & hm) * (size_t)hph)
: (BS + (size_t)i * (h + PAD));
TE *p1 = p0 + qsx, *p2 = p1 + qsx, *p3 = p2 + qsx, *p4 = p3 + qsx;
TE *c33r = cd[2] + (size_t)i * (h + PAD), *c44r = cd[3] + (size_t)i * (h + PAD);
int j = 0;
#define MKBC(J) { const int jj = (J); \
__m256i q11 = LDV(s11 + jj), q12 = LDV(s12 + jj), q21 = LDV(s21 + jj), q22 = LDV(s22 + jj); \
__m256i w11 = LDV(s11 + jj + 16), w12 = LDV(s12 + jj + 16), w21 = LDV(s21 + jj + 16), w22 = LDV(s22 + jj + 16); \
_mm256_stream_si256((__m256i *)(p0 + jj), _mm256_add_epi16(q11, q22)); \
_mm256_stream_si256((__m256i *)(p0 + jj + 16), _mm256_add_epi16(w11, w22)); \
_mm256_stream_si256((__m256i *)(p1 + jj), _mm256_sub_epi16(q12, q22)); \
_mm256_stream_si256((__m256i *)(p1 + jj + 16), _mm256_sub_epi16(w12, w22)); \
_mm256_stream_si256((__m256i *)(p2 + jj), _mm256_sub_epi16(q21, q11)); \
_mm256_stream_si256((__m256i *)(p2 + jj + 16), _mm256_sub_epi16(w21, w11)); \
_mm256_stream_si256((__m256i *)(p3 + jj), _mm256_add_epi16(q11, q12)); \
_mm256_stream_si256((__m256i *)(p3 + jj + 16), _mm256_add_epi16(w11, w12)); \
_mm256_stream_si256((__m256i *)(p4 + jj), _mm256_add_epi16(q21, q22)); \
_mm256_stream_si256((__m256i *)(p4 + jj + 16), _mm256_add_epi16(w21, w22)); \
_mm256_stream_si256((__m256i *)(c33r + jj), q11); \
_mm256_stream_si256((__m256i *)(c33r + jj + 16), w11); \
_mm256_stream_si256((__m256i *)(c44r + jj), q22); \
_mm256_stream_si256((__m256i *)(c44r + jj + 16), w22); }
#define MKBCS(J) { const int jj = (J); \
__m256i q11 = LDV(s11 + jj), q12 = LDV(s12 + jj), q21 = LDV(s21 + jj), q22 = LDV(s22 + jj); \
__m256i w11 = LDV(s11 + jj + 16), w12 = LDV(s12 + jj + 16), w21 = LDV(s21 + jj + 16), w22 = LDV(s22 + jj + 16); \
_mm256_stream_si256((__m256i *)(p0 + jj), _mm256_add_epi16(q11, q22)); \
_mm256_stream_si256((__m256i *)(p0 + jj + 16), _mm256_add_epi16(w11, w22)); \
_mm256_stream_si256((__m256i *)(p1 + jj), _mm256_sub_epi16(q12, q22)); \
_mm256_stream_si256((__m256i *)(p1 + jj + 16), _mm256_sub_epi16(w12, w22)); \
_mm256_stream_si256((__m256i *)(p2 + jj), _mm256_sub_epi16(q21, q11)); \
_mm256_stream_si256((__m256i *)(p2 + jj + 16), _mm256_sub_epi16(w21, w11)); \
_mm256_stream_si256((__m256i *)(p3 + jj), _mm256_add_epi16(q11, q12)); \
_mm256_stream_si256((__m256i *)(p3 + jj + 16), _mm256_add_epi16(w11, w12)); \
_mm256_stream_si256((__m256i *)(p4 + jj), _mm256_add_epi16(q21, q22)); \
_mm256_stream_si256((__m256i *)(p4 + jj + 16), _mm256_add_epi16(w21, w22)); }
if (copyB) {
if (bpw) { const int h1 = h >> 1;
for (; j + 32 <= h1; j += 32) { MKBC(j) }
p0 += qb - h1; p1 += qb - h1; p2 += qb - h1; p3 += qb - h1; p4 += qb - h1;
}
for (; j + 32 <= h; j += 32) { MKBC(j) }
} else {
if (bpw) { const int h1 = h >> 1;
for (; j + 32 <= h1; j += 32) { MKBCS(j) }
p0 += qb - h1; p1 += qb - h1; p2 += qb - h1; p3 += qb - h1; p4 += qb - h1;
}
for (; j + 32 <= h; j += 32) { MKBCS(j) }
}
#undef MKBCS
#undef MKBC
}
/* LIVENESS: a VERBATIM copy cannot be proven read by `chk` (its values are the same
either way), so this knob corrupts one element per 64 rows in each copy AFTER the
copy is written. chk then moves IFF the copy is on the read path -- and it must
NOT move for an arm that writes the copies and never reads them. */
if (g_lv) {
for (int i = 0; i < h; i += 64) {
cd[0][(size_t)i * (h + PAD)] ^= 1;
cd[1][(size_t)i * (h + PAD)] ^= 1;
cd[2][(size_t)i * (h + PAD)] ^= 1;
cd[3][(size_t)i * (h + PAD)] ^= 1;
}
}
return;
}
// pass 1: the five A sums (4 read streams + 5 write streams; the compiler no longer
// has to keep 18 streams live, which measurably costs ~9 ms at n=4096)
for (int i = 0; i < h; i++) {
const TE *r11 = a11 + (size_t)i * lda, *r12 = a12 + (size_t)i * lda;
const TE *r21 = a21 + (size_t)i * lda, *r22 = a22 + (size_t)i * lda;
TE *o0 = bpw ? (AS + ((size_t)((i >> hb) << 1)) * qb + (size_t)(i & hm) * (size_t)hph)
: (AS + (size_t)i * (h + PAD));
TE *o1 = o0 + qsx, *o2 = o1 + qsx, *o3 = o2 + qsx, *o4 = o3 + qsx;
int j = 0;
#define MKA(J) { const int jj = (J); \
__m256i v11 = LDV(r11 + jj), v12 = LDV(r12 + jj), v21 = LDV(r21 + jj), v22 = LDV(r22 + jj); \
_mm256_stream_si256((__m256i *)(o0 + jj), _mm256_add_epi16(v11, v22)); \
_mm256_stream_si256((__m256i *)(o1 + jj), _mm256_add_epi16(v21, v22)); \
_mm256_stream_si256((__m256i *)(o2 + jj), _mm256_add_epi16(v11, v12)); \
_mm256_stream_si256((__m256i *)(o3 + jj), _mm256_sub_epi16(v21, v11)); \
_mm256_stream_si256((__m256i *)(o4 + jj), _mm256_sub_epi16(v12, v22)); }
#if MK5O
#define MKA2(J) { const int jj = (J); \
__m256i q11 = LDV(r11 + jj), q12 = LDV(r12 + jj), q21 = LDV(r21 + jj), q22 = LDV(r22 + jj); \
__m256i w11 = LDV(r11 + jj + 16), w12 = LDV(r12 + jj + 16), w21 = LDV(r21 + jj + 16), w22 = LDV(r22 + jj + 16); \
_mm256_stream_si256((__m256i *)(o0 + jj), _mm256_add_epi16(q11, q22)); \
_mm256_stream_si256((__m256i *)(o0 + jj + 16), _mm256_add_epi16(w11, w22)); \
_mm256_stream_si256((__m256i *)(o1 + jj), _mm256_add_epi16(q21, q22)); \
_mm256_stream_si256((__m256i *)(o1 + jj + 16), _mm256_add_epi16(w21, w22)); \
_mm256_stream_si256((__m256i *)(o2 + jj), _mm256_add_epi16(q11, q12)); \
_mm256_stream_si256((__m256i *)(o2 + jj + 16), _mm256_add_epi16(w11, w12)); \
_mm256_stream_si256((__m256i *)(o3 + jj), _mm256_sub_epi16(q21, q11)); \
_mm256_stream_si256((__m256i *)(o3 + jj + 16), _mm256_sub_epi16(w21, w11)); \
_mm256_stream_si256((__m256i *)(o4 + jj), _mm256_sub_epi16(q12, q22)); \
_mm256_stream_si256((__m256i *)(o4 + jj + 16), _mm256_sub_epi16(w12, w22)); }
if (bpw) { const int h1 = h >> 1;
for (; j + 32 <= h1; j += 32) { MKA2(j) }
o0 += qb - h1; o1 += qb - h1; o2 += qb - h1; o3 += qb - h1; o4 += qb - h1;
}
for (; j + 32 <= h; j += 32) { MKA2(j) }
#else
for (; j + 32 <= h; j += 32) { MKA(j) MKA(j + 16) }
#endif
for (; j + 16 <= h; j += 16) { MKA(j) }
#undef MKA
for (; j < h; j++) {
TE x11 = r11[j], x12 = r12[j], x21 = r21[j], x22 = r22[j];
o0[j] = x11 + x22; o1[j] = x21 + x22; o2[j] = x11 + x12;
o3[j] = x21 - x11; o4[j] = x12 - x22;
}
}
// pass 2: the five B sums
for (int i = 0; i < h; i++) {
const TE *s11 = b11 + (size_t)i * ldb, *s12 = b12 + (size_t)i * ldb;
const TE *s21 = b21 + (size_t)i * ldb, *s22 = b22 + (size_t)i * ldb;
TE *p0 = bpw ? (BS + ((size_t)((i >> hb) << 1)) * qb + (size_t)(i & hm) * (size_t)hph)
: (BS + (size_t)i * (h + PAD));
TE *p1 = p0 + qsx, *p2 = p1 + qsx, *p3 = p2 + qsx, *p4 = p3 + qsx;
int j = 0;
#define MKB(J) { const int jj = (J); \
__m256i w11 = LDV(s11 + jj), w12 = LDV(s12 + jj), w21 = LDV(s21 + jj), w22 = LDV(s22 + jj); \
_mm256_stream_si256((__m256i *)(p0 + jj), _mm256_add_epi16(w11, w22)); \
_mm256_stream_si256((__m256i *)(p1 + jj), _mm256_sub_epi16(w12, w22)); \
_mm256_stream_si256((__m256i *)(p2 + jj), _mm256_sub_epi16(w21, w11)); \
_mm256_stream_si256((__m256i *)(p3 + jj), _mm256_add_epi16(w11, w12)); \
_mm256_stream_si256((__m256i *)(p4 + jj), _mm256_add_epi16(w21, w22)); }
#if MK5O
#define MKB2(J) { const int jj = (J); \
__m256i q11 = LDV(s11 + jj), q12 = LDV(s12 + jj), q21 = LDV(s21 + jj), q22 = LDV(s22 + jj); \
__m256i w11 = LDV(s11 + jj + 16), w12 = LDV(s12 + jj + 16), w21 = LDV(s21 + jj + 16), w22 = LDV(s22 + jj + 16); \
_mm256_stream_si256((__m256i *)(p0 + jj), _mm256_add_epi16(q11, q22)); \
_mm256_stream_si256((__m256i *)(p0 + jj + 16), _mm256_add_epi16(w11, w22)); \
_mm256_stream_si256((__m256i *)(p1 + jj), _mm256_sub_epi16(q12, q22)); \
_mm256_stream_si256((__m256i *)(p1 + jj + 16), _mm256_sub_epi16(w12, w22)); \
_mm256_stream_si256((__m256i *)(p2 + jj), _mm256_sub_epi16(q21, q11)); \
_mm256_stream_si256((__m256i *)(p2 + jj + 16), _mm256_sub_epi16(w21, w11)); \
_mm256_stream_si256((__m256i *)(p3 + jj), _mm256_add_epi16(q11, q12)); \
_mm256_stream_si256((__m256i *)(p3 + jj + 16), _mm256_add_epi16(w11, w12)); \
_mm256_stream_si256((__m256i *)(p4 + jj), _mm256_add_epi16(q21, q22)); \
_mm256_stream_si256((__m256i *)(p4 + jj + 16), _mm256_add_epi16(w21, w22)); }
if (bpw) { const int h1 = h >> 1;
for (; j + 32 <= h1; j += 32) { MKB2(j) }
p0 += qb - h1; p1 += qb - h1; p2 += qb - h1; p3 += qb - h1; p4 += qb - h1;
}
for (; j + 32 <= h; j += 32) { MKB2(j) }
#else
for (; j + 32 <= h; j += 32) { MKB(j) MKB(j + 16) }
#endif
for (; j + 16 <= h; j += 16) { MKB(j) }
#undef MKB
for (; j < h; j++) {
TE y11 = s11[j], y12 = s12[j], y21 = s21[j], y22 = s22[j];
p0[j] = y11 + y22; p1[j] = y12 - y22; p2[j] = y21 - y11;
p3[j] = y11 + y12; p4[j] = y21 + y22;
}
}
}
// THE SHIP: ordrevall + the 521 dead-code deletion.
static void combine(TE *C, int ldc, const TE *M1, const TE *M2, const TE *M3,
const TE *M4, const TE *M5, const TE *M6, const TE *M7, int h) {
volatile char g_cpad[128] __attribute__((aligned(64))); if (ldc < 0) g_cpad[0] = 0;
const size_t q = (size_t)(h + PAD) * h;
for (int i = 0; i < h; i++) {
const TE *m1 = M1 + (size_t)i * (h + PAD), *m2 = M2 + (size_t)i * (h + PAD), *m3 = M3 + (size_t)i * (h + PAD);
const TE *m4 = M4 + (size_t)i * (h + PAD), *m5 = M5 + (size_t)i * (h + PAD), *m6 = M6 + (size_t)i * (h + PAD);
const TE *m7 = M7 + (size_t)i * (h + PAD);
TE *c11 = C + (size_t)i * ldc, *c12 = c11 + h;
TE *c21 = c11 + (size_t)h * ldc, *c22 = c21 + h;
int j = 0;
#define CMB(J) { const int jj = (J); \
__m256i a1 = LDV(m1 + jj), a2 = LDV(m2 + jj), a3 = LDV(m3 + jj), a4 = LDV(m4 + jj); \
__m256i a5 = LDV(m5 + jj), a6 = LDV(m6 + jj), a7 = LDV(m7 + jj); \
_mm256_stream_si256((__m256i*)(c11 + jj), _mm256_add_epi16(_mm256_sub_epi16(_mm256_add_epi16(a1, a4), a5), a7)); \
_mm256_stream_si256((__m256i*)(c12 + jj), _mm256_add_epi16(a3, a5)); \
_mm256_stream_si256((__m256i*)(c21 + jj), _mm256_add_epi16(a2, a4)); \
_mm256_stream_si256((__m256i*)(c22 + jj), _mm256_add_epi16(_mm256_sub_epi16(_mm256_add_epi16(a1, a3), a2), a6)); }
#if ORD
// lane m4k_cb: E8ORD applied to the combine. SAME loads (14/32 cols), SAME adds (16),
// SAME stores (8), SAME addresses and values -- only the STORE ISSUE ORDER changes:
// each destination row's 64-B line is completed before the next destination is opened
// (the original emits 4 half-lines round-robin). Byte-exact by construction.
for (; j + 32 <= h; j += 32) {
const int jj = j;
if (i + 1 < h) { const int _nr = (h + PAD);
_mm_prefetch((const char*)(M1 + (size_t)(i+1)*_nr + jj), _MM_HINT_T1);
_mm_prefetch((const char*)(M2 + (size_t)(i+1)*_nr + jj), _MM_HINT_T1);
_mm_prefetch((const char*)(M3 + (size_t)(i+1)*_nr + jj), _MM_HINT_T1);
_mm_prefetch((const char*)(M4 + (size_t)(i+1)*_nr + jj), _MM_HINT_T1);
_mm_prefetch((const char*)(M5 + (size_t)(i+1)*_nr + jj), _MM_HINT_T1);
_mm_prefetch((const char*)(M6 + (size_t)(i+1)*_nr + jj), _MM_HINT_T1);
_mm_prefetch((const char*)(M7 + (size_t)(i+1)*_nr + jj), _MM_HINT_T1);
}
__m256i u1 = LDV(m1 + jj), u4 = LDV(m4 + jj), u5 = LDV(m5 + jj), u7 = LDV(m7 + jj);
__m256i v1 = LDV(m1 + jj + 16), v4 = LDV(m4 + jj + 16), v5 = LDV(m5 + jj + 16), v7 = LDV(m7 + jj + 16);
_mm256_stream_si256((__m256i *)(c11 + jj), _mm256_add_epi16(_mm256_sub_epi16(_mm256_add_epi16(u1, u4), u5), u7));
_mm256_stream_si256((__m256i *)(c11 + jj + 16), _mm256_add_epi16(_mm256_sub_epi16(_mm256_add_epi16(v1, v4), v5), v7));
__m256i u3 = LDV(m3 + jj), v3 = LDV(m3 + jj + 16);
_mm256_stream_si256((__m256i *)(c12 + jj), _mm256_add_epi16(u3, u5));
_mm256_stream_si256((__m256i *)(c12 + jj + 16), _mm256_add_epi16(v3, v5));
__m256i u2 = LDV(m2 + jj), v2 = LDV(m2 + jj + 16);
_mm256_stream_si256((__m256i *)(c21 + jj), _mm256_add_epi16(u2, u4));
_mm256_stream_si256((__m256i *)(c21 + jj + 16), _mm256_add_epi16(v2, v4));
__m256i u6 = LDV(m6 + jj), v6 = LDV(m6 + jj + 16);
_mm256_stream_si256((__m256i *)(c22 + jj), _mm256_add_epi16(_mm256_sub_epi16(_mm256_add_epi16(u1, u3), u2), u6));
_mm256_stream_si256((__m256i *)(c22 + jj + 16), _mm256_add_epi16(_mm256_sub_epi16(_mm256_add_epi16(v1, v3), v2), v6));
}
#else
for (; j + 32 <= h; j += 32) { CMB(j) CMB(j + 16) }
#endif
for (; j + 16 <= h; j += 16) { CMB(j) }
#undef CMB
for (; j < h; j++) {
c11[j] = m1[j] + m4[j] - m5[j] + m7[j];
c12[j] = m3[j] + m5[j];
c21[j] = m2[j] + m4[j];
c22[j] = m1[j] - m2[j] + m3[j] + m6[j];
}
}
(void)q;
}
// =====================================================================================
// mm(): C = A*B (Strassen), with the m4k_4cut FUSED-SUM THREADING.
// The node's operand is (A aop A2) -- i.e. a virtual matrix whose rows live in two
// real matrices. At n == 2*CUT the fusion consumes it directly (leaf's pair 2); at
// n == 4*CUT the node does NOT call mk5() and hands each child its sum as a COMPOSED
// operand down this same signature; every generic node has (A2 == 0, aop == OP_RAW).
// =====================================================================================
static void mm(TE *C, int ldc, const TE *A, int lda, const TE *A2, int lda2, int aop,
const TE *B, int ldb, const TE *B2, int ldb2, int bop, int n) {
if (n <= CUT) {
leaf(C, ldc, A, lda, A2, lda2, 0, 0, aop, OP_RAW, B, ldb, B2, ldb2, 0, 0, bop, OP_RAW, n);
return;
}
const int h = n >> 1;
volatile TE g_mpad[16]; if (ldc < 0) g_mpad[0] = 0; /* frame-matched arm: restore mm's 352 B */
const int hp = h + PAD;
const size_t q = (size_t)hp * h;
const size_t qs = q + 2048;
/* THE BLOCK-PACKED SLOT: at h == 2048 each sum slot is FOUR blocks of
(h/2) x (h/2 + PAD) = 4*qb elements. qsx does NOT depend on the knob, so every arm
allocates the same bytes and lands on the same addresses. */
const size_t qb = (size_t)((h >> 1) + PAD) * (h >> 1) + 2048;
const size_t qsx = (h == 2048 || h == 1024) ? 4 * qb : qs;
const int lds = (((h == 2048 && g_kbp >= 1) || (h == 1024 && g_kbp >= 2)) ? (int)((h >> 1) + PAD) : hp);
const int bipA = (g_kbp && h == 1024 && lda == hp);
const int bipB = (g_kbp && h == 1024 && ldb == hp);
const int bip4a = (g_kbp >= 2 && n == 1024 && lda == hp);
const int bip4b = (g_kbp >= 2 && n == 1024 && ldb == hp);
const int ismk5 = (n != 2 * CUT) && (n != 4 * CUT);
unsigned long long _mnt0 = 0;
const int deep = (n < g_topn);
size_t *const op = deep ? &g_soff : &g_off;
TE *AS = (deep ? g_scr : g_pool) + *op;
TE *BS = AS + 5 * qsx;
TE *M = ismk5 ? BS + 5 * qsx : AS;
const size_t alloc = (ismk5 ? 10 * qsx : 7 * qs);
*op += alloc;
const TE *A11, *A12, *A21, *A22, *B11, *B12, *B21, *B22;
if (bipA) { A11 = A; A12 = A + qs; A21 = A + 2 * qs; A22 = A + 3 * qs; }
else { A11 = A; A12 = A + h; A21 = A + (size_t)h * lda; A22 = A21 + h; }
if (bipB) { B11 = B; B12 = B + qs; B21 = B + 2 * qs; B22 = B + 3 * qs; }
else { B11 = B; B12 = B + h; B21 = B + (size_t)h * ldb; B22 = B21 + h; }
if (n == 2 * CUT) {
// deepest Strassen level: each sum feeds exactly one product, so fuse it into
// the consumer's packers instead of materialising 10 sum matrices.
// THE TWO SIDES USE DIFFERENT PATTERNS (that is the Strassen decomposition):
// A-side (X): M1 (0,3,+) M2 (2,3,+) M3 (0,-1) M4 (3,-1) M5 (0,1,+) M6 (2,0,-) M7 (1,3,-)
// B-side (Y): M1 (0,3,+) M2 (0,-1) M3 (1,3,-) M4 (2,0,-) M5 (3,-1) M6 (0,1,+) M7 (2,3,+)
// Under the m4k_4cut threading the mm-level operand V is itself a sum, so each
// leaf operand is re-expressed on the RAW blocks: first pair = sub-blocks of the
// mm-level operand (X, stride lda/ldb), second pair = the same sub-blocks of the
// threaded operand (X2, stride lda2/ldb2), joined by the mm-level op.
// With (A2 == 0, aop == OP_RAW) and (B2 == 0, bop == OP_RAW) this IS the shipped
// call sequence, instruction for instruction.
const TE *C11 = A2, *C21 = A2 ? A2 + (size_t)h * lda2 : 0, *C12 = A2 ? A2 + h : 0;
const TE *C22 = C21 ? C21 + h : 0;
const TE *D11 = B2, *D21 = B2 ? B2 + (size_t)h * ldb2 : 0, *D12 = B2 ? B2 + h : 0;
const TE *D22 = D21 ? D21 + h : 0;
leaf(M + 5 * qs, hp, A21, lda, A11, lda, C21, C11, OP_SUB, aop, B11, ldb, B12, ldb, D11, D12, OP_ADD, bop, h);
leaf(M + 2 * qs, hp, A11, lda, 0, lda, C11, 0, OP_RAW, aop, B12, ldb, B22, ldb, D12, D22, OP_SUB, bop, h);
leaf(M + 1 * qs, hp, A21, lda, A22, lda, C21, C22, OP_ADD, aop, B11, ldb, 0, ldb, D11, 0, OP_RAW, bop, h);
leaf(M + 0 * qs, hp, A11, lda, A22, lda, C11, C22, OP_ADD, aop, B11, ldb, B22, ldb, D11, D22, OP_ADD, bop, h);
leaf(M + 3 * qs, hp, A22, lda, 0, lda, C22, 0, OP_RAW, aop, B21, ldb, B11, ldb, D21, D11, OP_SUB, bop, h);
leaf(M + 6 * qs, hp, A12, lda, A22, lda, C12, C22, OP_SUB, aop, B21, ldb, B22, ldb, D21, D22, OP_ADD, bop, h);
leaf(M + 4 * qs, hp, A11, lda, A12, lda, C11, C12, OP_ADD, aop, B22, ldb, 0, ldb, D22, 0, OP_RAW, bop, h);
combine(C, ldc, M, M + qs, M + 2 * qs, M + 3 * qs, M + 4 * qs, M + 5 * qs, M + 6 * qs, h);
*op -= alloc;
if (n >= 512) MMNA(n, _mnt0);
return;
}
if (n == 4 * CUT) {
// ============ lane m4k_4cut: THE LIFTED FUSION ============
// No mk5(): the five A-sums and five B-sums of this node are NOT materialised.
// Each of the 7 children is a RECURSIVE mm() call whose operand is the sum
// (P_pa1 op P_pa2) of two raw h-blocks, threaded down to the leaves' second
// operand pair. Purely algebraic rewrite of the mk5 form; byte-exact.
const TE *ab[4], *bb[4];
if (bip4a) { ab[0] = A; ab[1] = A + qs; ab[2] = A + 2 * qs; ab[3] = A + 3 * qs; lda = hp; }
else { ab[0] = A; ab[1] = A + h; ab[2] = A + (size_t)h * lda; ab[3] = A + (size_t)h * lda + h; }
if (bip4b) { bb[0] = B; bb[1] = B + qs; bb[2] = B + 2 * qs; bb[3] = B + 3 * qs; ldb = hp; }
else { bb[0] = B; bb[1] = B + h; bb[2] = B + (size_t)h * ldb; bb[3] = B + (size_t)h * ldb + h; }
static const signed char CSQ[7] = {5,1,0,2,4,6,3};
static const signed char PA1[7] = {0, 2, 0, 3, 0, 2, 1};
static const signed char PA2[7] = {3, 3, -1, -1, 1, 0, 3};
static const signed char PAO[7] = {OP_ADD, OP_ADD, OP_RAW, OP_RAW, OP_ADD, OP_SUB, OP_SUB};
static const signed char PB1[7] = {0, 0, 1, 2, 3, 0, 2};
static const signed char PB2[7] = {3, -1, 3, 0, -1, 1, 3};
static const signed char PBO[7] = {OP_ADD, OP_RAW, OP_SUB, OP_SUB, OP_RAW, OP_ADD, OP_ADD};
for (int cc = 0; cc < 7; cc++) {
const int c = (int)CSQ[cc];
const TE *ap1 = ab[PA1[c]], *ap2 = (PA2[c] >= 0) ? ab[PA2[c]] : 0;
const TE *bp1 = bb[PB1[c]], *bp2 = (PB2[c] >= 0) ? bb[PB2[c]] : 0;
mm(M + (size_t)c * qs, hp, ap1, lda, ap2, lda, PAO[c], bp1, ldb, bp2, ldb, PBO[c], h);
}
combine(C, ldc, M, M + qs, M + 2 * qs, M + 3 * qs, M + 4 * qs, M + 5 * qs, M + 6 * qs, h);
*op -= alloc;
if (n >= 512) MMNA(n, _mnt0);
return;
}
MK5B();
mk5(AS, BS, A, lda, B, ldb, h, bipA, bipB);
MK5E(h);
/* ONLY the top node (h == g_ch) may take its raw operands from the copies: at any
other generic node A11/A22/B11/B22 are sub-blocks of THAT node's own operand. */
/* TWO SELECTOR LEVELS. The top node (h == g_ch) uses the SHIP's four copies; every
generic n/2 node (h == g_rh) uses THE RECURSIVE COPY buffer. `wmode` is the WA
ablation: it points the same slots at that node's OWN pool sum arrays instead. */
const int csel = (g_cok && (g_wmode == 3 || (g_wmode == 1 && h == g_ch) ||
(g_wmode == 2 && h != g_ch)));
const int cpsel2 = (g_cok && (h == g_ch || h == g_rh) &&
((h == g_ch) ? g_use : g_use2));
TE *const *const cc = (h == g_ch) ? g_cc : g_rc;
const int sA11 = csel && g_wA11, sA22 = csel && g_wA22;
const int sB11 = csel && g_wB11, sB22 = csel && g_wB22;
const TE *cpA11 = sA11 ? (AS + 3 * qsx) : (bipA ? A11 : (cpsel2 ? cc[0] : A11));
const TE *cpA22 = sA22 ? (AS + 2 * qsx) : (bipA ? A22 : (cpsel2 ? cc[1] : A22));
const TE *cpB11 = sB11 ? (BS + 0 * qsx) : (bipB ? B11 : (cpsel2 ? cc[2] : B11));
const TE *cpB22 = sB22 ? (BS + 4 * qsx) : (bipB ? B22 : (cpsel2 ? cc[3] : B22));
const int clA11 = (cpsel2 || sA11 || bipA) ? hp : lda;
const int clA22 = (cpsel2 || sA22 || bipA) ? hp : lda;
const int clB11 = (cpsel2 || sB11 || bipB) ? hp : ldb;
const int clB22 = (cpsel2 || sB22 || bipB) ? hp : ldb;
mm(AS + 1 * qsx, hp, AS + 1 * qsx, lds, 0, lds, OP_RAW, cpB11, clB11, 0, clB11, OP_RAW, h); // M2
mm(BS + 2 * qsx, hp, cpA22, clA22, 0, clA22, OP_RAW, BS + 2 * qsx, lds, 0, lds, OP_RAW, h); // M4
mm(AS + 3 * qsx, hp, AS + 3 * qsx, lds, 0, lds, OP_RAW, BS + 3 * qsx, lds, 0, lds, OP_RAW, h); // M6
mm(AS + 0 * qsx, hp, AS + 0 * qsx, lds, 0, lds, OP_RAW, BS + 0 * qsx, lds, 0, lds, OP_RAW, h); // M1
mm(BS + 1 * qsx, hp, cpA11, clA11, 0, clA11, OP_RAW, BS + 1 * qsx, lds, 0, lds, OP_RAW, h); // M3
mm(AS + 2 * qsx, hp, AS + 2 * qsx, lds, 0, lds, OP_RAW, cpB22, clB22, 0, clB22, OP_RAW, h); // M5
mm(AS + 4 * qsx, hp, AS + 4 * qsx, lds, 0, lds, OP_RAW, BS + 4 * qsx, lds, 0, lds, OP_RAW, h); // M7
combine(C, ldc, AS + 0 * qsx, AS + 1 * qsx, BS + 1 * qsx, BS + 2 * qsx, AS + 2 * qsx, AS + 3 * qsx, AS + 4 * qsx, h);
*op -= alloc;
if (n >= 512) MMNA(n, _mnt0);
}
static void cpy(TE *dst, int ldd, const TE *src, int lds, int n) {
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j += 16)
_mm256_store_si256((__m256i *)(dst + (size_t)i * ldd + j),
_mm256_loadu_si256((const __m256i *)(src + (size_t)i * lds + j)));
}
void matrix_multiply(int n, const short *A, const short *B, short *C) {
if ((n % 16) || ((unsigned long)C & 31) || ((unsigned long)A & 31) || ((unsigned long)B & 31)) {
for (int i = 0; i < n; i++) // never taken for n = 4096
for (int j = 0; j < n; j++) {
unsigned acc = 0;
for (int k = 0; k < n; k++) acc += (unsigned)(unsigned short)A[(size_t)i * n + k] *
(unsigned)(unsigned short)B[(size_t)k * n + j];
C[(size_t)i * n + j] = (TE)acc;
}
return;
}
volatile int g_mpad2[2]; /* frame-matched: restore matrix_multiply s 40 B */
if (C == 0) g_mpad2[0] = 0;
const int ld = n; // A/B used in place: n%16==0 and A,B,C 32B-aligned, so every lead
// dimension and column offset is still a multiple of 32 bytes.
size_t need = 64, need_scr = 0;
for (int m = n; m > CUT; m >>= 1) {
int hh = m >> 1;
size_t sl = (size_t)(hh + PAD) * hh + 2048;
if (hh == 2048 || hh == 1024) sl = 4 * ((size_t)((hh >> 1) + PAD) * (hh >> 1) + 2048); /* block */
size_t a = ((m == 2 * CUT || m == 4 * CUT) ? 7 : 10) * sl;
need += a;
if (m != n) need_scr += a; /* everything below the top node -> C's buffer */
}
/* the deep arena lives in C, whose usable length is n*n elements; require slack.
And require C's storage to be DISJOINT from A and B: the ROW only ever writes C at
the very end, so C == A would be legal for the row but is NOT legal here. */
const unsigned long ca = (unsigned long)(const void *)C;
const unsigned long aa = (unsigned long)(const void *)A, bb = (unsigned long)(const void *)B;
const unsigned long cb = (size_t)n * (size_t)n * sizeof(TE);
const int split = (need_scr > 0) &&
(need_scr + 4096 <= (size_t)n * (size_t)n) &&
(ca + cb <= aa || aa + cb <= ca) &&
(ca + cb <= bb || bb + cb <= ca);
if (need > 64) {
void *p = 0;
if (posix_memalign(&p, 4096, (split ? need - need_scr : need) * sizeof(TE))) return;
g_pool = (TE *)p;
g_off = 0;
g_topn = split ? n : -1;
{ const size_t _mx = (size_t)n * (size_t)n - need_scr;
const size_t _of = (size_t)3481664 < _mx ? (size_t)3481664 : _mx;
g_scr = split ? ((TE *)C + _of) : g_pool; }
/* lane m4k_layout COMPOSED: Bpanel +4096 B AND arena at 3481664 elems (clamped) */
g_soff = 0;
/* lane m4k_h1024_streams: the four padded copies, in their OWN allocation placed
AFTER the pool so the pool's size -- and therefore its mmap address -- is
UNCHANGED (a bigger pool lands lower under top-down mmap). */
void *pc = 0;
const size_t htop = (size_t)n >> 1;
const size_t qc = ((htop + PAD) * htop) + 2048; /* +1 page slot stagger */
g_ch = (int)htop;
g_cok = (htop % 32 == 0); /* the copy loops cover h%32==0 */
if (g_alloc && g_cok && posix_memalign(&pc, 4096, 4 * qc * sizeof(TE)) == 0) {
g_cc[0] = (TE *)pc; g_cc[1] = g_cc[0] + qc; g_cc[2] = g_cc[1] + qc; g_cc[3] = g_cc[2] + qc;
} else { g_cc[0] = g_cc[1] = g_cc[2] = g_cc[3] = 0; g_cok = 0; }
/* THE RECURSIVE COPY: ONE reused buffer for the four raw quadrants of every generic
n/2 node (h == g_rh). The nodes run one at a time (depth-first, siblings run
sequentially), so a single buffer is live at any moment; that is what keeps the
first-touch cost at 8.7 MB against the 58.7 MB of stores those seven nodes issue. */
void *pr = 0;
const size_t hr = htop >> 1;
const size_t qr = ((hr + PAD) * hr) + 2048;
g_rh = (int)hr;
if (g_alloc && g_rh > 0 && g_rh % 32 == 0 && posix_memalign(&pr, 4096, 4 * qr * sizeof(TE)) == 0) {
g_rc[0] = (TE *)pr; g_rc[1] = g_rc[0] + qr; g_rc[2] = g_rc[1] + qr; g_rc[3] = g_rc[2] + qr;
} else { g_rc[0] = g_rc[1] = g_rc[2] = g_rc[3] = 0; }
mm(C, n, A, ld, 0, ld, OP_RAW, B, ld, 0, ld, OP_RAW, n);
free(pr); g_rc[0] = g_rc[1] = g_rc[2] = g_rc[3] = 0; g_rh = 0;
free(pc); g_cc[0] = g_cc[1] = g_cc[2] = g_cc[3] = 0; g_ch = 0; g_cok = 0;
free(p);
} else {
leaf(C, n, A, ld, 0, ld, 0, 0, OP_RAW, OP_RAW, B, ld, 0, ld, 0, 0, OP_RAW, OP_RAW, n);
}
}
#pragma GCC pop_options
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 620.396 ms | 122 MB + 1016 KB | Accepted | Score: 100 | 显示更多 |