// lane m164 arm dmb64: artifact + 64 MB reserved before the arena
// =====================================================================================
// lane m4k_h1024_streams -- PAD THE INPUT: the four RAW operand blocks of the top node are
// delivered as packed copies, so that every read below the top node uses a 4-KB-break
// row stride instead of the judge's ld = n.
//
// WHY (measured, judge, in situ). The h=1024 level's own `mk5` bracket reads its four
// source streams at 133-145 ticks/MB when the operand is PACKED (row stride h+PAD) and at
// 175 ticks/MB when it is a 1024x1024 sub-block of the judge's A/B (lda = 4096 elements =
// 8192 B rows, i.e. 4-KB-CONGRUENT and 50 % dense). That is the row's own PAD law ("PAD=0
// +27.6 ms") applied to the one operand the row never pads. A WA ablation that replaced
// the top node's four RAW operands (A11, A22, B11, B22) with packed arrays of the same
// volume moved the ROW by -25.29 M ticks (-7.02 ms) while leaving the h=2048 bracket
// BIT-FOR-BIT unchanged -- an internal control -- so the effect is the operand LAYOUT, not
// volume.
//
// THIS IS BYTE-EXACT AND IT COSTS THE COPIES RATHER THAN A RE-READ: `mk5` at the top node
// already loads a11 = A11 and a22 = A22 (pass 1) and b11 = B11, b22 = B22 (pass 2) -- the
// values are in registers -- so it writes four packed h x (h+PAD) copies at ZERO extra
// reads. The top node's four raw-operand children then take their operand from the copy
// with ld = h+PAD; the copy's pad columns are never written and never read (every consumer
// takes its lead dimension from the array it is given), so the result is unchanged.
//
// GATE (judge `custom_test`, 3 interleaved rounds, `Pending 数量` 0 before every call,
// arms never overlapping, all three arms .text-BYTE-IDENTICAL -- md5 951a47d0a2, 523 635 B
// -- and differing only in ONE inert .rodata byte, the arm knob):
// plain (g_copy=0 g_use=0) 2 244 104 358 / 2 244 331 310 / 2 244 872 224 spread 0.034 %
// ctrl (g_copy=1 g_use=0) 2 254 690 162 / 2 254 667 990 / 2 255 061 754 spread 0.017 %
// cand (g_copy=1 g_use=1) 2 236 987 480 / 2 236 827 442 / 2 236 785 886 spread 0.009 %
// cand - ctrl = -17.88 M = the pure gain; ctrl - plain = +10.56 M = the copies' cost;
// cand - plain = -7.32 M TICKS = -2.03 ms, and cand's WORST round is 7.1 M below plain's
// BEST. chk 3684834600764090267 at n=4096 for ALL THREE ARMS.
// =====================================================================================
// =====================================================================================
// lane m4k_place -- THE ARENA'S SLOT STRIDE IS THE ONLY PLACE THE ROW'S ARENA LAYOUT IS
// FREE, AND +4096 B PER SLOT IS A MEASURED WIN: -0.917 ms (judge custom_test, 4
// interleaved rounds, control spread 0.010 %, candidate spread 0.017 %, 4/4 with a
// complete max-candidate < min-control separation).
// WHY: every sibling block in the arena (5 A-sums, 5 B-sums, 7 M) sits at a stride of
// exactly q = (h+PAD)*h elements, i.e. 2080 / 528 / 136 / 36 PAGES -- so every block is
// 4-KB-congruent with every other, and their L2/L3 set indices coincide. +2048 elements
// (+1 page) per slot moves each block's L2/L3 set index by 64 while LEAVING the L1 set
// index and the 4-KB-alias class untouched (4096 = 2^12) -- so this is a pure L2/L3
// set-mapping move, not an aliasing change.
// THE ALIASING AXIS ITSELF IS MEASURED NULL (a +64 B/slot arm preserves the phase
// separation but breaks only the 4K class: +0.032 ms), and the ROW-STRIDE pad is at its
// MEASURED OPTIMUM (PAD=32: PAD=0 +27.6 ms, PAD=16 +10.2, PAD=48 +10.4, PAD=64 +6.8,
// PAD=128 +? -- the 4-KB-alias-breaking pad is worth 27.6 ms and 32 is the best of them).
// BYTE-EXACT: the arena is private scratch; every block is fully written before it is
// read and no value depends on where a block sits. chk 3684834600764090267 at n=4096.
// =====================================================================================
// work/m4k_rescreen arm d_q_s46 : Q=[5, 2, 4, 0, 1, 3, 6] BAND LEADER Q ALONE (re-based)
// work/m4k_fp arm alFP : w45 + (i) the 7 products live in the 10 SUM slots (aliased pool)
// + (ii) m4k_joint's P order {1,3,5,0,2,4,6} at the generic mk5 nodes (JUDGE_MEM 551).
// =====================================================================================
// mmms4k -- LANE m4k_4cut: LIFT THE SUMS FUSION FROM n == 2*CUT TO n == 4*CUT.
//
// Base: problems/mmms4k/sol_pc_cp128_tk64.cpp (md5 c022f9ff741943a576d316ec65835058,
// sid 92062, board 656.076 ms) = the established layout row.
//
// THE MECHANISM. `mm()` at n == 2*CUT already fuses its 5 A-sums and 5 B-sums into
// the CONSUMER'S tr/packer builds -- it never materialises them -- and L4 is the only
// level with no sums and the only node whose footprint fits the 6 MB L3 (2.49 MB), and
// it measures 23 % DRAM-served against L1-L3's 92/82/100 %. This build lifts exactly
// that fusion one level: at n == 4*CUT the node skips `mk5()` entirely and threads the
// sum's operands (A2, lda2, aop) down through the recursive `mm()` to the leaves.
//
// THE ALGEBRA THAT MAKES IT EXACT (and why no 4-term machinery is needed beyond a
// second operand pair). Let the n==4*CUT node's A be split into 512-blocks A_p and
// each A_p into 256-blocks, and let a child's virtual operand be V = A_p1 aop A_p2.
// Its leaf Q's operand is V_q1 xop V_q2, and since all four ops are elementwise
// additions/subtractions in Z/2^16,
// (A_p1,q1 aop A_p2,q1) xop (A_p1,q2 aop A_p2,q2)
// = (A_p1,q1 xop A_p1,q2) aop (A_p2,q1 xop A_p2,q2)
// i.e. every fused leaf operand is (pair) aop (pair) with the SAME within-pair op xop
// on both pairs. So `leaf` needs ONE extra operand pair and ONE extra op, and when
// (A2 == 0, aop2 == OP_RAW) it is bit-for-bit the pre-existing 2-term build.
// operand_A = (A aop A2) aop2 (A3 aop A4) [aop2 == OP_RAW => pair 2 absent]
// A, A2 have row stride `lda`; A3, A4 have row stride `lda2`.
// The B side is symmetric: the packer gains the same second pair.
//
// WHAT IS REMOVED: 49 `mk5()` calls at h=512 per call (measured by DELETION at
// -93.7 M ticks, see rm_mk5_4cut.cpp / FINDINGS.txt).
// WHAT IS ADDED: the leaves now read raw A/B blocks with 2x the sub-block multiplicity
// (24 instead of 12 sub-block reads per child for the 5 sum-children).
// =====================================================================================
/* lane m4k_mk5: MK5O -- mk5's five NT sum writes complete each destination's 64-byte line
before opening the next destination, instead of emitting five half-lines round-robin
(which holds five write-combining buffers half-full at once). Same loads, adds,
stores, addresses and values; the store ISSUE ORDER is the only change, so it is
byte-exact by construction. Independent of lane m4k_place's slot-stride edit below/
above: that moves WHERE the slots are, this changes the ORDER they are written in.
Judge custom_test on sol_csq_cq1_tailmix: -9.64 M ticks = -2.68 ms, 3 interleaved
rounds, chk identical, arms never overlapped; board paid -3.006 ms (631.820->628.814). */
#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;
#ifndef DMB
#define DMB 64
#endif
#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;
__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;
/* lane m164 LAYOUT PROBE: reserve DMB megabytes (never touched) BEFORE the arena, so the
pool's virtual address moves. Pure placement; no value depends on it. */
if (DMB) { void *dv = 0; if (posix_memalign(&dv, 4096, (size_t)DMB << 20)) dv = 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 | 614.678 ms | 122 MB + 1020 KB | Accepted | Score: 100 | 显示更多 |