/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126197 (unit-interleaved scheduled integer kernel); all inherited citations and license notes retained below.
Idea: Bias the A cursor by 256 bytes and B cursor by 128 bytes, compensating every assembly load and prefetch displacement. All B data loads and twice as many A-row loads can then use signed-byte displacements; the loop's pointer increments and arithmetic remain identical.
Purpose: Reduce instruction bytes in both the initialized group and repeated group without adding loop work or changing the packed representation. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126173 (interleaved A representation and delayed-consumer loop); all original author and license references retained below.
Idea: Apply the main loop's delayed-consumer instruction schedule to the peeled initial group as well. Direct the first low and first high products into their accumulators and remove their corresponding additions, preserving direct initialization. The other three rotations now use the same latency-overlapped order as the fifteen main groups.
Purpose: Remove the old synchronous initialization schedule from one sixteenth of every microkernel without adding zeroing or changing exact modulo-32 arithmetic. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126165 baseline; full previous attributions and license information below.
Idea: Interleave the four A rows inside each eight-K replicated unit. The four 128-byte units become one adjacent 512-byte block instead of separate 2240-byte row streams. Remap each row displacement and advance A by 512 bytes per eight-K group; all low/high values and madd/add scheduling are identical.
Purpose: Improve locality and hardware stream detection while retaining every packed unit's alignment and exact modulo-32 arithmetic. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126131,
fixed-trip modular microkernel and Morton engine; inherited references below.
Idea: The microkernel assembly only reads packed A/B and returns four vectors
in declared output registers; C stores occur in ordinary intrinsics afterward.
Declare Apk and Bpk as memory inputs, retaining every register clobber, instead
of claiming that the assembly may overwrite arbitrary memory.
Purpose: Give the compiler accurate conservative dependencies so unrelated
addresses and transform state need not be invalidated around each dot kernel.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126049,
flat-Morton integer engine; all original references and licenses retained below.
Idea: For large h>=1024 in-place S/T transformations, issue each complete cache
line with streaming stores after its old data has been loaded. Retain cached
stores on smaller operands, contiguous traversal for equal layouts and full128
tile conversion for mixed layouts. Fence before recursive consumers can read it.
Purpose: Test reduced dirty-cache traffic for transformations whose paired
large operands exceed the shared cache, while preserving all modulo32 results.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126007,
complete-tile Morton traversal and fixed-trip kernel; inherited authors and
license notes retained below.
Idea: If the source and destination operands are both Morton packed, their
physical orders coincide for every scalar element. Flatten the entire addition
or mutation into one contiguous vector pass, eliminating tile decoding, row
bounds and strided prefetch addresses. Keep complete-tile conversion for mixed
row-major/Morton paths and retain the existing store-type threshold.
Purpose: Exploit identical layouts to remove all two-dimensional indexing from
operand transformations that need no actual layout conversion.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/125995,
current packed integer engine; inherited references and license notes below.
Our https://duck.ac/submission/125989 inspired using the kernel's fixed
iteration count to remove the dynamic prologue guard.
Idea: All call sites pass NPAIR=64. Keep the direct-initializing first group,
then execute exactly fifteen groups, with one countdown per group. Remove the
initial decrement and impossible short-panel branch without changing arithmetic.
Purpose: Reduce per-microkernel control work for the fixed KC=128 packing.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/125988:
tile-ordered Morton operand transformations, retaining original sources below.
No independent license declared.
Idea: Apply the complete-tile traversal to 256-square transformations as well.
Each tile has 64 KiB per operand; two sources plus a destination fit in the
256 KiB L2 when visited tile by tile, whereas the old per-row two-half
traversal interleaves both horizontal tiles. Preserve cached stores at this
size and the same Morton mapping and arithmetic.
Purpose: Extend locality gains to the frequent near-leaf transformations.
*/
/* References: saffah_cc_v41_agg1 https://duck.ac/submission/124010:
Morton-128 recursive engine; our https://duck.ac/submission/125847 supplies
compact panels and shared B blends. All original citations retained below;
no independent source license declared.
Idea: Traverse large Morton-destination operand additions and mutations one
complete 128x128 tile at a time. Compute tile bases once, then advance rows
with fixed pitches, replacing the row-first traversal that visits every tile
before returning to the next row. Preserve source layout interpretation,
operation signs, exact wraparound arithmetic, and existing NT-store threshold.
Purpose: Improve destination/TLB locality and remove repeated Morton index work.
*/
/* References: saffah_cc_v41_agg1 https://duck.ac/submission/124010:
reused the reordered madd/add kernel, Morton temporary layout, interleaved
products and recursive engine; full original attribution follows.
saffah_codex_6a_agg3 https://duck.ac/submission/121860 and
https://duck.ac/submission/122104: own 2240-byte A panels, shared intermediate
B blends, and optional huge-page advice on the owned arena. No separate license.
Idea: Combine the newer recursive layout with smaller A-panel spacing and
factor each four-output B diagonal permutation through four common blends,
removing eight blend instructions per group. Apply the existing huge-page hint
to this larger private arena, with unchanged correctness if unavailable.
Purpose: Reduce panel cache footprint, B-packing arithmetic, and TLB pressure.
*/
// ===== REFERENCES =====
// [1] duck.ac 本账号 saffah_cc_v41_agg1，提交 #120597 <https://duck.ac/submission/120597>
// 用途：**本件的正文本体**（本账号 mmmi4k 现役最好件 `mi4z_sub1.cpp`，官判 1723.084443 ms，
// 排行榜本账号第 1）。除下列唯一改动外，正文与该件逐字节相同。
// 继承其全部在档引用：[2] #120472（本账号，mmmi2k 现役引擎移植到 n=4096 的发放件）、
// [3] #117857（本账号，4 层 Strassen + `tiled2_128` 叶 + `vpmaddwd` 微内核 + `packA2`/`packB2`
// + `mi2ap_` 重指/缩映射的正文谱系）、[4] `saffah_codex_6a_agg3` #120382（n=4096 判据，
// = 我方 #117857 正文 + 恰 3 处 token）、[5] #101322/#101332（本账号，`mmap(MAP_POPULATE)` 否证）。
// [6] 判题机「立即体验」定价（本席 mi4r_，2026-10-03，`problems/mmmi4k/work/mi4r__tmp/mi4r_tr.py`，
// n=2048，3 轮轮转臂序，同码孪生臂 dup 作噪声底）：同批 8 臂逐位等价（chk=076457ea50ae29b5 全同），
// `ntab1026`（本发改动，h≥1024 才走 NT）= **−0.238%**（−2.08 M 核周期）；`ntab514` = −0.111%；
// `ntab`（全 h）= +1.276%；`ntc`（mikro 的 C 片改 NT）= +2.340%；`ntpk`（packB2 改 NT）= +4.455%。
// 用途：本发唯一改动的**直接实验依据**。
// [7] 判题机「立即体验」定价（本席 mi4ad_，2026-10-03，problems/mmmi4k/work/mi4ad__tmp/
// mi4ad_tr.py + mi4ad_gen.py，n=2048，轮转臂序，in-program `cyc` 口径，
// 每批同码孪生臂 dup 作噪声底 ≤0.024%）：本发的唯一改动 = `mutate_blk` 的
// 跨行源预取，**8 次独立批次的读数 −0.0051%(dup) / −0.1272% / −0.0975% /
// −0.0996% / −0.1275% / −0.0909% / −0.1192% / −0.1184%**（7 次同号，均值 ≈ −0.111%）；
// ★ 纯布局对照臂（同指令、同码位，只把目标地址换成「本行自身」= 无新流量）= **−0.0066%**
// ⇒ 效应可归因到预取本身，不是码位/对齐伪影。同批标定臂（把 NT 阈值抬到 h≥2048，
// 即在 n=2048 上关掉现役 h=1024 的 NT 支）= **+0.2400%**，与在册 mi4r_ 的 +0.238% 逐位吻合。
// [8] 判题机「立即体验」定价（本席 mi4av_，2026-10-05，`problems/mmmi4k/work/mi4av__tmp/`）：
// 本件唯一改动 = 在 #123988 之上再加一条**逐位等价的 packB2 建对改写**（见下 思路①），
// 全程序 n=2048 通道轮转 3 rep：base 867,719,592 ↔ 本件 864,317,312 中位 = **-0.392%**；
// 同批 zapint（= #123988）865,252,560 = -0.284% ⇒ 本发改动的**边际 = -0.108%**（3 rep 内极差
// 仅 0.001%，同批 chk 逐位同 076457ea50ae29b5）。
// [9] 在册 `zapint` 合成件：`work/mi4aq__tmp/u_zapint.cpp`（mi4aq_ 封条 (2077) 建、mi4ar_ (2080)
// 原样重建、mi4au_ (2098) 逐类拆账）。本件正文 = 该件**去掉仪器后的正文** + 本席的 packB2 改写；
// 本席以同一 patch 链重建，`work/mi4aq__tmp/u_subzapint.cpp` 与重建件 `cmp -s` 逐字节相同。
// [10] #123988（本账号，官判 1711.963739 ms）= 把在册 `zapint` 发射出去的实测读数：判题机 n=4096
// 只兑现 **-0.128%**，而在册免费通道（子问题 str=4096 与全程序 n=2048）读数是 -0.38 ~ -0.41%
// ⇒ **本题免费通道对「打包布局类」改动系统性高估 ≈2.2x**，后席定价必按此折算。
// [11] 判题机「立即体验」免费通道定价（本席 mi4cd_，2026-10-04，工件
// `problems/mmmi4k/work/mi4cd__tmp/{mi4cd_gen.py,mi4cd_u3.cpp,mi4cd_art_gen.py,
// mi4cd_art_R.cpp,mi4cd_art_P.cpp}`）：**本发唯一改动的实验依据**。
// ① 装置口径（判题机，轮转臂序 64 轮，同码孪生分置两处位置）：
// `R`（在役 mikro）= 34898 / 34937，`PR`（本发重排）= 34515 / 34500
// ⇒ **−1.17%**（孪生互差 0.04%）；同批证书 `mix` 33500 / `mixDep` 35187 / `a96` 33437
// 复现 (2199)(2861) 在册值。
// ② 真件口径（判题机，**整个引擎** n=2048，自计时 min，参考臂 M 定标）：
// `art_R`（在役）858 542 998 / 858 562 870 / 858 658 558；
// `art_P`（只做读后置）856 289 452 / 856 315 200 / 856 315 432 ⇒ −0.269%；
// `art_W`（**本件** = 读后置 + bp 载入前置）**855 310 270 / 855 314 428 ⇒ −0.377%**；
// 参考臂 M 8220/8224（差 0.05%）⇒ 机器状态同；三件 C 的 FNV **逐位相同**
// （7de45bf4ee0d3f0c）⇒ 逐位等价。
// ③ 逐位闸门：本机四档 FNV 与在役同；判题机 mikro 值闸门 14397103 两臂同。
// ===== 思路 =====
// 【mi4cd_ 席 · 本发 = 唯一改动：**`mikro` 稳态环体 96 条 ALU 指令的纯次序重排**（pushR）】
// ① 机理（A：读后置）：环体 4 个 chunk 各 2 载入 + 12 `vpmaddwd` + 12 `vpaddd`。每条 `vpaddd` 与它所要读的
// `vpmaddwd` 在原件里只相隔 4 条指令，而 vpmaddwd 的延迟是 5 拍 ⇒ 每条 add 都在 RS 里空等，
// 而 p0+p1（48 madd + 16 add = 64 槽 / 2 每拍）与 p5（32 add / 1 每拍）**同时 100% 饱和**、
// 零余量 ⇒ 任何空等等不回来。本发把每条 add **推迟到「同一临时寄存器 ymm10..15 的下一次
// 写入」之前**：读取仍发生在被读值写入之后、覆写之前 ⇒ 数据流逐条不变；madd 序、8 条
// `vmovdqu`、`prefetcht0`、环控指令**一字不动**。madd→add 距离 4 → 平均 **11.1**（max 13）。
// ② 逐位等价：模 2^32 加法可交换结合 ⇒ 每个累加器收到的加数集合不变、只次序变；另有本席的
// 寄存器数据流模拟器逐条核对（每条 add 读到的「被读值身份」+ 目标累加器与原件逐条一致），
// 且指令多重集逐条守恒（指令数、码长、`.text` 大小全不变：两件 `size` 同为 32469/664/92608）。
// ③ 价：装置口径 −1.17%（mikro），真件口径 −0.265%/−0.273%（整程序 n=2048）——两者自洽
// （n=2048 时 mikro 约占整程序 23%）。
// ④ 机理（B：bp 载入前置）：8 条 `vmovdqu` 的「载入→首次使用」距离在原件里是 2/9，
// 即 chunk 首条 madd 要等一条 L1 载入。本发把每条 `vmovdqu` 搬到「它之前最后一次读该
// 寄存器之后」的最早合法位置 ⇒ 纯提前，写仍在所有旧读之后、在新读之前。距离 2 → 19。
// ⑤ 判题机真件口径（n=2048，两次独立发，参考臂 M 定标）：R 858 542 998 ↔ 本件
// 855 310 270/855 314 428 ⇒ **−0.377%**（无 B 的 PR 版为 856 289 452/856 315 200 ⇒ −0.269%；
// B 的边际 = **−0.108%**）；C 的 FNV 两件逐位同 7de45bf4ee0d3f0c。
// ⑥ 否证面（本席与在册）：mikro 环内 15 件重排/字距臂全负；本发与它们的差别是「**不动寄存器
// 分配、不改 uop 混比，只增 madd→add 距离**」。同批剂量证书把距离拉到 21（`d4T11` = 34531）
// 仍只与本发（34515/34500）同级 ⇒ **距离收益本发已吃满**；"溢出累加器换寄存器"被载入口
// 卡死（56 条折叠载入已占 28/32 拍，每条溢出累加器每 body 再加 ≥4 条载入），无余量。
// 【mi4av_ 席 · 本发 = 在 #123988（`zapint`：Z128 A/B 打包象限 + depth-1 M 区行交织）之上，
// 再加唯一一条改动：**`packB2` 建对改用 `vpblendw`**（`vpshufb` 掩码恒等式，逐位等价）。】
// ① 本发唯一改动（`packB2` 的 `F11_PAIR`）：原式
// LO = or(shuf(x0,BLA), shuf(x1,BLB)); HI = or(shuf(x0+half,BHA), shuf(x1+half,BHB));
// 经微测逐字节验证的恒等式改写为
// LO = vpblendw(x0, vpslld(x1,16), 0xAA); HI = vpblendw(vpsrld(x0+half,16), x1+half, 0xAA);
// （BLA/BLB/BHA/BHB 四张表本来就只是「取每个 dword 的低 16 位再交错两行」，故恒等；
// 已用独立微测 LO/HI 两式 memcmp 全等 ✓）。4 条 `vpshufb` + 2 条 `vpor`（其中 4 条占 p5）
// 每对 → 0 条 p5 uop，每对少 2 条 uop ⇒ 每 32 B 输出块 **-8 uop、-16 条端口 5 uop**。
// 8 条 store 与 24 条 `vpblendd` 一条不动。
// ② 判题机价（全程序 n=2048，轮转 3 rep，同批 chk 逐位同）：base 867,719,592 · zapint 865,252,560 ·
// 本件 864,317,312（中位）⇒ 本改动**边际 -0.108%**（vs zapint），本件合计 -0.392%（vs base）。
// `objdump`：`packB2` 342 → 321 条指令，**端口 5 类（vpshufb/vpermd/vpermq/vpunpck）32 → 0**。
// ③ 逐位等价（本席本机，四档全元素 FNV 与基座同值）：
// n=512/1024/2048/4096 = f045914c11ac8a83 / 99c52671a6c29e65 / 076457ea50ae29b5 / 471d8c5eaedc75a1。
// ④ 否证面（在册，勿重走）：`packA2` 已贴端口地板（每 32 元素 = 16 store-data(p4) + 16 p5）；
// `packB2` 那 24 条 `vpblendd` 的 4 源逐道 round-robin 是**三下界**（L 块 k 槽 m 必须取
// (m+k) mod 4 行对，一次 32 B 载入只能给一行）⇒ 本发**不碰 blend**，只削建对的 uop 与端口压力。
// ⑤ 风险面：纯恒等改写，不增删算术；四档 FNV + 同批 chk + 微测三重背书。
// 【mi4ad_ 席 · 本发 = 唯一改动：`mutate_blk` 的行环内加一条「跨行源预取」
// `__builtin_prefetch(src + i*lds + lds, 0, 3)`，其余正文与 #123136 逐字节相同。】
// ① 动机（本题新量成的「strided 读是延迟型」事实）：`rec_w` 的 4 条 `mutate_blk`
// （S2/T2/S4/T4，判题机口径 ≈50 ms）读的是 A11/A12/B22/B21 四个**16 KB 跨步的象限行**
// （h=256 级一行仅 1 KB），而它是**全程序唯一没有预取的跨步读通道**：`addblk` 每行有 2 条
// `__builtin_prefetch`，`tiled2_128` 有 8 条，mikro 的 B 面板有 `prefetcht0 768`
// ——只有 `mutate_blk` 是零预取（在册「逐循环跨页访存条数 vs prefetch 条数」表内明列 0）。
// 本席把它补上：**每行 1 条、提前 1 行**（= 提前 lds = 16 KB 地址，h=256 时 ≈146 拍处理时间的深度）。
// ② 判题机定价（`cyc` 口径 ±0.01%，见引用 [7]）：**7 批同号、均值 −0.111%**；对照臂
// `mut2`（提前 2 行）= −0.1018%、`mutloc2`（locality 2）= −0.1011%（均略差于提前 1 行 locality 3）
// ⇒ 深度已扫，取「提前 1 行」。
// ③ 否证面（本席同批实测，全部判负，勿重走）：预取族的**其余站点全部无正号** ——
// `addblk` 提高 2 行 +0.0232% · 删 `addblk` 的 Y 预取 +0.0632% · `packB2` 加 1 行 +0.0141% /
// 加整组 8 行 +0.0979% · `comb_win` 加 7 条下一行 +0.4183%（locality 1 同样 +0.4097%）·
// `comb_all` 加 7 条 +0.1398% · `comb_win_interleaved` 加 1 条 +0.0164% · `tiled2_128` 把
// A 面板预取提前 2 个 ic 块 +0.1226% / locality 降 1 −0.0029% · 叶片 C 行预取 +0.2441%
// ⇒ **本题的预取余额只在这一处**。
// ④ 逐位等价：n=512/1024/2048/4096 四档全元素 FNV 与 #123136 逐位同
// （f045914c11ac8a83 / 99c52671a6c29e65 / 076457ea50ae29b5 / 471d8c5eaedc75a1）；
// `.text` 14316 B 与基座同、`mutate_blk` 两个克隆各 +3 条（1 预取 + 2 地址）、
// `tiled2_128`/mikro 的 `.p2align 5` 热环未动。风险面：预取只读、不改语义；
// 末行 (i=h-1) 预取地址落在母矩阵内（lda ≥ h），非越界。
#include <immintrin.h>
#include <string.h>
#include <stdint.h>
#include <stdlib.h>
#include <sys/mman.h>
#define MR 4
#define NCG 1
#define KC 128
#define MC 4
#define NC 128
#define NPAIR (KC / 2)
#define APSTR 560 /* 2048 B per row = 16 units x 128 B */
#define NR (8 * NCG)
#define BASE 128
typedef int32_t T;
#define AVX2 __attribute__((target("avx2")))
#define LD(p) _mm256_load_si256((const __m256i *)(p))
#define LDA(p) _mm256_load_si256((const __m256i *)(p))
#define STA(p, v) _mm256_store_si256((__m256i *)(p), (v))
#define ST(p, v) _mm256_store_si256((__m256i *)(p), (v))
#define ZQ 1
static inline size_t zoff(int h, int i, int j) {
size_t o = 0;
if (h >= 2048) { o += ((size_t)((i>>10)*2 + (j>>10))) << 20; i &= 1023; j &= 1023; }
if (h >= 1024) { o += ((size_t)((i>>9)*2 + (j>>9))) << 18; i &= 511; j &= 511; }
if (h >= 512) { o += ((size_t)((i>>8)*2 + (j>>8))) << 16; i &= 255; j &= 255; }
if (h >= 256) { o += ((size_t)((i>>7)*2 + (j>>7))) << 14; i &= 127; j &= 127; }
return o + (size_t)(i<<7) + (size_t)j;
}
static int32_t Apk[(MC + 4) * APSTR + 64] __attribute__((aligned(64)));
static int32_t Bpk[16896] __attribute__((aligned(64))); /* [jg][m][B0|B1] */
/* ===== mi2ap_：P 与前两级递归帧重指到调用方输出缓冲 C，并把私有 arena 的映射一并缩小 ===== */
static T *g_pool = 0; /* == C：深度-0 的重指池 */
static T *g_deep = 0; /* 最深一帧（h == BASE），留在 scr 内 */
static int g_rel = 0; /* 1 iff 本次调用启用重指与缩映射（仅顶层 n == 2048） */
/* MIKRO-BEGIN */
#define AVX2 __attribute__((target("avx2")))
AVX2 static inline __attribute__((always_inline)) void mikro(int npair, const int32_t *ap, const int32_t *bp,
int32_t *c, int ldc, int first) {
(void)npair;
static_assert(KC==128, "fixed packing kernel");
int groups=KC/8-1;
ap+=64; bp+=32;
__m256i q0, q1, q2, q3;
__asm__ volatile(
"vmovdqu -128(%[bp]), %%ymm8\n\t"
"vmovdqu 0(%[bp]), %%ymm9\n\t"
"vpmaddwd -192(%[ap]), %%ymm8, %%ymm4\n\t"
"vpmaddwd -64(%[ap]), %%ymm8, %%ymm5\n\t"
"vpmaddwd 64(%[ap]), %%ymm8, %%ymm6\n\t"
"vpmaddwd 192(%[ap]), %%ymm8, %%ymm7\n\t"
"vpmaddwd -256(%[ap]), %%ymm9, %%ymm12\n\t"
"vpmaddwd -128(%[ap]), %%ymm9, %%ymm13\n\t"
"vpmaddwd 0(%[ap]), %%ymm9, %%ymm14\n\t"
"vpmaddwd 128(%[ap]), %%ymm9, %%ymm15\n\t"
"vmovdqu 64(%[bp]), %%ymm9\n\t"
"vpmaddwd -256(%[ap]), %%ymm8, %%ymm0\n\t"
"vpmaddwd -128(%[ap]), %%ymm8, %%ymm1\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 0(%[ap]), %%ymm8, %%ymm2\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 128(%[ap]), %%ymm8, %%ymm3\n\t"
"vmovdqu -64(%[bp]), %%ymm8\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpmaddwd -184(%[ap]), %%ymm8, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpmaddwd -56(%[ap]), %%ymm8, %%ymm15\n\t"
"vpmaddwd 72(%[ap]), %%ymm8, %%ymm10\n\t"
"vpmaddwd 200(%[ap]), %%ymm8, %%ymm11\n\t"
"vpmaddwd -248(%[ap]), %%ymm9, %%ymm12\n\t"
"vpmaddwd -120(%[ap]), %%ymm9, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 8(%[ap]), %%ymm9, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 136(%[ap]), %%ymm9, %%ymm15\n\t"
"vmovdqu 32(%[bp]), %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpmaddwd -248(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\n\t"
"vpmaddwd -120(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 8(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 136(%[ap]), %%ymm8, %%ymm13\n\t"
"vmovdqu -96(%[bp]), %%ymm8\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpmaddwd -188(%[ap]), %%ymm8, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpmaddwd -60(%[ap]), %%ymm8, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm0, %%ymm0\n\t"
"vpmaddwd 68(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm1, %%ymm1\n\t"
"vpmaddwd 196(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
"vpmaddwd -252(%[ap]), %%ymm9, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
"vpmaddwd -124(%[ap]), %%ymm9, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 4(%[ap]), %%ymm9, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 132(%[ap]), %%ymm9, %%ymm15\n\t"
"vmovdqu 96(%[bp]), %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpmaddwd -252(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\n\t"
"vpmaddwd -124(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 4(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 132(%[ap]), %%ymm8, %%ymm13\n\t"
"vmovdqu -32(%[bp]), %%ymm8\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpmaddwd -180(%[ap]), %%ymm8, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpmaddwd -52(%[ap]), %%ymm8, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm0, %%ymm0\n\t"
"vpmaddwd 76(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm1, %%ymm1\n\t"
"vpmaddwd 204(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
"vpmaddwd -244(%[ap]), %%ymm9, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
"vpmaddwd -116(%[ap]), %%ymm9, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 12(%[ap]), %%ymm9, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 140(%[ap]), %%ymm9, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpmaddwd -244(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\n\t"
"vpmaddwd -116(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 12(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 140(%[ap]), %%ymm8, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm10, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm11, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
"add $512, %[ap]\n\t"
"add $256, %[bp]\n\t"
".p2align 5\n\t"
"1:\n\t"
"prefetcht0 640(%[bp])\n\t"
"vmovdqu -128(%[bp]), %%ymm8\n\t"
"vmovdqu 0(%[bp]), %%ymm9\n\t"
"vpmaddwd -192(%[ap]), %%ymm8, %%ymm14\n\t"
"vpmaddwd -64(%[ap]), %%ymm8, %%ymm15\n\t"
"vpmaddwd 64(%[ap]), %%ymm8, %%ymm10\n\t"
"vpmaddwd 192(%[ap]), %%ymm8, %%ymm11\n\t"
"vpmaddwd -256(%[ap]), %%ymm9, %%ymm12\n\t"
"vpmaddwd -128(%[ap]), %%ymm9, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 0(%[ap]), %%ymm9, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 128(%[ap]), %%ymm9, %%ymm15\n\t"
"vmovdqu 64(%[bp]), %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpmaddwd -256(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\n\t"
"vpmaddwd -128(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 0(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 128(%[ap]), %%ymm8, %%ymm13\n\t"
"vmovdqu -64(%[bp]), %%ymm8\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpmaddwd -184(%[ap]), %%ymm8, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpmaddwd -56(%[ap]), %%ymm8, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm0, %%ymm0\n\t"
"vpmaddwd 72(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm1, %%ymm1\n\t"
"vpmaddwd 200(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
"vpmaddwd -248(%[ap]), %%ymm9, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
"vpmaddwd -120(%[ap]), %%ymm9, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 8(%[ap]), %%ymm9, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 136(%[ap]), %%ymm9, %%ymm15\n\t"
"vmovdqu 32(%[bp]), %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpmaddwd -248(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\n\t"
"vpmaddwd -120(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 8(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 136(%[ap]), %%ymm8, %%ymm13\n\t"
"vmovdqu -96(%[bp]), %%ymm8\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpmaddwd -188(%[ap]), %%ymm8, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpmaddwd -60(%[ap]), %%ymm8, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm0, %%ymm0\n\t"
"vpmaddwd 68(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm1, %%ymm1\n\t"
"vpmaddwd 196(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
"vpmaddwd -252(%[ap]), %%ymm9, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
"vpmaddwd -124(%[ap]), %%ymm9, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 4(%[ap]), %%ymm9, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 132(%[ap]), %%ymm9, %%ymm15\n\t"
"vmovdqu 96(%[bp]), %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpmaddwd -252(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\n\t"
"vpmaddwd -124(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 4(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 132(%[ap]), %%ymm8, %%ymm13\n\t"
"vmovdqu -32(%[bp]), %%ymm8\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpmaddwd -180(%[ap]), %%ymm8, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpmaddwd -52(%[ap]), %%ymm8, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm0, %%ymm0\n\t"
"vpmaddwd 76(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm1, %%ymm1\n\t"
"vpmaddwd 204(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
"vpmaddwd -244(%[ap]), %%ymm9, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
"vpmaddwd -116(%[ap]), %%ymm9, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 12(%[ap]), %%ymm9, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 140(%[ap]), %%ymm9, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpmaddwd -244(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\n\t"
"vpmaddwd -116(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 12(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 140(%[ap]), %%ymm8, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm10, %%ymm0, %%ymm0\n\t"
"vpaddd %%ymm11, %%ymm1, %%ymm1\n\t"
"vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
"vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
"add $512, %[ap]\n\t"
"add $256, %[bp]\n\t"
"sub $1, %[cnt]\n\t"
"jnz 1b\n\t"
"2:\n\t"
"vpslld $16, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm4, %%ymm0, %[q0]\n\t"
"vpslld $16, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm5, %%ymm1, %[q1]\n\t"
"vpslld $16, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm6, %%ymm2, %[q2]\n\t"
"vpslld $16, %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm7, %%ymm3, %[q3]\n\t"
: [ap] "+r" (ap), [bp] "+r" (bp), [cnt] "+r" (groups),
[q0] "=x" (q0), [q1] "=x" (q1), [q2] "=x" (q2), [q3] "=x" (q3)
: [amem] "m"(Apk), [bmem] "m"(Bpk)
: "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7", "ymm8", "ymm9", "ymm10", "ymm15"); /* ymm15 kept DECLARED on purpose -- see notes */
{ int32_t *cp = c;
if (first) {
_mm256_store_si256((__m256i *)(cp + (size_t)3 * ldc), q3);
_mm256_store_si256((__m256i *)(cp + (size_t)2 * ldc), q2);
_mm256_store_si256((__m256i *)(cp + (size_t)1 * ldc), q1);
_mm256_store_si256((__m256i *)(cp + (size_t)0 * ldc), q0);
} else {
int32_t *cpr3 = cp + (size_t)3 * ldc;
_mm256_store_si256((__m256i *)cpr3, _mm256_add_epi32(_mm256_load_si256((const __m256i *)cpr3), q3));
int32_t *cpr2 = cp + (size_t)2 * ldc;
_mm256_store_si256((__m256i *)cpr2, _mm256_add_epi32(_mm256_load_si256((const __m256i *)cpr2), q2));
int32_t *cpr1 = cp + (size_t)1 * ldc;
_mm256_store_si256((__m256i *)cpr1, _mm256_add_epi32(_mm256_load_si256((const __m256i *)cpr1), q1));
int32_t *cpr0 = cp + (size_t)0 * ldc;
_mm256_store_si256((__m256i *)cpr0, _mm256_add_epi32(_mm256_load_si256((const __m256i *)cpr0), q0));
} }
}
/* MIKRO-END */
/* ---- ld-aware packing + tiled macro-kernel (base case) ---- */
#pragma GCC push_options
#pragma GCC optimize("O2")
AVX2 static inline void packA2(int ic, int kk, int n, int lda, const T * __restrict A, int sg, const T * __restrict A2) {
const __m256i m16 = _mm256_set1_epi32(0xFFFF);
__m256i half = _mm256_set1_epi32(0x8000);
__asm__ volatile("" : "+x"(half));
/* 常量索引 -> vpermq 立即数 0x88 */
__m256i SHLO = _mm256_setr_epi8(0,1,4,5,8,9,12,13,0,1,4,5,8,9,12,13,0,1,4,5,8,9,12,13,0,1,4,5,8,9,12,13);
__m256i SHHI = _mm256_setr_epi8(2,3,6,7,10,11,14,15,2,3,6,7,10,11,14,15,2,3,6,7,10,11,14,15,2,3,6,7,10,11,14,15);
__asm__ volatile("" : "+x"(SHLO));
__asm__ volatile("" : "+x"(SHHI));
(void)n;
/* UNSWITCHED: the (A2!=0) and (sg>0) decisions are hoisted to the top so each
inner loop is branch-free and no source pointer has to live in an array.
Same stores, same order of arithmetic -> bit-identical panels. */
const size_t st = (size_t)lda;
const T *s0 = A + (size_t)ic * st + kk, *s1 = s0 + st, *s2 = s1 + st, *s3 = s2 + st;
#define PA_EMIT(RO, XV) \
do { __m256i rl = _mm256_permute4x64_epi64(_mm256_shuffle_epi8((XV), SHLO), 0x88); \
__m256i rh = _mm256_permute4x64_epi64( \
_mm256_shuffle_epi8(_mm256_add_epi32((XV), half), SHHI), 0x88); \
T *d = Apk + (size_t)(RO) * 32 + 128 * u; \
ST(d + 0, rl); _mm_store_si128((__m128i *)(d + 8), _mm256_castsi256_si128(rl)); \
ST(d + 16, rh); _mm_store_si128((__m128i *)(d + 24), _mm256_castsi256_si128(rh)); } while (0)
if (A2) {
const T *t0 = A2 + (size_t)ic * st + kk, *t1 = t0 + st, *t2 = t1 + st, *t3 = t2 + st;
if (sg > 0) {
for (int u = 0; u < KC / 8; u++) {
PA_EMIT(0, _mm256_add_epi32(LD(s0 + 8 * u), LD(t0 + 8 * u)));
PA_EMIT(1, _mm256_add_epi32(LD(s1 + 8 * u), LD(t1 + 8 * u)));
PA_EMIT(2, _mm256_add_epi32(LD(s2 + 8 * u), LD(t2 + 8 * u)));
PA_EMIT(3, _mm256_add_epi32(LD(s3 + 8 * u), LD(t3 + 8 * u)));
}
} else {
for (int u = 0; u < KC / 8; u++) {
PA_EMIT(0, _mm256_sub_epi32(LD(s0 + 8 * u), LD(t0 + 8 * u)));
PA_EMIT(1, _mm256_sub_epi32(LD(s1 + 8 * u), LD(t1 + 8 * u)));
PA_EMIT(2, _mm256_sub_epi32(LD(s2 + 8 * u), LD(t2 + 8 * u)));
PA_EMIT(3, _mm256_sub_epi32(LD(s3 + 8 * u), LD(t3 + 8 * u)));
}
}
} else {
for (int u = 0; u < KC / 8; u++) {
PA_EMIT(0, LD(s0 + 8 * u)); PA_EMIT(1, LD(s1 + 8 * u));
PA_EMIT(2, LD(s2 + 8 * u)); PA_EMIT(3, LD(s3 + 8 * u));
}
}
#undef PA_EMIT
}
#pragma GCC pop_options
#pragma GCC push_options
#pragma GCC optimize("O3")
AVX2 static inline void packB2(int jc, int kk, int n, int ldb, const T * __restrict B, int sg, const T * __restrict B2) {
const __m256i m16 = _mm256_set1_epi32(0xFFFF);
const __m256i half = _mm256_set1_epi32(0x8000);
(void)n;
for (int u = 0; u < NPAIR / 4; u++) {
const T *s = B + (size_t)(kk + 8 * u) * ldb + jc;
const T *t = B2 ? (B2 + (size_t)(kk + 8 * u) * ldb + jc) : 0;
for (int j = 0; j < NC; j += 8) {
T *d = Bpk + (size_t)(j / 8) * (NPAIR * 16 + 24) + (size_t)u * 64;
__m256i L0, L1, L2, L3, H0, H1, H2, H3;
/* [mi4av_] bit-exact, verified on a micro-test (LO/HI both match byte-for-byte):
the two byte-shuffle tables only assemble the low-16 of each dword, i.e.
L_P = (x0[i]&0xffff) | (x1[i]&0xffff)<<16 = vpblendw(x0, vpslld(x1,16), 0xAA);
H_P = high16(x+0x8000) on both rows = vpblendw(vpsrld(x0+half,16), x1+half, 0xAA).
4 vpshufb + 2 vpor -- four of them port-5 uops -- per pair become 0 port-5 uops
and 2 fewer uops per pair: -8 uops and -16 port-5 uops per 32 B output block.
The 8 stores and the 24 vpblendd are untouched. */
#define F11_PAIR(P, LO, HI) \
{ const T *sp = s + (size_t)(2 * P) * ldb + j; \
__m256i x0 = LD(sp), x1 = LD(sp + ldb); \
if (t) { const T *tp = t + (size_t)(2 * P) * ldb + j; \
__m256i y0 = LD(tp), y1 = LD(tp + ldb); \
if (sg > 0) { x0 = _mm256_add_epi32(x0, y0); x1 = _mm256_add_epi32(x1, y1); } \
else { x0 = _mm256_sub_epi32(x0, y0); x1 = _mm256_sub_epi32(x1, y1); } } \
LO = _mm256_blend_epi16(x0, _mm256_slli_epi32(x1, 16), 0xAA); \
HI = _mm256_blend_epi16(_mm256_srli_epi32(_mm256_add_epi32(x0, half), 16), \
_mm256_add_epi32(x1, half), 0xAA); }
F11_PAIR(0, L0, H0)
F11_PAIR(1, L1, H1)
F11_PAIR(2, L2, H2)
F11_PAIR(3, L3, H3)
#undef F11_PAIR
{
__m256i a=_mm256_blend_epi32(L0,L2,0xcc),b=_mm256_blend_epi32(L0,L2,0x33);
__m256i c=_mm256_blend_epi32(L1,L3,0xcc),e=_mm256_blend_epi32(L1,L3,0x33);
ST(d+0,_mm256_blend_epi32(a,c,0xaa));
ST(d+8,_mm256_blend_epi32(c,b,0xaa));
ST(d+16,_mm256_blend_epi32(b,e,0xaa));
ST(d+24,_mm256_blend_epi32(e,a,0xaa));
}
{
__m256i a=_mm256_blend_epi32(H0,H2,0xcc),b=_mm256_blend_epi32(H0,H2,0x33);
__m256i c=_mm256_blend_epi32(H1,H3,0xcc),e=_mm256_blend_epi32(H1,H3,0x33);
ST(d+32,_mm256_blend_epi32(a,c,0xaa));
ST(d+40,_mm256_blend_epi32(c,b,0xaa));
ST(d+48,_mm256_blend_epi32(b,e,0xaa));
ST(d+56,_mm256_blend_epi32(e,a,0xaa));
}
}
}
}
#pragma GCC pop_options
AVX2 static inline __attribute__((always_inline)) void block_mikro(int ic, int jc, int n, int ldc, T *C, int first) {
T *cp0 = C + (size_t)ic * ldc + jc;
for (int r = 0; r < MC; r += MR) {
const int32_t *ap = Apk + (size_t)r * APSTR;
T *cp = cp0 + (size_t)r * ldc;
const int32_t *bp = Bpk;
for (int jg = 0; jg < NC; jg += NR) {
mikro(NPAIR, ap, bp, cp + jg, ldc, first);
bp += (NPAIR * 16 + 24);
}
}
}
/* LEAF SPECIALISED: with BASE==KC==NC==128 every tiled2 leaf call has h==128 and
ldc==h==128, so both become compile-time constants -> the jc/kk/ic loops fully fold. */
AVX2 static void tiled2_128(const T *A1, const T *A2, int sgA, int lda,
const T *B1, const T *B2, int sgB, int ldb, T *C) {
const int n = 128;
const int ldc = 128;
for (int jc = 0; jc < n; jc += NC) {
for (int kk = 0; kk < n; kk += KC) {
packB2(jc, kk, n, ldb, B1, sgB, B2);
for (int ic = n - MC; ic >= 0; ic -= MC) {
packA2(ic, kk, n, lda, A1, sgA, A2);
if (ic - MC >= 0) {
const T *n1 = A1 + (size_t)(ic - MC) * lda + kk;
const T *n2 = A2 ? (A2 + (size_t)(ic - MC) * lda + kk) : 0;
for (int r = 0; r < MC; r++) {
const T *p1 = n1 + (size_t)r * lda;
const T *p2 = n2 ? (n2 + (size_t)r * lda) : 0;
__builtin_prefetch(p1, 0, 2);
if (p2) __builtin_prefetch(p2, 0, 2);
}
}
block_mikro(ic, jc, n, ldc, C, kk == 0);
}
}
}
}
/* correct for any n, just slow */
AVX2 static void generic(int n, const T *A, int lda, const T *B, int ldb, T *C, int ldc) {
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++) {
uint32_t s = 0;
for (int k = 0; k < n; k++)
s += (uint32_t)A[(size_t)i * lda + k] * (uint32_t)B[(size_t)k * ldb + j];
C[(size_t)i * ldc + j] = (T)s;
}
}
/* tiled base case: requires n % NC == 0, n % KC == 0, n >= MC */
AVX2 static void tiled(int n, const T *A, int lda, const T *B, int ldb, T *C, int ldc) {
/* C is fully written by the kk==0 pass (micro stores instead of RMWs) */
for (int jc = 0; jc < n; jc += NC) {
for (int kk = 0; kk < n; kk += KC) {
packB2(jc, kk, n, ldb, B, 1, (const T *)0);
for (int ic = 0; ic < n; ic += MC) {
packA2(ic, kk, n, lda, A, 1, (const T *)0);
block_mikro(ic, jc, n, ldc, C, kk == 0);
}
}
}
}
AVX2 static int ok_tiled(int n) {
return (n % NC) == 0 && (n % KC) == 0 && (n % MC) == 0 && (MC % MR) == 0 && (NC % NR) == 0;
}
/* ---- Strassen over Z/2^32. Exact: all ops are ring ops of Z/2^32. ---- */
#define STR_DEPTH 6
/* P = X + s*Y in ONE pass over h x h (replaces copyblk + axpy: one write instead of two) */
#pragma GCC push_options
#pragma GCC optimize("O3")
AVX2 static void addblk(int h, T * __restrict P, const T * __restrict X, int ldx, const T * __restrict Y, int ldy, int s, int zp, int sp) {
if(zp && sp){
size_t z=(size_t)h*h;
if(h>=1024){
if(s>0)for(size_t j=0;j<z;j+=8)_mm256_stream_si256((__m256i*)(P+j),_mm256_add_epi32(LD(X+j),LD(Y+j)));
else for(size_t j=0;j<z;j+=8)_mm256_stream_si256((__m256i*)(P+j),_mm256_sub_epi32(LD(X+j),LD(Y+j)));
}else{
if(s>0)for(size_t j=0;j<z;j+=8)ST(P+j,_mm256_add_epi32(LD(X+j),LD(Y+j)));
else for(size_t j=0;j<z;j+=8)ST(P+j,_mm256_sub_epi32(LD(X+j),LD(Y+j)));
}
return;
}
if(zp && h>=256){
for(int rb=0;rb<h;rb+=128)for(int cb=0;cb<h;cb+=128){
T*db=P+zoff(h,rb,cb);
const T*xb=sp?X+zoff(h,rb,cb):X+(size_t)rb*ldx+cb;
const T*yb=sp?Y+zoff(h,rb,cb):Y+(size_t)rb*ldy+cb;
const int sx=sp?128:ldx,sy=sp?128:ldy;
for(int r=0;r<128;r++){
T*d=db+(size_t)r*128;const T*x=xb+(size_t)r*sx,*y=yb+(size_t)r*sy;
__builtin_prefetch(x+sx,0,3);__builtin_prefetch(y+sy,0,3);
if(h>=1024){
if(s>0)for(int j=0;j<128;j+=8)_mm256_stream_si256((__m256i*)(d+j),_mm256_add_epi32(LD(x+j),LD(y+j)));
else for(int j=0;j<128;j+=8)_mm256_stream_si256((__m256i*)(d+j),_mm256_sub_epi32(LD(x+j),LD(y+j)));
}else{
if(s>0)for(int j=0;j<128;j+=8)ST(d+j,_mm256_add_epi32(LD(x+j),LD(y+j)));
else for(int j=0;j<128;j+=8)ST(d+j,_mm256_sub_epi32(LD(x+j),LD(y+j)));
}
}
}
return;
}
for (int i = 0; i < h; i++) {
const T *xp = X + (size_t)i * ldx, *yp = Y + (size_t)i * ldy;
__builtin_prefetch(X + (size_t)(i+1)*ldx, 0, 3); __builtin_prefetch(Y + (size_t)(i+1)*ldy, 0, 3);
if (!zp) {
T *pp = P + (size_t)i * h;
if (h >= 1024) {
if (s > 0) for (int j = 0; j < h; j += 8) _mm256_stream_si256((__m256i *)(pp + j), _mm256_add_epi32(LDA(xp + j), LDA(yp + j)));
else for (int j = 0; j < h; j += 8) _mm256_stream_si256((__m256i *)(pp + j), _mm256_sub_epi32(LDA(xp + j), LDA(yp + j)));
} else {
if (s > 0) for (int j = 0; j < h; j += 8) STA(pp + j, _mm256_add_epi32(LDA(xp + j), LDA(yp + j)));
else for (int j = 0; j < h; j += 8) STA(pp + j, _mm256_sub_epi32(LDA(xp + j), LDA(yp + j)));
}
} else if (h == 256) {
T *p0 = P + (size_t)(i >> 7) * 32768 + (size_t)(i & 127) * 128;
T *p1 = p0 + 16384;
const T *x0 = sp ? X + zoff(256, i, 0) : xp, *x1 = sp ? X + zoff(256, i, 128) : xp + 128;
const T *y0 = sp ? Y + zoff(256, i, 0) : yp, *y1 = sp ? Y + zoff(256, i, 128) : yp + 128;
if (s > 0) {
for (int j = 0; j < 128; j += 8) STA(p0 + j, _mm256_add_epi32(LDA(x0 + j), LDA(y0 + j)));
for (int j = 0; j < 128; j += 8) STA(p1 + j, _mm256_add_epi32(LDA(x1 + j), LDA(y1 + j)));
} else {
for (int j = 0; j < 128; j += 8) STA(p0 + j, _mm256_sub_epi32(LDA(x0 + j), LDA(y0 + j)));
for (int j = 0; j < 128; j += 8) STA(p1 + j, _mm256_sub_epi32(LDA(x1 + j), LDA(y1 + j)));
}
} else {
for (int c = 0; c < h; c += 128) {
T *pp = P + zoff(h, i, c);
const T *xc = sp ? X + zoff(h, i, c) : xp + c, *yc = sp ? Y + zoff(h, i, c) : yp + c;
if (h >= 1024) {
if (s > 0) for (int j = 0; j < 128; j += 8) _mm256_stream_si256((__m256i*)(pp + j), _mm256_add_epi32(LDA(xc + j), LDA(yc + j)));
else for (int j = 0; j < 128; j += 8) _mm256_stream_si256((__m256i*)(pp + j), _mm256_sub_epi32(LDA(xc + j), LDA(yc + j)));
} else {
if (s > 0) for (int j = 0; j < 128; j += 8) STA(pp + j, _mm256_add_epi32(LDA(xc + j), LDA(yc + j)));
else for (int j = 0; j < 128; j += 8) STA(pp + j, _mm256_sub_epi32(LDA(xc + j), LDA(yc + j)));
}
}
}
}
}
#pragma GCC pop_options
#pragma GCC push_options
#pragma GCC optimize("O3")
template<int SA, int SB>
AVX2 static void addblk_pair(int h, T * __restrict P,
const T * __restrict A, int lda, const T * __restrict A2,
T * __restrict Q,
const T * __restrict B, int ldb, const T * __restrict B2) {
for (int i=0;i<h;i++) {
const T *a=A+(size_t)i*lda,*a2=A2+(size_t)i*lda;
const T *b=B+(size_t)i*ldb,*b2=B2+(size_t)i*ldb;
T *p=P+(size_t)i*h,*q=Q+(size_t)i*h;
__builtin_prefetch(a+lda,0,3); __builtin_prefetch(a2+lda,0,3);
__builtin_prefetch(b+ldb,0,3); __builtin_prefetch(b2+ldb,0,3);
for(int j=0;j<h;j+=8) {
__m256i x=LD(a+j), y=LD(a2+j), z=LD(b+j), w=LD(b2+j);
STA(p+j,SA>0?_mm256_add_epi32(x,y):_mm256_sub_epi32(x,y));
STA(q+j,SB>0?_mm256_add_epi32(z,w):_mm256_sub_epi32(z,w));
}
}
}
#pragma GCC pop_options
/* dst = s1*m1 + s2*m2 , ONE pass (M blocks are contiguous, row stride h) */
/* dst = s1*m1 + s2*m2 + s3*m3 + s4*m4 , ONE pass */
#define B11 (A)
#define B12 (ap ? (A + h2) : (A + h))
#define B21 (ap ? (A + 2*h2) : (A + (size_t)h * lda))
#define B22 (ap ? (A + 3*h2) : (A + (size_t)h * lda + h))
#define D11 (B)
#define D12 (bp ? (B + h2) : (B + h))
#define D21 (bp ? (B + 2*h2) : (B + (size_t)h * ldb))
#define D22 (bp ? (B + 3*h2) : (B + (size_t)h * ldb + h))
/* ---- FUSED combine: one pass instead of four (l_sub1pct) --------------------
The shipped rec() makes FOUR passes over the seven M blocks (12 block-reads
for 4 outputs: M1..M5 are each read twice). Holding the seven row vectors in
registers lets all four quadrants be formed from 7 loads, so the combine's
issue cost per 8 outputs goes 24 -> 19 uops and its read volume 12h^2 -> 7h^2.
Same stores, same NT form, byte-exact. */
#pragma GCC push_options
#pragma GCC optimize("O3")
AVX2 static void comb_all(int h, T * __restrict C, int ldc,
const T * __restrict M1, const T * __restrict M2, const T * __restrict M3, const T * __restrict M4,
const T * __restrict M5, const T * __restrict M6, const T * __restrict M7) {
for (int i = 0; i < h; i++) {
const T *a = M1 + (size_t)i * h, *b = M2 + (size_t)i * h, *c = M3 + (size_t)i * h,
*d = M4 + (size_t)i * h, *e = M5 + (size_t)i * h, *f = M6 + (size_t)i * h,
*g = M7 + (size_t)i * h;
T *d1 = C + (size_t)i * ldc, *d2 = d1 + h,
*d3 = d1 + (size_t)h * ldc, *d4 = d3 + h;
int j = 0;
for (; j + 16 <= h; j += 16) { /* [i2y_] 同上 */
{ const int jj = j;
__m256i m1=LDA(a+jj),m2=LDA(b+jj),m3=LDA(c+jj),m4=LDA(d+jj),
m5=LDA(e+jj),m6=LDA(f+jj),m7=LDA(g+jj);
__m256i r1=LDA(a+jj+8),r2=LDA(b+jj+8),r3=LDA(c+jj+8),r4=LDA(d+jj+8),
r5=LDA(e+jj+8),r6=LDA(f+jj+8),r7=LDA(g+jj+8);
_mm256_stream_si256((__m256i *)(d1+jj), _mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(m1,m4),m5),m7));
_mm256_stream_si256((__m256i *)(d1+jj+8), _mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(r1,r4),r5),r7));
_mm256_stream_si256((__m256i *)(d3+jj), _mm256_add_epi32(m2,m4));
_mm256_stream_si256((__m256i *)(d3+jj+8), _mm256_add_epi32(r2,r4));
_mm256_stream_si256((__m256i *)(d2+jj), _mm256_add_epi32(m3,m5));
_mm256_stream_si256((__m256i *)(d2+jj+8), _mm256_add_epi32(r3,r5));
_mm256_stream_si256((__m256i *)(d4+jj), _mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(m1,m3),m2),m6));
_mm256_stream_si256((__m256i *)(d4+jj+8), _mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(r1,r3),r2),r6));
}
}
for (; j + 8 <= h; j += 8) {
__m256i m1 = LDA(a + j), m2 = LDA(b + j), m3 = LDA(c + j), m4 = LDA(d + j),
m5 = LDA(e + j), m6 = LDA(f + j), m7 = LDA(g + j);
_mm256_stream_si256((__m256i *)(d1 + j),
_mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(m1, m4), m5), m7));
_mm256_stream_si256((__m256i *)(d2 + j), _mm256_add_epi32(m3, m5));
_mm256_stream_si256((__m256i *)(d3 + j), _mm256_add_epi32(m2, m4));
_mm256_stream_si256((__m256i *)(d4 + j),
_mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(m1, m3), m2), m6));
}
/* scalar tail removed: same invariant as above. */
}
}
#pragma GCC pop_options
AVX2 static void rec(int n, const T *A, int lda, const T *B, int ldb, T *C, int ldc,
T *scr, int depth, int bp, int ap) {
if (depth >= STR_DEPTH || (n & 1) || n <= BASE || !ok_tiled(n)) {
if (ok_tiled(n)) tiled(n, A, lda, B, ldb, C, ldc);
else generic(n, A, lda, B, ldb, C, ldc);
return;
}
const int h = n / 2;
const size_t h2 = (size_t)h * h;
T *S = scr; /* 7 result blocks M1..M7 then 2 scratch S blocks */
T *M1 = S, *M2 = M1 + h2 + 4608, *M3 = M2 + h2 + 4608, *M4 = M3 + h2 + 4608,
*M5 = M4 + h2 + 4608, *M6 = M5 + h2 + 4608, *M7 = M6 + h2 + 4608;
T *P = M7 + h2 + 4608; /* operand scratch block 1 (h x h, ldp = h) */
T *Q = P + h2 + 4608; /* operand scratch block 2 */
T *sub = Q + h2 + 4608; /* recursion arena */
/* Each of the 7 products uses an operand that is A1 (+/- A2) on the A side and
B1 (+/- B2) on the B side. When the children are LEAVES (h <= BASE) the summed
operand is consumed only by the packers, so we fuse the sum INTO the packer and
never materialise P/Q: that removes 2 of every 4 accesses per summed operand
element (addblk's write + the packer's read of it, replaced by a second load of
the original block, which is exactly what addblk did anyway). */
#define LOP(a1, a2, sgA, b1, b2, sgB, Mi) \
do { \
if (h <= BASE) { \
tiled2_128((a1), (a2), (sgA), (ap ? 128 : lda), \
(b1), (b2), (sgB), (bp ? 128 : ldb), (Mi)); \
} else { \
const T *pa = (a1); int la = lda; \
const T *pb = (b1); int lb = ldb; \
if ((a2) && (b2)) { \
addblk_pair<(sgA),(sgB)>(h,P,(a1),lda,(a2),Q,(b1),ldb,(b2)); \
pa=P; la=h; pb=Q; lb=h; \
} else { \
if (a2) { addblk(h,P,(a1),lda,(a2),lda,(sgA),1,ap); pa=P; la=h; } \
if (b2) { addblk(h,Q,(b1),ldb,(b2),ldb,(sgB),ZQ,bp); pb=Q; lb=h; } \
} \
rec(h, pa, la, pb, lb, (Mi), h, sub, depth + 1, bp, ap); \
} \
} while (0)
LOP(B21, B11, -1, D11, D12, +1, M6); /* M6 = (A21-A11)(B11+B12) */
LOP(B21, B22, +1, D11, 0, +1, M2); /* M2 = (A21+A22) B11 */
LOP(B22, 0, +1, D21, D11, -1, M4); /* M4 = A22 (B21-B11) */
LOP(B11, B22, +1, D11, D22, +1, M1); /* M1 = (A11+A22)(B11+B22) */
LOP(B11, 0, +1, D12, D22, -1, M3); /* M3 = A11 (B12-B22) */
LOP(B11, B12, +1, D22, 0, +1, M5); /* M5 = (A11+A12) B22 */
LOP(B12, B22, -1, D21, D22, +1, M7); /* M7 = (A12-A22)(B21+B22) */
#undef LOP
/* C11 = M1 + M4 - M5 + M7 ; C12 = M3 + M5 ; C21 = M2 + M4 ; C22 = M1 - M2 + M3 + M6 */
comb_all(h, C, ldc, M1, M2, M3, M4, M5, M6, M7);
}
#undef B11
#undef B12
#undef B21
#undef B22
#undef D11
#undef D12
#undef D21
#undef D22
/* Update an already materialized S or T block without a separate copy.
[7] verbatim. */
#pragma GCC push_options
#pragma GCC optimize("O3")
AVX2 static void mutate_blk(int h, T *dst, const T *src, int lds, int reverse, int zp, int sp) {
if(zp && h>=1024){
if(sp){
size_t z=(size_t)h*h;
for(size_t j=0;j<z;j+=16){
__m256i d0=LD(dst+j),d1=LD(dst+j+8),s0=LD(src+j),s1=LD(src+j+8);
_mm256_stream_si256((__m256i*)(dst+j),reverse?_mm256_sub_epi32(s0,d0):_mm256_sub_epi32(d0,s0));
_mm256_stream_si256((__m256i*)(dst+j+8),reverse?_mm256_sub_epi32(s1,d1):_mm256_sub_epi32(d1,s1));
}
}else for(int rb=0;rb<h;rb+=128)for(int cb=0;cb<h;cb+=128){
T*db=dst+zoff(h,rb,cb);const T*sb=src+(size_t)rb*lds+cb;
for(int r=0;r<128;r++){
T*d=db+(size_t)r*128;const T*s=sb+(size_t)r*lds;__builtin_prefetch(s+lds,0,3);
for(int j=0;j<128;j+=16){
__m256i d0=LD(d+j),d1=LD(d+j+8),s0=LD(s+j),s1=LD(s+j+8);
_mm256_stream_si256((__m256i*)(d+j),reverse?_mm256_sub_epi32(s0,d0):_mm256_sub_epi32(d0,s0));
_mm256_stream_si256((__m256i*)(d+j+8),reverse?_mm256_sub_epi32(s1,d1):_mm256_sub_epi32(d1,s1));
}
}
}
_mm_sfence();return;
}
if(zp && sp){
size_t z=(size_t)h*h;
if(reverse)for(size_t j=0;j<z;j+=8)ST(dst+j,_mm256_sub_epi32(LD(src+j),LD(dst+j)));
else for(size_t j=0;j<z;j+=8)ST(dst+j,_mm256_sub_epi32(LD(dst+j),LD(src+j)));
return;
}
if(zp && h>=256){
for(int rb=0;rb<h;rb+=128)for(int cb=0;cb<h;cb+=128){
T*db=dst+zoff(h,rb,cb);
const T*sb=sp?src+zoff(h,rb,cb):src+(size_t)rb*lds+cb;
int ss=sp?128:lds;
for(int r=0;r<128;r++){
T*d=db+(size_t)r*128;const T*s=sb+(size_t)r*ss;__builtin_prefetch(s+ss,0,3);
if(reverse)for(int j=0;j<128;j+=8)ST(d+j,_mm256_sub_epi32(LD(s+j),LD(d+j)));
else for(int j=0;j<128;j+=8)ST(d+j,_mm256_sub_epi32(LD(d+j),LD(s+j)));
}
}
return;
}
for (int i=0;i<h;i++) {
const T *s = src + (size_t)i*lds;
__builtin_prefetch(s + lds, 0, 3);
if (!zp) {
T *d = dst + (size_t)i*h;
if (reverse) for(int j=0;j<h;j+=8) STA(d+j,_mm256_sub_epi32(LDA(s+j),LDA(d+j)));
else for(int j=0;j<h;j+=8) STA(d+j,_mm256_sub_epi32(LDA(d+j),LDA(s+j)));
} else if (h == 256) {
T *d0 = dst + zoff(256, i, 0), *d1 = d0 + 16384;
const T *s0 = sp ? src + zoff(256, i, 0) : s, *s1 = sp ? src + zoff(256, i, 128) : s + 128;
if (reverse) {
for(int j=0;j<128;j+=8) STA(d0+j,_mm256_sub_epi32(LDA(s0+j),LDA(d0+j)));
for(int j=0;j<128;j+=8) STA(d1+j,_mm256_sub_epi32(LDA(s1+j),LDA(d1+j)));
} else {
for(int j=0;j<128;j+=8) STA(d0+j,_mm256_sub_epi32(LDA(d0+j),LDA(s0+j)));
for(int j=0;j<128;j+=8) STA(d1+j,_mm256_sub_epi32(LDA(d1+j),LDA(s1+j)));
}
} else {
for (int c = 0; c < h; c += 128) {
T *d = dst + zoff(h, i, c);
const T *sc = sp ? src + zoff(h, i, c) : s + c;
if (reverse) for(int j=0;j<128;j+=8) STA(d+j,_mm256_sub_epi32(LDA(sc+j),LDA(d+j)));
else for(int j=0;j<128;j+=8) STA(d+j,_mm256_sub_epi32(LDA(d+j),LDA(sc+j)));
}
}
}
}
AVX2 static void comb_win(int h, T *C, int ldc, const T *M1, const T *M2,
const T *M3, const T *M4, const T *M5,
const T *M6, const T *M7) {
for(int i=0;i<h;i++) {
const T *a=M1+(size_t)i*h,*b=M2+(size_t)i*h,*c=M3+(size_t)i*h,
*d=M4+(size_t)i*h,*e=M5+(size_t)i*h,*f=M6+(size_t)i*h,*g=M7+(size_t)i*h;
T *o11=C+(size_t)i*ldc,*o12=o11+h,*o21=o11+(size_t)h*ldc,*o22=o21+h;
int j=0;
for(;j+16<=h;j+=16) { /* [i2y_] 整条 64 B 线一次发完（每象限两条 32 B 存背靠背） */
{ const int jj=j;
__m256i p1=LDA(a+jj),p2=LDA(b+jj),p3=LDA(c+jj),p4=LDA(d+jj),
p5=LDA(e+jj),p6=LDA(f+jj),p7=LDA(g+jj);
__m256i q1=LDA(a+jj+8),q2=LDA(b+jj+8),q3=LDA(c+jj+8),q4=LDA(d+jj+8),
q5=LDA(e+jj+8),q6=LDA(f+jj+8),q7=LDA(g+jj+8);
__m256i u2=_mm256_add_epi32(p1,p6), u3=_mm256_add_epi32(u2,p7), u4=_mm256_add_epi32(u2,p5);
__m256i v2=_mm256_add_epi32(q1,q6), v3=_mm256_add_epi32(v2,q7), v4=_mm256_add_epi32(v2,q5);
_mm256_stream_si256((__m256i*)(o11+jj),_mm256_add_epi32(p1,p2));
_mm256_stream_si256((__m256i*)(o11+jj+8),_mm256_add_epi32(q1,q2));
_mm256_stream_si256((__m256i*)(o21+jj),_mm256_sub_epi32(u3,p4));
_mm256_stream_si256((__m256i*)(o21+jj+8),_mm256_sub_epi32(v3,q4));
_mm256_stream_si256((__m256i*)(o12+jj),_mm256_add_epi32(u4,p3));
_mm256_stream_si256((__m256i*)(o12+jj+8),_mm256_add_epi32(v4,q3));
_mm256_stream_si256((__m256i*)(o22+jj),_mm256_add_epi32(u3,p5));
_mm256_stream_si256((__m256i*)(o22+jj+8),_mm256_add_epi32(v3,q5));
}
}
for(;j<h;j+=8) {
__m256i p1=LDA(a+j),p2=LDA(b+j),p3=LDA(c+j),p4=LDA(d+j),
p5=LDA(e+j),p6=LDA(f+j),p7=LDA(g+j);
__m256i u2=_mm256_add_epi32(p1,p6);
__m256i u3=_mm256_add_epi32(u2,p7);
__m256i u4=_mm256_add_epi32(u2,p5);
_mm256_stream_si256((__m256i*)(o11+j),_mm256_add_epi32(p1,p2));
_mm256_stream_si256((__m256i*)(o12+j),_mm256_add_epi32(u4,p3));
_mm256_stream_si256((__m256i*)(o21+j),_mm256_sub_epi32(u3,p4));
_mm256_stream_si256((__m256i*)(o22+j),_mm256_add_epi32(u3,p5));
}
}
}
AVX2 static void comb_win_interleaved(int h, T *C, int ldc, const T *M1, const T *M2,
const T *M3, const T *M4, const T *M5,
const T *M6, const T *M7) {
for(int i=0;i<h;i++) {
const T *a=M1+(size_t)i*(7*h),*b=M2+(size_t)i*(7*h),*c=M3+(size_t)i*(7*h),
*d=M4+(size_t)i*(7*h),*e=M5+(size_t)i*(7*h),*f=M6+(size_t)i*(7*h),*g=M7+(size_t)i*(7*h);
T *o11=C+(size_t)i*ldc,*o12=o11+h,*o21=o11+(size_t)h*ldc,*o22=o21+h;
int j=0;
for(;j+16<=h;j+=16) { /* [i2y_] 整条 64 B 线一次发完（每象限两条 32 B 存背靠背） */
{ const int jj=j;
__m256i p1=LDA(a+jj),p2=LDA(b+jj),p3=LDA(c+jj),p4=LDA(d+jj),
p5=LDA(e+jj),p6=LDA(f+jj),p7=LDA(g+jj);
__m256i q1=LDA(a+jj+8),q2=LDA(b+jj+8),q3=LDA(c+jj+8),q4=LDA(d+jj+8),
q5=LDA(e+jj+8),q6=LDA(f+jj+8),q7=LDA(g+jj+8);
__m256i u2=_mm256_add_epi32(p1,p6), u3=_mm256_add_epi32(u2,p7), u4=_mm256_add_epi32(u2,p5);
__m256i v2=_mm256_add_epi32(q1,q6), v3=_mm256_add_epi32(v2,q7), v4=_mm256_add_epi32(v2,q5);
_mm256_stream_si256((__m256i*)(o11+jj),_mm256_add_epi32(p1,p2));
_mm256_stream_si256((__m256i*)(o11+jj+8),_mm256_add_epi32(q1,q2));
_mm256_stream_si256((__m256i*)(o21+jj),_mm256_sub_epi32(u3,p4));
_mm256_stream_si256((__m256i*)(o21+jj+8),_mm256_sub_epi32(v3,q4));
_mm256_stream_si256((__m256i*)(o12+jj),_mm256_add_epi32(u4,p3));
_mm256_stream_si256((__m256i*)(o12+jj+8),_mm256_add_epi32(v4,q3));
_mm256_stream_si256((__m256i*)(o22+jj),_mm256_add_epi32(u3,p5));
_mm256_stream_si256((__m256i*)(o22+jj+8),_mm256_add_epi32(v3,q5));
}
}
for(;j<h;j+=8) {
__m256i p1=LDA(a+j),p2=LDA(b+j),p3=LDA(c+j),p4=LDA(d+j),
p5=LDA(e+j),p6=LDA(f+j),p7=LDA(g+j);
__m256i u2=_mm256_add_epi32(p1,p6);
__m256i u3=_mm256_add_epi32(u2,p7);
__m256i u4=_mm256_add_epi32(u2,p5);
_mm256_stream_si256((__m256i*)(o11+j),_mm256_add_epi32(p1,p2));
_mm256_stream_si256((__m256i*)(o12+j),_mm256_add_epi32(u4,p3));
_mm256_stream_si256((__m256i*)(o21+j),_mm256_sub_epi32(u3,p4));
_mm256_stream_si256((__m256i*)(o22+j),_mm256_add_epi32(u3,p5));
}
}
}
#pragma GCC pop_options
AVX2 static void rec_w(int n, const T *A, int lda, const T *B, int ldb,
T *C, int ldc, T *scr, int depth, int bp, int ap) {
if (n<=256 || depth>=STR_DEPTH || (n&1) || !ok_tiled(n)) {
rec(n,A,lda,B,ldb,C,ldc,scr,depth,bp,ap); return;
}
const int h=n/2; const size_t h2=(size_t)h*h;
T *M1=scr,*M2=M1+h2+4608,*M3=M2+h2+4608,*M4=M3+h2+4608,
*M5=M4+h2+4608,*M6=M5+h2+4608,*M7=M6+h2+4608;
T *P=M7+h2+4608,*Q=P+h2+4608,*sub=Q+h2+4608;
const bool interleaved = (depth == 0 && n == 4096) || (depth == 1 && n == 2048);
if (interleaved) {
M1=scr; M2=M1+h; M3=M2+h; M4=M3+h; M5=M4+h; M6=M5+h; M7=M6+h;
P=scr+7*h2+7*4608; Q=P+h2+4608; sub=Q+h2+4608;
if (g_rel && depth == 0) { /* mi2aq_: relocation is depth-0 only */
P = g_pool; Q = scr + 7*h2; sub = g_pool + h2;
}
} else if (g_rel && depth == 3) {
sub = g_deep;
}
const int mld=interleaved?7*h:h;
const size_t q2=(size_t)h*h;
const T *A11=A,*A12=ap?A+q2:A+h,*A21=ap?A+2*q2:A+(size_t)h*lda,*A22=ap?A+3*q2:A21+h;
const T *B11=B,*B12=bp?B+q2:B+h,*B21=bp?B+2*q2:B+(size_t)h*ldb,*B22=bp?B+3*q2:B21+h;
rec_w(h,A12,lda,B21,ldb,M2,mld,sub,depth+1,bp,ap); // P2
rec_w(h,A11,lda,B11,ldb,M1,mld,sub,depth+1,bp,ap); // P1
addblk(h,Q,B12,ldb,B11,ldb,-1,ZQ,bp); // T1 [mi2au_ ord-1: 第一对对调]
addblk(h,P,A21,lda,A22,lda,+1,1,ap); // S1
rec_w(h,P,h,Q,h,M5,mld,sub,depth+1,1,1); // P5
mutate_blk(h,P,A11,lda,0,1,ap); // S2 = S1-A11
mutate_blk(h,Q,B22,ldb,1,ZQ,bp); // T2 = B22-T1
rec_w(h,P,h,Q,h,M6,mld,sub,depth+1,1,1); // P6
mutate_blk(h,P,A12,lda,1,1,ap); // S4 = A12-S2
mutate_blk(h,Q,B21,ldb,0,ZQ,bp); // T4 = T2-B21
rec_w(h,A22,lda,Q,h,M4,mld,sub,depth+1,1,ap); // P4
rec_w(h,P,h,B22,ldb,M3,mld,sub,depth+1,bp,1); // P3
addblk(h,Q,B22,ldb,B12,ldb,-1,ZQ,bp); // T3
addblk(h,P,A11,lda,A21,lda,-1,1,ap); // S3
rec_w(h,P,h,Q,h,M7,mld,sub,depth+1,1,1); // P7
if(interleaved) comb_win_interleaved(h,C,ldc,M1,M2,M3,M4,M5,M6,M7);
else comb_win(h,C,ldc,M1,M2,M3,M4,M5,M6,M7);
}
AVX2 void matrix_multiply(int n, const T *A, const T *B, T *C) {
unsigned long long al = (unsigned long long)(size_t)A | (unsigned long long)(size_t)B | (unsigned long long)(size_t)C;
if (n <= BASE || !ok_tiled(n) || (al & 31)) {
if (ok_tiled(n) && !(al & 31)) tiled(n, A, n, B, n, C, n);
else generic(n, A, n, B, n, C, n);
return;
}
/* depth-first arena: 9 h^2 blocks at each level along the path */
size_t need = 0, h = (size_t)n;
while (h > (size_t)BASE) { h /= 2; need += 9 * (h * h + 4608); }
const int rel = (n == 4096);
size_t need_alloc = need;
if (rel) { /* mi2ap_: [7h^2][Q][最深一帧] */
const size_t h0 = (size_t)n / 2;
need_alloc = 7*h0*h0 + (h0*h0 + 4608) + 9*((size_t)BASE*BASE + 4608);
}
static T *scr = 0; static size_t cap = 0;
if (need_alloc > cap) { if (scr) free(scr); cap = need_alloc + 1;
void *pv = 0; if (posix_memalign(&pv, (size_t)2 << 20, cap * sizeof(T))) return; scr = (T *)pv; (void)madvise(scr,(cap*sizeof(T))&~size_t(4095),MADV_HUGEPAGE); }
g_rel = rel;
if (rel) {
const size_t h0 = (size_t)n / 2;
g_pool = C;
g_deep = scr + 7*h0*h0 + (h0*h0 + 4608);
}
rec_w(n, A, n, B, n, C, n, scr, 0, 0, 0);
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 1.695 s | 192 MB + 564 KB | Accepted | Score: 100 | 显示更多 |