/* 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;
__m256i q0, q1, q2, q3;
__asm__ volatile(
"vmovdqu 0(%[bp]), %%ymm8\n\t" "vmovdqu 128(%[bp]), %%ymm9\n\t" "vpmaddwd 0(%[ap]), %%ymm8, %%ymm0\n\t" "vpmaddwd 2240(%[ap]), %%ymm8, %%ymm1\n\t" "vpmaddwd 4480(%[ap]), %%ymm8, %%ymm2\n\t" "vpmaddwd 6720(%[ap]), %%ymm8, %%ymm3\n\t" "vpmaddwd 64(%[ap]), %%ymm8, %%ymm4\n\t" "vpmaddwd 2304(%[ap]), %%ymm8, %%ymm5\n\t" "vpmaddwd 4544(%[ap]), %%ymm8, %%ymm6\n\t" "vpmaddwd 6784(%[ap]), %%ymm8, %%ymm7\n\t" "vpmaddwd 0(%[ap]), %%ymm9, %%ymm10\n\t" "vpmaddwd 2240(%[ap]), %%ymm9, %%ymm11\n\t" "vpmaddwd 4480(%[ap]), %%ymm9, %%ymm12\n\t" "vpmaddwd 6720(%[ap]), %%ymm9, %%ymm13\n\t" "vpaddd %%ymm10, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm11, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm12, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm13, %%ymm7, %%ymm7\n\t" "vmovdqu 32(%[bp]), %%ymm8\n\t" "vmovdqu 160(%[bp]), %%ymm9\n\t" "vpmaddwd 4(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 2244(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 4484(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 6724(%[ap]), %%ymm8, %%ymm13\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" "vpmaddwd 68(%[ap]), %%ymm8, %%ymm14\n\t" "vpmaddwd 2308(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 4548(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 6788(%[ap]), %%ymm8, %%ymm11\n\t" "vpaddd %%ymm14, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm15, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm10, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm11, %%ymm7, %%ymm7\n\t" "vpmaddwd 4(%[ap]), %%ymm9, %%ymm12\n\t" "vpmaddwd 2244(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 4484(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 6724(%[ap]), %%ymm9, %%ymm15\n\t" "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm14, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t" "vmovdqu 64(%[bp]), %%ymm8\n\t" "vmovdqu 192(%[bp]), %%ymm9\n\t" "vpmaddwd 8(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 2248(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 4488(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 6728(%[ap]), %%ymm8, %%ymm13\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" "vpmaddwd 72(%[ap]), %%ymm8, %%ymm14\n\t" "vpmaddwd 2312(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 4552(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 6792(%[ap]), %%ymm8, %%ymm11\n\t" "vpaddd %%ymm14, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm15, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm10, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm11, %%ymm7, %%ymm7\n\t" "vpmaddwd 8(%[ap]), %%ymm9, %%ymm12\n\t" "vpmaddwd 2248(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 4488(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 6728(%[ap]), %%ymm9, %%ymm15\n\t" "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm14, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t" "vmovdqu 96(%[bp]), %%ymm8\n\t" "vmovdqu 224(%[bp]), %%ymm9\n\t" "vpmaddwd 12(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 2252(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 4492(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 6732(%[ap]), %%ymm8, %%ymm13\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" "vpmaddwd 76(%[ap]), %%ymm8, %%ymm14\n\t" "vpmaddwd 2316(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 4556(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 6796(%[ap]), %%ymm8, %%ymm11\n\t" "vpaddd %%ymm14, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm15, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm10, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm11, %%ymm7, %%ymm7\n\t" "vpmaddwd 12(%[ap]), %%ymm9, %%ymm12\n\t" "vpmaddwd 2252(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 4492(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 6732(%[ap]), %%ymm9, %%ymm15\n\t" "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm14, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t" "add $128, %[ap]\n\t"
"add $256, %[bp]\n\t"
".p2align 5\n\t"
"1:\n\t" "prefetcht0 768(%[bp])\n\t"
"vmovdqu 0(%[bp]), %%ymm8\n\t"
"vmovdqu 128(%[bp]), %%ymm9\n\t"
"vpmaddwd 64(%[ap]), %%ymm8, %%ymm14\n\t"
"vpmaddwd 2304(%[ap]), %%ymm8, %%ymm15\n\t"
"vpmaddwd 4544(%[ap]), %%ymm8, %%ymm10\n\t"
"vpmaddwd 6784(%[ap]), %%ymm8, %%ymm11\n\t"
"vpmaddwd 0(%[ap]), %%ymm9, %%ymm12\n\t"
"vpmaddwd 2240(%[ap]), %%ymm9, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 4480(%[ap]), %%ymm9, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 6720(%[ap]), %%ymm9, %%ymm15\n\t"
"vmovdqu 192(%[bp]), %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpmaddwd 0(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\n\t"
"vpmaddwd 2240(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 4480(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 6720(%[ap]), %%ymm8, %%ymm13\n\t"
"vmovdqu 64(%[bp]), %%ymm8\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpmaddwd 72(%[ap]), %%ymm8, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpmaddwd 2312(%[ap]), %%ymm8, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm0, %%ymm0\n\t"
"vpmaddwd 4552(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm1, %%ymm1\n\t"
"vpmaddwd 6792(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
"vpmaddwd 8(%[ap]), %%ymm9, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
"vpmaddwd 2248(%[ap]), %%ymm9, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 4488(%[ap]), %%ymm9, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 6728(%[ap]), %%ymm9, %%ymm15\n\t"
"vmovdqu 160(%[bp]), %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpmaddwd 8(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\n\t"
"vpmaddwd 2248(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 4488(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 6728(%[ap]), %%ymm8, %%ymm13\n\t"
"vmovdqu 32(%[bp]), %%ymm8\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpmaddwd 68(%[ap]), %%ymm8, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpmaddwd 2308(%[ap]), %%ymm8, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm0, %%ymm0\n\t"
"vpmaddwd 4548(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm1, %%ymm1\n\t"
"vpmaddwd 6788(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
"vpmaddwd 4(%[ap]), %%ymm9, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
"vpmaddwd 2244(%[ap]), %%ymm9, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 4484(%[ap]), %%ymm9, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 6724(%[ap]), %%ymm9, %%ymm15\n\t"
"vmovdqu 224(%[bp]), %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpmaddwd 4(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\n\t"
"vpmaddwd 2244(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 4484(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 6724(%[ap]), %%ymm8, %%ymm13\n\t"
"vmovdqu 96(%[bp]), %%ymm8\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpmaddwd 76(%[ap]), %%ymm8, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpmaddwd 2316(%[ap]), %%ymm8, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm0, %%ymm0\n\t"
"vpmaddwd 4556(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm1, %%ymm1\n\t"
"vpmaddwd 6796(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
"vpmaddwd 12(%[ap]), %%ymm9, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
"vpmaddwd 2252(%[ap]), %%ymm9, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 4492(%[ap]), %%ymm9, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 6732(%[ap]), %%ymm9, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpmaddwd 12(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\n\t"
"vpmaddwd 2252(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 4492(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 6732(%[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 $128, %[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)
:
: "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7", "ymm8", "ymm9", "ymm10", "ymm15", "memory"); /* 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) * APSTR + 32 * 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 && 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>=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.693 s | 192 MB + 564 KB | Accepted | Score: 100 | 显示更多 |