/*
References:
- saffah_codex_6s_agg2, https://duck.ac/submission/110268 : copied our
accepted mmml1k engine and all inherited citations below.
- saffah_codex_6s_agg2, work/mmml1k/direct_ac.cpp : adapted our new
direct-register output technique in this workspace.
- saffah_cc_v41_agg1, https://duck.ac/submission/109595 : inherited
Winograd/AVX2 arithmetic; no separate license notice appeared.
Approach:
- Keep six products in YMM registers after the common six-row dot product.
Before writing C1, read all C2 lines that need an add/subtract; this
preserves the original avoidance of a 4 KiB aliasing store-load pair.
Purpose:
- Official Try It comparison with the direct-register output version.
*/
/*
References:
- saffah_codex_6s_agg2, https://duck.ac/submission/109634 : copied our
accepted mmml1k Winograd/AVX2 engine with all inherited citations below.
- saffah_cc_v41_agg1, https://duck.ac/submission/109595 : inherited kernel
and packing design; no separate license notice appeared.
Approach:
- Replace the unused A-panel lookahead broadcast at byte offset 4 with a
prefetch at offset 64. All consumed A loads and arithmetic stay exact.
Purpose:
- Official Try It search for earlier L1 arrival without increasing the
microkernel's effective data work. Validate the unchanged output hash.
*/
/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/109595: copied its public accepted mmml1k engine and retained all its inherited citations. The source displays no separate license; this attributed adaptation is submitted to the same contest site.
Approach:
- For full 128 by 128 B panels, dispatch by operand count and sign before packing. The one-term path uses straight vector copies, and the two-term paths use straight vector additions or subtractions, omitting an operand-count and sign check from every 8-column group. The existing packer remains the fallback for other shapes.
Purpose:
- Experimentally measure full-panel B packing specialization on the official judge.
*/
// ===== REFERENCES (l1y_) =====
// [L1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 **#105870**
// <https://duck.ac/submission/105870>（112.513150 ms，我方当值）= 本件的**直接基座**。
// 正文逐字节取自 `problems/mmml1k/work/m1l_ncthp.cpp`，只做下面一处改动。基座自身的引用链见其头部。
// [L2] `problems/mmml1k/notes.md` `y1b_` 席 §五（**判题机 in-situ 相位账 + 剂量臂**）：
// ① packA = **10.65 Mticks / 3773 次调用**，占全局 **2.7%**；
// ② 剂量臂（每 k 迭代插 4 条 `nop`）实测 **+0.246%**（= 定量预言的 70%）⇒ **~70% 发射/前端受限**；
// ③ 反汇编：每 4 个元素 = 4×`vmovdqu` 载入 + 3×`vpaddq/vpsubq` + 1×`vmovdqu` 存储
// **+ 3 条循环不变量的分支测试**（`test %eax,%eax` / `cmp $0x3,%r12d` / `test %r13d`
// —— **nt≤2 与 nt≥3 两路被编译器合并进同一循环**）+ 地址/循环开销 ≈ **15 uops**
// ⇒ 前端地板 ≈ 3.7 拍/4 元素，**实测 7.58 拍/4 元素 = 地板的 2.0×**。
// [L3] ★ 本席**当场量到的判据**（`work/l1y_ntprobe.cpp`，n=1024 实跑打 nt 直方图）：
// **packA {nt1=1078, nt2=2695} · packB {nt1=98, nt2=245} —— `nt>=3` 一次都没有执行**
// ⇒ `packA`/`packB` 里的 `nt>=3` 通用分支在本题（n=1024）是**死代码**，而它正被编译器
// 合进热循环、且把 nt 判决留在循环体内。
// ======================
// ===== 做法（l1y_，2026-09-28）=====
// **单变量**：只改 `packA` 的**控制形态**，不动任何数值语义与内存布局 ——
// ① 把 `nt` 判决**提到 r 环之外**，按 `nt==1` / `nt==2` 分出**两条完全特化的热循环**
// （循环体内 **0 条** nt 测试、0 条 `cmp $0x3` 分支、不生成 4 载入/3 加法的通用体）；
// ② `nt>=3` 通用路**原样保留**为冷分支（`mc-` 之外），以满足一般 n 的正确性
// —— 本席**没有删它**，只是把它移出热循环。
// 依据：`[L2]③` 的反汇编显示 15 uops/4 元素里有 3 条是**循环不变量**测试；`[L3]` 证明
// 本题热路上它们**恒为同一结果**（nt 只有 1 与 2）⇒ 移除后热循环应逼近 `[L2]③` 的 3.7 拍/4 元素前端地板。
// 正确性：三条路逐元素与基座**同式同号**（nt==1 是纯拷贝；nt==2 是 `s0 ± s1`；冷路逐字保留）⇒ 输出应逐位相同。
// ===== 目的（试验性提交）=====
// 本题缺口 1.124639 ms（1.000% = 4.05 Mticks）。若 packA 从 7.58 拍/4 元素压到 ~4.5~5，即
// **−30~40% of packA = −3.2~4.3 Mticks（全局 0.8~1.1%）**，正好压在缺口上。
// 本发**只动 packA**；`packB` 的特化留作下一发（一发一变量，先看边际）。
// ==============================================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #105839 <https://duck.ac/submission/105839>
// 用途：**本文件的正文基座**（`work/m1l_nc128.cpp` = #105805 + 对手 #105791 的 NC 减半；112.541260 ms）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #105769 <https://duck.ac/submission/105769>
// 用途：本文件叠加的**第二把已验证刀**：`wg_mid` 的 `aligned(1<<21)` + `MADV_HUGEPAGE`（判题机 −25.3 µs）。
// 本条在 [1] 的世系里**缺失**（[1] 的 `wg_mid` 仍是 `aligned(64)`，无 madvise）。
// ======================
// ===== 思路 =====
// **两把各自已被判题机定价过的刀叠加**（机理互不重叠）：
// (a) [1] 的打包宽度减半：`packB` 的 `jb` 循环上界 = NC，而 `nt≤2` 时 `jb0=0` ⇒ 原 NC=256 把打包量做成 2 倍
// （对手 #105791 用它拿到 −2.35%，本账号 #105839 复刻 ✓）；
// (b) [2] 的 `wg_mid` 页表刀：4.7 MB 和式缓冲 64 B → 2 MB 对齐 + 新增构造函数对其 madvise
// （既有的 `sm_thp_init` 在它之前，名字不在作用域内）。
// 两者作用面不同（循环上界 vs 页映射）⇒ 按本族"已验证刀可叠"的惯例合发。
// 预期：112.541260（[1] 的读数）− ~25 µs（[2] 在同世系的同类读数）≈ **112.516**。
// 正确性：只改打包宽度与页映射，数值路径一字未动 ⇒ `ml7_vfy` 双引擎应与 [1] **逐位相同**。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #105805 <https://duck.ac/submission/105805>
// 用途：**本文件的正文基座**（`work/sm1_st3072_sub.cpp`，本账号最好件 115.209099 ms）。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 #105791 <https://duck.ac/submission/105791>
// 用途：本发复刻他**唯一**的增量。取证（`tools/fetch_rival.py` + 逐行 diff
// `ref2/rival_105770.cpp → ref2/rival_105791.cpp`）：他 #105770(115.220316) → #105791(**112.512637**，−2.35%)
// 的全部代码差异就是 **`#define NC 256 → 128`** 一行（他的自述："narrower reduction panels"）。
// ======================
// ===== 思路 =====
// **单变量 `#define NC 256 → 128`**（其余一字未动）。
// 机理（★ 更正我此前在 notes 里写错的"NC 是 no-op"）：`NC` 不只是 `Bpanel` 的声明尺寸 ——
// `packB` 里 `for (int jb = jb0; jb < NC; jb += 4)` 对 **nt ≤ 2**（叶的常态）取 `jb0 = 0`，
// 于是它按 **NC 列**打包，而叶实际只需要 `nc`（≤ n = CUTL = 128）列 ⇒ NC=256 时**把打包量做成了需要的 2 倍**
// （多打的那一半从不被 mikro 读走，但源读、和式、备份存储都照付）⇒ 判题机 −2.35%。
// 安全性：`nc ≤ NC` 恒成立（`jc` 循环只有 1 次，n=128 ≤ NC）⇒ 打包范围仍覆盖全部被读的列 ✓；
// 对可能出现的 `leaf(>128)`（nt>2 路径）也成立（`nc = min(n-jc, NC) ≤ NC` ✓）。
// 正确性：只改打包**宽度**（未读的列本来就不参与）⇒ `ml7_vfy` 双引擎应与基座**逐位相同**。
// ================
// ===== REFERENCES =====
// [1] 本账号（saffah_cc_v41_agg1）提交 #105779 <https://duck.ac/submission/105779>
// （115.355225 ms，我方最好）—— 本文件 = 该件正文**只改顶层 9 个 512² 平面的间距**。
// [2] duck.ac 用户 saffah_codex_6s_agg2，提交 #105740 <https://duck.ac/submission/105740>
// <https://duck.ac/submission/105740> —— "Winograd 平面间距"这条轴的来源：
// 该件整份复制我方 #105635，唯一改动就是把 **mid 级**间距 32 词提到 3072 词（自报 −0.51%）。
// 本发把同一条规律推广到**顶层**平面（9 个 512² 平面，原为 32 词），并在判题机上用
// 运行期旋钮台独立定价（见"思路"）。
// [3] 本账号 #105613 <https://duck.ac/submission/105613> 与用户 saffah_codex_6s_agg2 #105565
// <https://duck.ac/submission/105565> —— [1] 的基座（第三层由本账号 #105068 出）。
// ======================
// ===== 思路 =====
// 现状：mine = 115.355225（#105779），T = 115.420218，严支线 0.99*T+1 = 114.267016 ⇒ 缺 1.088 ms（0.94%）。
//
// ★ 唯一改动：把**顶层** `wg_buf` 里 9 个 512² 平面（M1..M7/P/Q）的间距从 **32 词（256 B = 4 条 line）**
// 提到 **3072 词（24576 B = 384 条 line）**，缓冲相应放大；纯地址改动，数值路径一字未动。
// 机理与 [2] 那条同源：**256 B 间距把 9 个平面在 L2 的 set 映射上挤在同一小组**
// （L2 = 256 KB/4-way ⇒ 1024 组；4 条 line 的相移几乎不改变组号），而 384 条 line 的相移把它们散开；
// 4K 别名相位不变（24576 = 6×4096）⇒ 与"破别名"那一族不是同一轴。
//
// ★ 定价（判题机 in-situ，**运行期旋钮 + 交错臂**，故无"整份二进制布局"混淆；mid 级固定为 3072）：
// `g_st` = 32, 3072, 32, 3072, 8192（同一二进制 5 次 n=1024）：
// st=32: 409.79 M st=3072: 408.41 / 408.46 M（两次相差 0.01%） st=8192: 409.15 M
// ⇒ **3072 比 32 快 0.33%**，比 8192 快 0.17% ⇒ 3072 是本轴极值（与 mid 级那条独立得到同一档）。
// ★ 同时复核了 mid 级在该点附近是**平台**（2048/3072/3584/4096 四档相差 ≤0.07%）⇒ 3072 可放心固定。
// 未动：mikro 汇编体、packA/packB、面板尺寸/行距、KC/MC/NC/CUTL、Strassen/Winograd 的分解与合并式、
// 预取配置。正确性：与 [1] 同进程逐元素对拍 n=512/1024 **bad=0**；旋钮台每个臂的 FNV 校验和
// 均等于基座的 9ae99c（⇒ 纯地址改动、逐位同值）。
// 附：本轴是**整程序级**的（改的是所有相位的 cache 组映射），不是"循环内调度"级 ⇒ 按本工作区
// 已验证的迁移规律（NT 写 C −0.268%、THP −0.175% 都成功迁移；循环内相位台读数则不能外推）
// **值得正式提交定价**。
// ================
/* mmml: C = A*B mod 2^64, long long, FUNCTION-style C++ linkage.
5-ALU-uop mikro (vpmuludq + vpaddq + vpmulld + vpaddd per k x row x 4 cols, no p5 uop)
plus Strassen with fused operand combination and dual-destination stores
(no extra memory). MR=6 KC=256 MC=24 NC=256 CUTL=256 U=1 BSW=0 EXP=0 */
#include <immintrin.h>
#include <stddef.h>
#include <stdint.h>
#include <sys/mman.h>
#ifndef MADV_HUGEPAGE
#define MADV_HUGEPAGE 14
#endif
/* force single-instruction 32-byte loads/stores (gcc 9.3 -O2 splits some
_mm256_loadu_si256 into vmovdqu xmm + vinserti128 = 2 extra uops) */
__attribute__((target("avx2"))) static inline __m256i LDU(const void *p) {
__m256i v; __asm__("vmovdqu %1, %0" : "=x"(v) : "m"(*(const __m256i *)p)); return v;
}
__attribute__((target("avx2"))) static inline __m256i LDA(const void *p) {
__m256i v; __asm__("vmovdqa %1, %0" : "=x"(v) : "m"(*(const __m256i *)p)); return v;
}
__attribute__((target("avx2"))) static inline void STU(void *p, __m256i v) {
__asm__("vmovdqu %0, %1" ::"x"(v), "m"(*(__m256i *)p) : "memory");
}
/* sm_nt: 32-byte non-temporal store (write-allocate bypass) */
static inline void STNT(void *p, __m256i v) {
__asm__ volatile("vmovntdq %0, %1" ::"x"(v), "m"(*(__m256i *)p) : "memory");
}
/* [1] #101958 的 mikro 函数体里用它自己的 32 字节存取宏，逐字节带回这里 */
#define LDQ(p) ({ __m256i _ldq_v; __asm__("vmovdqu %1,%0" : "=x"(_ldq_v) : "m"(*(const __m256i *)(p)) : "memory"); _ldq_v; })
#define STQ(p, v) __asm__ volatile("vmovdqu %0,%1" :: "x"(v), "m"(*(__m256i *)(p)) : "memory")
#pragma GCC optimize("O3","unroll-loops")
#pragma GCC push_options
#pragma GCC target("avx2")
typedef long long TE;
#define MR 6
#define KC 128
#define MC 12
#define NC 128
#define UU 1
#define CUTL 128
static TE Apanel[(size_t)MC * KC * 1 + 64] __attribute__((aligned(4096)));
static TE Bpanel[(size_t)2 * (NC / 4 + 1) * KC * 8 + 64] __attribute__((aligned(4096)));
static __m256i AC[12] __attribute__((aligned(64)));
/* op: 0 = store, 1 = add, 2 = sub; c2 may be 0 (single destination) */
static inline void mikro(long cnt, const long long *ap, const long long *bp,
long long *c1, int ldc1, int op1,
long long *c2, int ldc2, int op2, int rows) {
long idx = -(cnt << 3);
long long *apE = (long long *)ap + cnt;
long long *bpE = (long long *)bp + (cnt << 2);
if (rows == 6) {
register __m256i R0 __asm__("ymm1");
register __m256i R1 __asm__("ymm3");
register __m256i R2 __asm__("ymm5");
register __m256i R3 __asm__("ymm7");
register __m256i R4 __asm__("ymm9");
register __m256i R5 __asm__("ymm11");
__asm__ volatile(
"vpxor %%ymm0, %%ymm0, %%ymm0\n\t"
"vmovdqa %%ymm0, %%ymm1\n\t"
"vmovdqa %%ymm0, %%ymm2\n\t"
"vmovdqa %%ymm0, %%ymm3\n\t"
"vmovdqa %%ymm0, %%ymm4\n\t"
"vmovdqa %%ymm0, %%ymm5\n\t"
"vmovdqa %%ymm0, %%ymm6\n\t"
"vmovdqa %%ymm0, %%ymm7\n\t"
"vmovdqa %%ymm0, %%ymm8\n\t"
"vmovdqa %%ymm0, %%ymm9\n\t"
"vmovdqa %%ymm0, %%ymm10\n\t"
"vmovdqa %%ymm0, %%ymm11\n\t"
"test %[idx], %[idx]\n\t"
"jz 2f\n\t"
"vmovdqu (%[bp],%[idx],4), %%ymm12\n\t"
"1:\n\t"
"vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
"vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
"vbroadcastsd 1024(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
"vbroadcastsd 2048(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
"vbroadcastsd 3072(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
"vbroadcastsd 4096(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm6, %%ymm6\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm9, %%ymm9\n\t"
"vbroadcastsd 5120(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
"prefetcht0 64(%[ap],%[idx],1)\n\t"
"vpaddq %%ymm15, %%ymm10, %%ymm10\n\t"
"vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
"add $8, %[idx]\n\t"
"jnz 1b\n\t"
"2:\n\t"
"vpsrlq $32, %%ymm1, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
"vpsllq $32, %%ymm1, %%ymm1\n\t"
"vpaddq %%ymm0, %%ymm1, %%ymm1\n\t"
"vpsrlq $32, %%ymm3, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
"vpsllq $32, %%ymm3, %%ymm3\n\t"
"vpaddq %%ymm2, %%ymm3, %%ymm3\n\t"
"vpsrlq $32, %%ymm5, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpsllq $32, %%ymm5, %%ymm5\n\t"
"vpaddq %%ymm4, %%ymm5, %%ymm5\n\t"
"vpsrlq $32, %%ymm7, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpsllq $32, %%ymm7, %%ymm7\n\t"
"vpaddq %%ymm6, %%ymm7, %%ymm7\n\t"
"vpsrlq $32, %%ymm9, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpsllq $32, %%ymm9, %%ymm9\n\t"
"vpaddq %%ymm8, %%ymm9, %%ymm9\n\t"
"vpsrlq $32, %%ymm11, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm11, %%ymm11\n\t"
"vpsllq $32, %%ymm11, %%ymm11\n\t"
"vpaddq %%ymm10, %%ymm11, %%ymm11\n\t"
""
: [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx),
[R0] "=x"(R0), [R1] "=x"(R1), [R2] "=x"(R2), [R3] "=x"(R3), [R4] "=x"(R4), [R5] "=x"(R5)
:
: "ymm0", "ymm2", "ymm4", "ymm6",
"ymm8", "ymm10", "ymm12", "ymm13", "ymm14", "ymm15",
"cc", "memory");
__m256i D0,D1,D2,D3,D4,D5;
if(c2 && op2 != 0){
D0=LDQ((const __m256i*)(c2+(size_t)0*ldc2));
D1=LDQ((const __m256i*)(c2+(size_t)1*ldc2));
D2=LDQ((const __m256i*)(c2+(size_t)2*ldc2));
D3=LDQ((const __m256i*)(c2+(size_t)3*ldc2));
D4=LDQ((const __m256i*)(c2+(size_t)4*ldc2));
D5=LDQ((const __m256i*)(c2+(size_t)5*ldc2));
}
if (op1 == 0) {
STQ((__m256i*)(c1 + (size_t)0*ldc1), R0);
STQ((__m256i*)(c1 + (size_t)1*ldc1), R1);
STQ((__m256i*)(c1 + (size_t)2*ldc1), R2);
STQ((__m256i*)(c1 + (size_t)3*ldc1), R3);
STQ((__m256i*)(c1 + (size_t)4*ldc1), R4);
STQ((__m256i*)(c1 + (size_t)5*ldc1), R5);
} else {
__m256i E1_0=LDQ((const __m256i*)(c1+(size_t)0*ldc1));
__m256i E1_1=LDQ((const __m256i*)(c1+(size_t)1*ldc1));
__m256i E1_2=LDQ((const __m256i*)(c1+(size_t)2*ldc1));
__m256i E1_3=LDQ((const __m256i*)(c1+(size_t)3*ldc1));
__m256i E1_4=LDQ((const __m256i*)(c1+(size_t)4*ldc1));
__m256i E1_5=LDQ((const __m256i*)(c1+(size_t)5*ldc1));
if (op1 == 1) {
STQ((__m256i*)(c1+(size_t)0*ldc1), _mm256_add_epi64(E1_0,R0));
STQ((__m256i*)(c1+(size_t)1*ldc1), _mm256_add_epi64(E1_1,R1));
STQ((__m256i*)(c1+(size_t)2*ldc1), _mm256_add_epi64(E1_2,R2));
STQ((__m256i*)(c1+(size_t)3*ldc1), _mm256_add_epi64(E1_3,R3));
STQ((__m256i*)(c1+(size_t)4*ldc1), _mm256_add_epi64(E1_4,R4));
STQ((__m256i*)(c1+(size_t)5*ldc1), _mm256_add_epi64(E1_5,R5));
} else {
STQ((__m256i*)(c1+(size_t)0*ldc1), _mm256_sub_epi64(E1_0,R0));
STQ((__m256i*)(c1+(size_t)1*ldc1), _mm256_sub_epi64(E1_1,R1));
STQ((__m256i*)(c1+(size_t)2*ldc1), _mm256_sub_epi64(E1_2,R2));
STQ((__m256i*)(c1+(size_t)3*ldc1), _mm256_sub_epi64(E1_3,R3));
STQ((__m256i*)(c1+(size_t)4*ldc1), _mm256_sub_epi64(E1_4,R4));
STQ((__m256i*)(c1+(size_t)5*ldc1), _mm256_sub_epi64(E1_5,R5));
}
}
if(c2){
if(op2==0){
STQ((__m256i*)(c2+(size_t)0*ldc2),R0);
STQ((__m256i*)(c2+(size_t)1*ldc2),R1);
STQ((__m256i*)(c2+(size_t)2*ldc2),R2);
STQ((__m256i*)(c2+(size_t)3*ldc2),R3);
STQ((__m256i*)(c2+(size_t)4*ldc2),R4);
STQ((__m256i*)(c2+(size_t)5*ldc2),R5);
}else if(op2==1){
STQ((__m256i*)(c2+(size_t)0*ldc2),_mm256_add_epi64(D0,R0));
STQ((__m256i*)(c2+(size_t)1*ldc2),_mm256_add_epi64(D1,R1));
STQ((__m256i*)(c2+(size_t)2*ldc2),_mm256_add_epi64(D2,R2));
STQ((__m256i*)(c2+(size_t)3*ldc2),_mm256_add_epi64(D3,R3));
STQ((__m256i*)(c2+(size_t)4*ldc2),_mm256_add_epi64(D4,R4));
STQ((__m256i*)(c2+(size_t)5*ldc2),_mm256_add_epi64(D5,R5));
}else{
STQ((__m256i*)(c2+(size_t)0*ldc2),_mm256_sub_epi64(D0,R0));
STQ((__m256i*)(c2+(size_t)1*ldc2),_mm256_sub_epi64(D1,R1));
STQ((__m256i*)(c2+(size_t)2*ldc2),_mm256_sub_epi64(D2,R2));
STQ((__m256i*)(c2+(size_t)3*ldc2),_mm256_sub_epi64(D3,R3));
STQ((__m256i*)(c2+(size_t)4*ldc2),_mm256_sub_epi64(D4,R4));
STQ((__m256i*)(c2+(size_t)5*ldc2),_mm256_sub_epi64(D5,R5));
}
}
return;
} else
if (rows == 4) {
__asm__ volatile(
"vpxor %%ymm0, %%ymm0, %%ymm0\n\t"
"vmovdqa %%ymm0, %%ymm1\n\t"
"vmovdqa %%ymm0, %%ymm2\n\t"
"vmovdqa %%ymm0, %%ymm3\n\t"
"vmovdqa %%ymm0, %%ymm4\n\t"
"vmovdqa %%ymm0, %%ymm5\n\t"
"vmovdqa %%ymm0, %%ymm6\n\t"
"vmovdqa %%ymm0, %%ymm7\n\t"
"test %[idx], %[idx]\n\t"
"jz 2f\n\t"
"vmovdqu (%[bp],%[idx],4), %%ymm12\n\t"
"1:\n\t"
"vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
"vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
"vbroadcastsd 1024(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
"vbroadcastsd 2048(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
"vbroadcastsd 3072(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
"prefetcht0 64(%[ap],%[idx],1)\n\t"
"vpaddq %%ymm15, %%ymm6, %%ymm6\n\t"
"vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
"add $8, %[idx]\n\t"
"jnz 1b\n\t"
"2:\n\t"
"vpsrlq $32, %%ymm1, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
"vpsllq $32, %%ymm1, %%ymm1\n\t"
"vpaddq %%ymm0, %%ymm1, %%ymm1\n\t"
"vmovdqa %%ymm1, 0(%[ac])\n\t"
"vpsrlq $32, %%ymm3, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
"vpsllq $32, %%ymm3, %%ymm3\n\t"
"vpaddq %%ymm2, %%ymm3, %%ymm3\n\t"
"vmovdqa %%ymm3, 32(%[ac])\n\t"
"vpsrlq $32, %%ymm5, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpsllq $32, %%ymm5, %%ymm5\n\t"
"vpaddq %%ymm4, %%ymm5, %%ymm5\n\t"
"vmovdqa %%ymm5, 64(%[ac])\n\t"
"vpsrlq $32, %%ymm7, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpsllq $32, %%ymm7, %%ymm7\n\t"
"vpaddq %%ymm6, %%ymm7, %%ymm7\n\t"
"vmovdqa %%ymm7, 96(%[ac])\n\t"
""
: [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx)
: [ac] "r"(&AC[0])
: "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7",
"ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15",
"cc", "memory");
} else
if (rows == 2) {
__asm__ volatile(
"vpxor %%ymm0, %%ymm0, %%ymm0\n\t"
"vmovdqa %%ymm0, %%ymm1\n\t"
"vmovdqa %%ymm0, %%ymm2\n\t"
"vmovdqa %%ymm0, %%ymm3\n\t"
"test %[idx], %[idx]\n\t"
"jz 2f\n\t"
"vmovdqu (%[bp],%[idx],4), %%ymm12\n\t"
"1:\n\t"
"vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
"vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
"vbroadcastsd 1024(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
"prefetcht0 64(%[ap],%[idx],1)\n\t"
"vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
"vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
"add $8, %[idx]\n\t"
"jnz 1b\n\t"
"2:\n\t"
"vpsrlq $32, %%ymm1, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
"vpsllq $32, %%ymm1, %%ymm1\n\t"
"vpaddq %%ymm0, %%ymm1, %%ymm1\n\t"
"vmovdqa %%ymm1, 0(%[ac])\n\t"
"vpsrlq $32, %%ymm3, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
"vpsllq $32, %%ymm3, %%ymm3\n\t"
"vpaddq %%ymm2, %%ymm3, %%ymm3\n\t"
"vmovdqa %%ymm3, 32(%[ac])\n\t"
""
: [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx)
: [ac] "r"(&AC[0])
: "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7",
"ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15",
"cc", "memory");
} else
{
__asm__ volatile(
"vpxor %%ymm0, %%ymm0, %%ymm0\n\t"
"vmovdqa %%ymm0, %%ymm1\n\t"
"vmovdqa %%ymm0, %%ymm2\n\t"
"vmovdqa %%ymm0, %%ymm3\n\t"
"vmovdqa %%ymm0, %%ymm4\n\t"
"vmovdqa %%ymm0, %%ymm5\n\t"
"vmovdqa %%ymm0, %%ymm6\n\t"
"vmovdqa %%ymm0, %%ymm7\n\t"
"vmovdqa %%ymm0, %%ymm8\n\t"
"vmovdqa %%ymm0, %%ymm9\n\t"
"vmovdqa %%ymm0, %%ymm10\n\t"
"vmovdqa %%ymm0, %%ymm11\n\t"
"test %[idx], %[idx]\n\t"
"jz 2f\n\t"
"vmovdqu (%[bp],%[idx],4), %%ymm12\n\t"
"1:\n\t"
"vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
"vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
"vbroadcastsd 1024(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
"vbroadcastsd 2048(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
"vbroadcastsd 3072(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm4, %%ymm4\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
"vbroadcastsd 4096(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm6, %%ymm6\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm9, %%ymm9\n\t"
"vbroadcastsd 5120(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %%ymm8, %%ymm8\n\t"
"vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
"vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
"prefetcht0 64(%[ap],%[idx],1)\n\t"
"vpaddq %%ymm15, %%ymm10, %%ymm10\n\t"
"vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
"add $8, %[idx]\n\t"
"jnz 1b\n\t"
"2:\n\t"
"vpsrlq $32, %%ymm1, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
"vpsllq $32, %%ymm1, %%ymm1\n\t"
"vpaddq %%ymm0, %%ymm1, %%ymm1\n\t"
"vmovdqa %%ymm1, 0(%[ac])\n\t"
"vpsrlq $32, %%ymm3, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
"vpsllq $32, %%ymm3, %%ymm3\n\t"
"vpaddq %%ymm2, %%ymm3, %%ymm3\n\t"
"vmovdqa %%ymm3, 32(%[ac])\n\t"
"vpsrlq $32, %%ymm5, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpsllq $32, %%ymm5, %%ymm5\n\t"
"vpaddq %%ymm4, %%ymm5, %%ymm5\n\t"
"vmovdqa %%ymm5, 64(%[ac])\n\t"
"vpsrlq $32, %%ymm7, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpsllq $32, %%ymm7, %%ymm7\n\t"
"vpaddq %%ymm6, %%ymm7, %%ymm7\n\t"
"vmovdqa %%ymm7, 96(%[ac])\n\t"
"vpsrlq $32, %%ymm9, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
"vpsllq $32, %%ymm9, %%ymm9\n\t"
"vpaddq %%ymm8, %%ymm9, %%ymm9\n\t"
"vmovdqa %%ymm9, 128(%[ac])\n\t"
"vpsrlq $32, %%ymm11, %%ymm15\n\t"
"vpaddd %%ymm15, %%ymm11, %%ymm11\n\t"
"vpsllq $32, %%ymm11, %%ymm11\n\t"
"vpaddq %%ymm10, %%ymm11, %%ymm11\n\t"
"vmovdqa %%ymm11, 160(%[ac])\n\t"
""
: [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx)
: [ac] "r"(&AC[0])
: "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7",
"ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15",
"cc", "memory");
}
if (rows > 0) {
/* batch: all C reads first, then all stores. c1 and c2 are 4 MB apart for
the dual-destination M products, i.e. congruent mod 4096, so an
interleaved load/store pairs them into the 4K-aliasing stall. */
__m256i E1[6], E2[6];
if (op1) for (int r = 0; r < rows; r++) E1[r] = LDQ((const __m256i *)(c1 + (size_t)r * ldc1));
if (c2 && op2) for (int r = 0; r < rows; r++) E2[r] = LDQ((const __m256i *)(c2 + (size_t)r * ldc2));
for (int r = 0; r < rows; r++) {
__m256i res = AC[r];
long long *cp = c1 + (size_t)r * ldc1;
__m256i v = (op1 == 0) ? res : (op1 == 1 ? _mm256_add_epi64(E1[r], res) : _mm256_sub_epi64(E1[r], res));
STQ((__m256i *)cp, v);
if (c2) { cp = c2 + (size_t)r * ldc2;
v = (op2 == 0) ? res : (op2 == 1 ? _mm256_add_epi64(E2[r], res) : _mm256_sub_epi64(E2[r], res));
STQ((__m256i *)cp, v); }
}
}
}
struct Opd { /* value(i,j) = sum_t s[t]*p[(r[t]+i)*ld + c[t]+j] (nt <= 4) */
const long long *p; int ld; int nt; int r[4], c[4], s[4];
};
static inline Opd oq(const Opd *o, int qi, int qj, int h) {
Opd d = *o;
for (int t = 0; t < o->nt; t++) { d.r[t] += qi * h; d.c[t] += qj * h; }
return d;
}
static inline Opd oadd(const Opd *a, const Opd *b, int sgn) {
Opd d = *a;
for (int t = 0; t < b->nt; t++) { d.r[d.nt] = b->r[t]; d.c[d.nt] = b->c[t]; d.s[d.nt] = b->s[t] < 0 ? -sgn : sgn; d.nt++; }
return d;
}
static inline Opd o1(const long long *p, int ld) {
Opd d; d.p = p; d.ld = ld; d.nt = 1; d.r[0] = 0; d.c[0] = 0; d.s[0] = 1; return d;
}
/* A panel: rows [ic, ic+mc), k in [pc, pc+kc) of the operand */
static void packA(const Opd *o, int ic, int mc, int pc, int kc) {
long long *p = Apanel;
const int nt = o->nt;
const long long *P = o->p; const int ld = o->ld;
/* [l1y_] nt 特化：把 nt 判决提到 r 环之外，热循环体内 0 条 nt 测试。
nt>=3 的通用路在下面**原样保留**（冷分支），以满足一般 n 的正确性。 */
if (nt == 1) {
const long long *s0b = P + (size_t)o->c[0] + pc;
const long r0 = o->r[0];
for (int r = 0; r < MC; r++) {
int k = 0;
if (r < mc) {
const long long *s0 = s0b + (size_t)(r0 + ic + r) * ld;
for (; k + 4 <= kc; k += 4) { STU(p, LDU(s0 + k)); p += 4; }
for (; k < kc; k++) *p++ = s0[k];
}
for (; k < KC; k++) *p++ = 0;
}
return;
}
if (nt == 2) {
const long long *s0b = P + (size_t)o->c[0] + pc;
const long long *s1b = P + (size_t)o->c[1] + pc;
const long r0 = o->r[0], r1 = o->r[1];
if (o->s[1] > 0) {
for (int r = 0; r < MC; r++) {
int k = 0;
if (r < mc) {
const long long *s0 = s0b + (size_t)(r0 + ic + r) * ld;
const long long *s1 = s1b + (size_t)(r1 + ic + r) * ld;
for (; k + 4 <= kc; k += 4)
{ STU(p, _mm256_add_epi64(LDU(s0 + k), LDU(s1 + k))); p += 4; }
for (; k < kc; k++) *p++ = s0[k] + s1[k];
}
for (; k < KC; k++) *p++ = 0;
}
} else {
for (int r = 0; r < MC; r++) {
int k = 0;
if (r < mc) {
const long long *s0 = s0b + (size_t)(r0 + ic + r) * ld;
const long long *s1 = s1b + (size_t)(r1 + ic + r) * ld;
for (; k + 4 <= kc; k += 4)
{ STU(p, _mm256_sub_epi64(LDU(s0 + k), LDU(s1 + k))); p += 4; }
for (; k < kc; k++) *p++ = s0[k] - s1[k];
}
for (; k < KC; k++) *p++ = 0;
}
}
return;
}
for (int r = 0; r < MC; r++) {
int k = 0;
if (r < mc) {
const long long *s0 = P + (size_t)(o->r[0] + ic + r) * ld + o->c[0] + pc;
if (nt == 1 || nt == 2) {
const long long *s1 = (nt == 2) ? P + (size_t)(o->r[1] + ic + r) * ld + o->c[1] + pc : 0;
int sgn1 = (nt == 2) ? o->s[1] : 1;
for (; k + 4 <= kc; k += 4) {
__m256i v = LDU(s0 + k);
if (nt == 2) {
__m256i w = LDU(s1 + k);
v = (sgn1 > 0) ? _mm256_add_epi64(v, w) : _mm256_sub_epi64(v, w);
}
STU(p, v);
p += 4;
}
for (; k < kc; k++) { long long v = s0[k];
if (nt == 2) v = (sgn1 > 0) ? v + s1[k] : v - s1[k];
*p++ = (v); }
} else {
/* nt >= 3 (mmxR 2-layer composition): same 4-wide vector path as nt == 2,
just more terms. Signs are per-term constants of the Opd. */
const long long *sa = P + (size_t)(o->r[1] + ic + r) * ld + o->c[1] + pc;
const long long *sb = P + (size_t)(o->r[2] + ic + r) * ld + o->c[2] + pc;
const long long *sc = (nt >= 4) ? P + (size_t)(o->r[3] + ic + r) * ld + o->c[3] + pc : 0;
const int ga = o->s[1] > 0, gb = o->s[2] > 0, gc = (nt >= 4) ? (o->s[3] > 0) : 1;
for (; k + 4 <= kc; k += 4) {
__m256i v = LDU(s0 + k);
__m256i wa = LDU(sa + k);
v = ga ? _mm256_add_epi64(v, wa) : _mm256_sub_epi64(v, wa);
__m256i wb = LDU(sb + k);
v = gb ? _mm256_add_epi64(v, wb) : _mm256_sub_epi64(v, wb);
if (nt >= 4) {
__m256i wc = LDU(sc + k);
v = gc ? _mm256_add_epi64(v, wc) : _mm256_sub_epi64(v, wc);
}
STU(p, v);
p += 4;
}
for (; k < kc; k++) {
long long v = 0;
for (int t = 0; t < nt; t++)
v += o->s[t] < 0 ? -P[(size_t)(o->r[t] + ic + r) * ld + o->c[t] + pc + k]
: P[(size_t)(o->r[t] + ic + r) * ld + o->c[t] + pc + k];
*p++ = (v);
}
}
}
for (; k < KC; k++) { *p++ = 0; }
}
}
static void packB(const Opd *o, int jc, int nc, int pc, int kc) {
const int nt = o->nt;
const long long *P = o->p; const int ld = o->ld;
if (nc == NC && kc == KC && nt <= 2) {
const long long *s0 = P + (size_t)o->r[0] * ld + o->c[0] + jc;
if (nt == 1) {
for (int k = 0; k < KC; k += 2) {
const long long *r0 = s0 + (size_t)(pc + k) * ld;
const long long *r0b = r0 + ld;
long long *qb = Bpanel + (size_t)k * 4;
for (int g = 0; g < NC / 8; ++g) {
long long *q = qb + (size_t)g * 2 * KC * 4;
STU(q, LDU(r0 + g * 8));
STU(q + 4, LDU(r0b + g * 8));
STU(q + (size_t)KC * 4, LDU(r0 + g * 8 + 4));
STU(q + (size_t)KC * 4 + 4, LDU(r0b + g * 8 + 4));
}
}
} else {
const long long *s1 = P + (size_t)o->r[1] * ld + o->c[1] + jc;
if (o->s[1] > 0) {
for (int k = 0; k < KC; k += 2) {
const long long *r0 = s0 + (size_t)(pc + k) * ld;
const long long *r1 = s1 + (size_t)(pc + k) * ld;
long long *qb = Bpanel + (size_t)k * 4;
for (int g = 0; g < NC / 8; ++g) {
long long *q = qb + (size_t)g * 2 * KC * 4;
STU(q, _mm256_add_epi64(LDU(r0 + g*8), LDU(r1 + g*8)));
STU(q + 4, _mm256_add_epi64(LDU(r0 + ld + g*8), LDU(r1 + ld + g*8)));
STU(q + (size_t)KC*4, _mm256_add_epi64(LDU(r0 + g*8+4), LDU(r1 + g*8+4)));
STU(q + (size_t)KC*4 + 4, _mm256_add_epi64(LDU(r0 + ld + g*8+4), LDU(r1 + ld + g*8+4)));
}
}
} else {
for (int k = 0; k < KC; k += 2) {
const long long *r0 = s0 + (size_t)(pc + k) * ld;
const long long *r1 = s1 + (size_t)(pc + k) * ld;
long long *qb = Bpanel + (size_t)k * 4;
for (int g = 0; g < NC / 8; ++g) {
long long *q = qb + (size_t)g * 2 * KC * 4;
STU(q, _mm256_sub_epi64(LDU(r0 + g*8), LDU(r1 + g*8)));
STU(q + 4, _mm256_sub_epi64(LDU(r0 + ld + g*8), LDU(r1 + ld + g*8)));
STU(q + (size_t)KC*4, _mm256_sub_epi64(LDU(r0 + g*8+4), LDU(r1 + g*8+4)));
STU(q + (size_t)KC*4 + 4, _mm256_sub_epi64(LDU(r0 + ld + g*8+4), LDU(r1 + ld + g*8+4)));
}
}
}
}
return;
}
/* 8-column fast path: read a full 64-byte cache line of B per row per pass.
The 4-column path used only 32 B of each 64 B line while the row stride is
ld*8 bytes, so each line was fetched twice. Output layout is unchanged. */
int jb0 = 0;
if (nt >= 3) {
/* mmxR 2-layer composition: 3-4 term sources, same layout as the nt == 2 path */
const long long *s0 = P + (size_t)o->r[0] * ld + o->c[0] + jc;
const long long *s1 = P + (size_t)o->r[1] * ld + o->c[1] + jc;
const long long *s2 = P + (size_t)o->r[2] * ld + o->c[2] + jc;
const long long *s3 = (nt >= 4) ? P + (size_t)o->r[3] * ld + o->c[3] + jc : 0;
const int g1 = o->s[1] > 0, g2 = o->s[2] > 0, g3 = (nt >= 4) ? (o->s[3] > 0) : 1;
int ngrp = 0;
for (int jb = 0; jb + 8 <= nc; jb += 8) ngrp++;
int k = 0;
for (; k + 2 <= kc; k += 2) {
const long long *r0 = s0 + (size_t)(pc + k) * ld, *r0b = s0 + (size_t)(pc + k + 1) * ld;
const long long *r1 = s1 + (size_t)(pc + k) * ld, *r1b = s1 + (size_t)(pc + k + 1) * ld;
const long long *r2 = s2 + (size_t)(pc + k) * ld, *r2b = s2 + (size_t)(pc + k + 1) * ld;
const long long *r3 = (nt >= 4) ? s3 + (size_t)(pc + k) * ld : 0;
const long long *r3b = (nt >= 4) ? s3 + (size_t)(pc + k + 1) * ld : 0;
long long *qb = Bpanel + (size_t)k * 4;
for (int g = 0; g < ngrp; g++) {
long long *q0 = qb + (size_t)g * 2 * KC * 4;
__m256i a0 = LDU(r0 + g * 8), a1 = LDU(r0 + g * 8 + 4);
__m256i c0 = LDU(r0b + g * 8), c1 = LDU(r0b + g * 8 + 4);
{ __m256i b0 = LDU(r1 + g * 8), b1 = LDU(r1 + g * 8 + 4);
__m256i d0 = LDU(r1b + g * 8), d1 = LDU(r1b + g * 8 + 4);
if (g1) { a0=_mm256_add_epi64(a0,b0); a1=_mm256_add_epi64(a1,b1);
c0=_mm256_add_epi64(c0,d0); c1=_mm256_add_epi64(c1,d1); }
else { a0=_mm256_sub_epi64(a0,b0); a1=_mm256_sub_epi64(a1,b1);
c0=_mm256_sub_epi64(c0,d0); c1=_mm256_sub_epi64(c1,d1); } }
{ __m256i b0 = LDU(r2 + g * 8), b1 = LDU(r2 + g * 8 + 4);
__m256i d0 = LDU(r2b + g * 8), d1 = LDU(r2b + g * 8 + 4);
if (g2) { a0=_mm256_add_epi64(a0,b0); a1=_mm256_add_epi64(a1,b1);
c0=_mm256_add_epi64(c0,d0); c1=_mm256_add_epi64(c1,d1); }
else { a0=_mm256_sub_epi64(a0,b0); a1=_mm256_sub_epi64(a1,b1);
c0=_mm256_sub_epi64(c0,d0); c1=_mm256_sub_epi64(c1,d1); } }
if (nt >= 4) {
__m256i b0 = LDU(r3 + g * 8), b1 = LDU(r3 + g * 8 + 4);
__m256i d0 = LDU(r3b + g * 8), d1 = LDU(r3b + g * 8 + 4);
if (g3) { a0=_mm256_add_epi64(a0,b0); a1=_mm256_add_epi64(a1,b1);
c0=_mm256_add_epi64(c0,d0); c1=_mm256_add_epi64(c1,d1); }
else { a0=_mm256_sub_epi64(a0,b0); a1=_mm256_sub_epi64(a1,b1);
c0=_mm256_sub_epi64(c0,d0); c1=_mm256_sub_epi64(c1,d1); }
}
STU(q0, a0);
STU(q0 + 4, c0);
STU(q0 + (size_t)KC * 4, a1);
STU(q0 + (size_t)KC * 4 + 4, c1);
}
}
for (; k < kc; k++) {
const long long *r0 = s0 + (size_t)(pc + k) * ld;
const long long *r1 = s1 + (size_t)(pc + k) * ld;
const long long *r2 = s2 + (size_t)(pc + k) * ld;
const long long *r3 = (nt >= 4) ? s3 + (size_t)(pc + k) * ld : 0;
long long *qb = Bpanel + (size_t)k * 4;
for (int g = 0; g < ngrp; g++) {
const int jb = g * 8;
long long *q0 = qb + (size_t)(jb / 4) * KC * 4;
__m256i a0 = LDU(r0 + jb), a1 = LDU(r0 + jb + 4);
{ __m256i b0 = LDU(r1 + jb), b1 = LDU(r1 + jb + 4);
if (g1) { a0=_mm256_add_epi64(a0,b0); a1=_mm256_add_epi64(a1,b1); }
else { a0=_mm256_sub_epi64(a0,b0); a1=_mm256_sub_epi64(a1,b1); } }
{ __m256i b0 = LDU(r2 + jb), b1 = LDU(r2 + jb + 4);
if (g2) { a0=_mm256_add_epi64(a0,b0); a1=_mm256_add_epi64(a1,b1); }
else { a0=_mm256_sub_epi64(a0,b0); a1=_mm256_sub_epi64(a1,b1); } }
if (nt >= 4) { __m256i b0 = LDU(r3 + jb), b1 = LDU(r3 + jb + 4);
if (g3) { a0=_mm256_add_epi64(a0,b0); a1=_mm256_add_epi64(a1,b1); }
else { a0=_mm256_sub_epi64(a0,b0); a1=_mm256_sub_epi64(a1,b1); } }
STU(q0, a0);
STU(q0 + (size_t)KC * 4, a1);
}
}
jb0 = ngrp * 8;
for (int k = kc; k < KC; k++) {
long long *qb = Bpanel + (size_t)k * 4;
for (int g = 0; g < ngrp; g++)
{ long long *q0 = qb + (size_t)g * 8 * KC; STU(q0, _mm256_setzero_si256()); STU(q0 + (size_t)KC * 4, _mm256_setzero_si256()); }
}
}
if (nt <= 2) {
const long long *s0 = P + (size_t)o->r[0] * ld + o->c[0] + jc;
const long long *s1 = (nt == 2) ? P + (size_t)o->r[1] * ld + o->c[1] + jc : 0;
int ngrp = 0;
for (int jb = 0; jb + 8 <= nc; jb += 8) ngrp++;
int k = 0;
/* 两行一组：同一块里 k 与 k+1 的记录相邻 32 B = 同一条 64 B line，
两个 32 B 写背靠背 ⇒ 整行写。布局/输出逐字节不变。 */
for (; k + 2 <= kc; k += 2) {
const long long *r0 = s0 + (size_t)(pc + k) * ld;
const long long *r0b = s0 + (size_t)(pc + k + 1) * ld;
const long long *r1 = (nt == 2) ? s1 + (size_t)(pc + k) * ld : 0;
const long long *r1b = (nt == 2) ? s1 + (size_t)(pc + k + 1) * ld : 0;
long long *qb = Bpanel + (size_t)k * 4;
for (int g = 0; g < ngrp; g++) {
long long *q0 = qb + (size_t)g * 2 * KC * 4;
__m256i a0 = LDU(r0 + g * 8), a1 = LDU(r0 + g * 8 + 4);
__m256i c0 = LDU(r0b + g * 8), c1 = LDU(r0b + g * 8 + 4);
if (nt == 2) {
__m256i b0 = LDU(r1 + g * 8), b1 = LDU(r1 + g * 8 + 4);
__m256i d0 = LDU(r1b + g * 8), d1 = LDU(r1b + g * 8 + 4);
if (o->s[1] > 0) { a0=_mm256_add_epi64(a0,b0); a1=_mm256_add_epi64(a1,b1);
c0=_mm256_add_epi64(c0,d0); c1=_mm256_add_epi64(c1,d1); }
else { a0=_mm256_sub_epi64(a0,b0); a1=_mm256_sub_epi64(a1,b1);
c0=_mm256_sub_epi64(c0,d0); c1=_mm256_sub_epi64(c1,d1); }
}
STU(q0, a0);
STU(q0 + 4, c0);
STU(q0 + (size_t)KC * 4, a1);
STU(q0 + (size_t)KC * 4 + 4, c1);
}
}
for (; k < kc; k++) {
const long long *r0 = s0 + (size_t)(pc + k) * ld;
const long long *r1 = (nt == 2) ? s1 + (size_t)(pc + k) * ld : 0;
long long *qb = Bpanel + (size_t)k * 4;
for (int g = 0; g < ngrp; g++) {
const int jb = g * 8;
long long *q0 = qb + (size_t)(jb / 4) * KC * 4;
__m256i a0 = LDU(r0 + jb), a1 = LDU(r0 + jb + 4);
if (nt == 2) {
__m256i b0 = LDU(r1 + jb), b1 = LDU(r1 + jb + 4);
if (o->s[1] > 0) { a0 = _mm256_add_epi64(a0, b0); a1 = _mm256_add_epi64(a1, b1); }
else { a0 = _mm256_sub_epi64(a0, b0); a1 = _mm256_sub_epi64(a1, b1); }
}
STU(q0, a0);
STU(q0 + (size_t)KC * 4, a1);
}
}
jb0 = ngrp * 8;
for (int k = kc; k < KC; k++) {
long long *qb = Bpanel + (size_t)k * 4;
for (int g = 0; g < ngrp; g++)
{ long long *q0 = qb + (size_t)g * 8 * KC; STU(q0, _mm256_setzero_si256()); STU(q0 + (size_t)KC * 4, _mm256_setzero_si256()); }
}
}
if (nt >= 3) jb0 = NC;
for (int jb = jb0; jb < NC; jb += 4) {
long long *q = Bpanel + (size_t)(jb / 4) * KC * 4;
int full = (jb + 4 <= nc);
int lim = nc - jb; if (lim > 4) lim = 4; if (lim < 0) lim = 0;
const long long *sq = P + (size_t)o->r[0] * ld + o->c[0] + jc + jb;
const long long *sq1 = (nt >= 2) ? P + (size_t)o->r[1] * ld + o->c[1] + jc + jb : 0;
if (full && nt <= 2) {
for (int k = 0; k < kc; k++) {
__m256i v = LDU(sq + (size_t)(pc + k) * ld);
if (nt == 2) {
__m256i w = LDU(sq1 + (size_t)(pc + k) * ld);
v = (o->s[1] > 0) ? _mm256_add_epi64(v, w) : _mm256_sub_epi64(v, w);
}
STU(q + (size_t)k * 4, v);
}
for (int k = kc; k < KC; k++) STU(q + (size_t)k * 4, _mm256_setzero_si256());
} else {
for (int k = 0; k < KC; k++) {
long long t[4] = {0, 0, 0, 0};
if (k < kc) {
int ncol = full ? 4 : lim;
for (int c = 0; c < ncol; c++) {
long long v = 0;
for (int tt = 0; tt < nt; tt++)
v += o->s[tt] < 0 ? -P[(size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb + c]
: P[(size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb + c];
t[c] = v;
}
}
__m256i v = LDU(t);
STU(q + (size_t)k * 4, v);
}
}
}
}
static inline int effop(int base, int term) { /* term: 1 = +M, 2 = -M ; result: 0 store 1 add 2 sub */
if (base == 0) return term == 1 ? 0 : 2;
if (base == 1) return term == 2 ? 2 : 1;
return term == 2 ? 1 : 2;
}
/* C1 op1 (+- second destination C2 op2, either may be absent) = Ao * Bo, block size n */
static void leaf(int n, const Opd *Ao, const Opd *Bo,
long long *C1, int ldc1, int op1, long long *C2, int ldc2, int op2) {
int nb4 = n - (n & 3);
for (int jc = 0; jc < nb4; jc += NC) {
int nc = (nb4 - jc < NC) ? (nb4 - jc) : NC;
for (int pc = 0; pc < n; pc += KC) {
int kc = (n - pc < KC) ? (n - pc) : KC;
packB(Bo, jc, nc, pc, kc);
int e1 = (pc == 0) ? op1 : (op1 == 2 ? 2 : 1);
int e2 = (pc == 0) ? op2 : (op2 == 2 ? 2 : 1);
for (int ic = 0; ic < n; ic += MC) {
int mc = (n - ic < MC) ? (n - ic) : MC;
packA(Ao, ic, mc, pc, kc);
for (int ir = 0; ir < mc; ir += MR) {
int rows = (mc - ir < MR) ? (mc - ir) : MR;
const long long *ap = Apanel + (size_t)ir * KC * 1;
for (int jg = 0; jg < nc; jg += 4) {
const long long *bp = Bpanel + (size_t)(jg / 4) * KC * 4;
mikro(KC / UU, ap, bp, C1 + (size_t)(ic + ir) * ldc1 + jc + jg, ldc1, e1,
C2 ? C2 + (size_t)(ic + ir) * ldc2 + jc + jg : 0, ldc2, e2, rows);
}
}
}
}
}
for (int j = nb4; j < n; j++)
for (int i = 0; i < n; i++) {
unsigned long long acc = 0;
for (int k = 0; k < n; k++) {
long long av = 0, bv = 0;
for (int t = 0; t < Ao->nt; t++)
av += Ao->s[t] < 0 ? -Ao->p[(size_t)(Ao->r[t] + i) * Ao->ld + Ao->c[t] + k]
: Ao->p[(size_t)(Ao->r[t] + i) * Ao->ld + Ao->c[t] + k];
for (int t = 0; t < Bo->nt; t++)
bv += Bo->s[t] < 0 ? -Bo->p[(size_t)(Bo->r[t] + k) * Bo->ld + Bo->c[t] + j]
: Bo->p[(size_t)(Bo->r[t] + k) * Bo->ld + Bo->c[t] + j];
acc += (unsigned long long)av * (unsigned long long)bv;
}
long long *cp = C1 + (size_t)i * ldc1 + j;
if (op1 == 0) *cp = (long long)acc;
else if (op1 == 1) *cp += (long long)acc;
else *cp -= (long long)acc;
if (C2) {
long long *cq = C2 + (size_t)i * ldc2 + j;
if (op2 == 0) *cq = (long long)acc;
else if (op2 == 1) *cq += (long long)acc;
else *cq -= (long long)acc;
}
}
}
/* ------------------------------------------------------------------ *
* 本发：单目的地递归 Strassen（移植自同族 mmml4k 的 `mmxR`）。
* 旧 `mmx` 的守卫里有 `d1`（第二目的地）⇒ **M1..M5 一进来就落 leaf**，
* 于是"递归"与"双目的地"绑死：本题 n=1024 实际只有 1 层 Strassen。
* `mmxR` 只有一个目的地，第二目的地用 `cmb`（先算进 scratch 再加减）实现
* ⇒ 每一层都能继续递归。`Opd` 的 nt <= 4 正好够：顶层 nt=1 → 512 层 nt=2 →
* 256 层 nt=4 ⇒ 在 256 层落 leaf（守卫 nt > 2）。
* ------------------------------------------------------------------ */
/* scratch：n=1024 / CUTL=256 时，顶层 mmxR(1024) 需要 T 与 T2 = T + 512*512
⇒ 2 * 512^2 个 long。n<=1024 时递归深度最多用到这一对。 */
static long long g_scratch[2 * (512 * 512) + 4096] __attribute__((aligned(64)));
/* dst op= sign * T , (h x ldt) -> (h x ldd) , sign = +-1 , op in {0,1,2} */
static void cmb(long long *dst, int ldd, const long long *T, int ldt, int h, int sign, int op) {
for (int i = 0; i < h; i++) {
const long long *t = T + (size_t)i * ldt;
long long *c = dst + (size_t)i * ldd;
int j = 0;
if (op == 0) {
for (; j + 4 <= h; j += 4) {
__m256i v = LDQ((const __m256i *)(t + j));
if (sign < 0) v = _mm256_sub_epi64(_mm256_setzero_si256(), v);
STQ((__m256i *)(c + j), v);
}
for (; j < h; j++) c[j] = sign < 0 ? -t[j] : t[j];
} else {
for (; j + 4 <= h; j += 4) {
__m256i v = LDQ((const __m256i *)(t + j));
__m256i w = LDQ((const __m256i *)(c + j));
if ((op == 1) != (sign < 0)) STQ((__m256i *)(c + j), _mm256_add_epi64(w, v));
else STQ((__m256i *)(c + j), _mm256_sub_epi64(w, v));
}
for (; j < h; j++) { if ((op == 1) != (sign < 0)) c[j] += t[j]; else c[j] -= t[j]; }
}
}
}
/* single-destination recursive Strassen: d op= Ao*Bo , T = scratch block (>= 2*h*h) */
static void mmxR(int n, const Opd *Ao, const Opd *Bo, long long *d, int ld, int op, long long *T) {
if (n <= CUTL || (n & 1) || Ao->nt > 2 || Bo->nt > 2) {
leaf(n, Ao, Bo, d, ld, op, 0, 0, 0); return; }
const int h = n >> 1;
Opd A11 = oq(Ao, 0, 0, h), A12 = oq(Ao, 0, 1, h), A21 = oq(Ao, 1, 0, h), A22 = oq(Ao, 1, 1, h);
Opd B11 = oq(Bo, 0, 0, h), B12 = oq(Bo, 0, 1, h), B21 = oq(Bo, 1, 0, h), B22 = oq(Bo, 1, 1, h);
long long *Q[4];
Q[0] = d; Q[1] = d + h; Q[2] = d + (size_t)h * ld; Q[3] = Q[2] + h;
long long *T2 = T + (size_t)h * h;
int wr[4] = {0, 0, 0, 0};
#define MK(qi, si, qj, sj, Ma, Mb) do { \
int oa = wr[qi] ? effop(1, (si)) : effop(op, (si)); \
int ob = wr[qj] ? effop(1, (sj)) : effop(op, (sj)); \
wr[qi] = 1; wr[qj] = 1; \
if (h <= CUTL) { \
/* sub-call IS a leaf: hand it BOTH destination quadrants directly \
(one mikro pass, no scratch, no cmb) */ \
leaf(h, &(Ma), &(Mb), Q[qi], ld, oa, Q[qj], ld, ob); \
} else if (oa == 0) { \
mmxR(h, &(Ma), &(Mb), Q[qi], ld, 0, T); \
cmb(Q[qj], ld, Q[qi], ld, h, (sj), ob); \
} else { \
mmxR(h, &(Ma), &(Mb), T, h, 0, T2); \
cmb(Q[qi], ld, T, h, h, (si), oa); \
cmb(Q[qj], ld, T, h, h, (sj), ob); \
} \
} while (0)
#define MK1(qi, si, Ma, Mb) do { \
int oa = wr[qi] ? effop(1, (si)) : effop(op, (si)); \
wr[qi] = 1; \
mmxR(h, &(Ma), &(Mb), Q[qi], ld, oa, T); \
} while (0)
{ /* M1 = (A11+A22)(B11+B22) -> +C11, +C22 */
Opd Ma = oadd(&A11, &A22, 1), Mb = oadd(&B11, &B22, 1);
MK(0, 1, 3, 1, Ma, Mb);
}
{ /* M2 = (A21+A22)B11 -> +C21, -C22 */
Opd Ma = oadd(&A21, &A22, 1), Mb = B11;
MK(2, 1, 3, 2, Ma, Mb);
}
{ /* M3 = A11(B12-B22) -> +C12, +C22 */
Opd Ma = A11, Mb = oadd(&B12, &B22, -1);
MK(1, 1, 3, 1, Ma, Mb);
}
{ /* M4 = A22(B21-B11) -> +C11, +C21 */
Opd Ma = A22, Mb = oadd(&B21, &B11, -1);
MK(0, 1, 2, 1, Ma, Mb);
}
{ /* M5 = (A11+A12)B22 -> -C11, +C12 */
Opd Ma = oadd(&A11, &A12, 1), Mb = B22;
MK(0, 2, 1, 1, Ma, Mb);
}
{ /* M6 = (A21-A11)(B11+B12) -> +C22 */
Opd Ma = oadd(&A21, &A11, -1), Mb = oadd(&B11, &B12, 1);
MK1(3, 1, Ma, Mb);
}
{ /* M7 = (A12-A22)(B21+B22) -> +C11 */
Opd Ma = oadd(&A12, &A22, -1), Mb = oadd(&B21, &B22, 1);
MK1(0, 1, Ma, Mb);
}
#undef MK
#undef MK1
}
/* ------------------------------------------------------------------ *
* 顶层改成 **物化版 Winograd 7 乘积**：
* 每个乘积的操作数是**单项**（nt=1）的子矩阵/物化和式 ⇒ 512 级递归里的 256 叶
* 看到的 nt 从 <=4 降到 <=2（pack 的源读项数减半），而 7 个和式（P/Q）在 4 MB 缓冲里
* **增量改写**、每个只算一次；最后 7 个乘积在**一趟向量**里合并进 C（不再有 cmb 的 RMW 链）。
* ------------------------------------------------------------------ */
/* sm_thp: 2 MB-aligned backing store + MADV_HUGEPAGE, so the 18.9 MB of Winograd
buffers costs ~9 faults instead of ~4700. */
#define SM_WG_ELEMS (9 * 512 * 512 + 9 * 32)
static long long wg_buf_store[9 * (512 * 512 + 3072 + 16)] __attribute__((aligned(1 << 21)));
static long long *wg_buf = wg_buf_store;
__attribute__((constructor)) static void sm_thp_init(void) {
uintptr_t a = (uintptr_t)wg_buf_store;
uintptr_t e = a + (uintptr_t)sizeof(wg_buf_store);
uintptr_t ha = (a + ((uintptr_t)1 << 21) - 1) & ~(((uintptr_t)1 << 21) - 1);
uintptr_t he = e & ~(((uintptr_t)1 << 21) - 1);
if (he > ha) madvise((void *)ha, (size_t)(he - ha), MADV_HUGEPAGE);
}
static void wg_add(int h, long long *dst, const long long *a, int lda,
const long long *b, int ldb, int sg) {
for (int i = 0; i < h; i++) {
long long *d = dst + (size_t)i*h;
const long long *x = a + (size_t)i*lda, *y = b + (size_t)i*ldb;
for (int j = 0; j < h; j += 4) {
__m256i vx=LDU(x+j), vy=LDU(y+j);
STU(d+j, sg > 0 ? _mm256_add_epi64(vx,vy) : _mm256_sub_epi64(vx,vy));
}
}
}
static void wg_mutate(int h, long long *dst, const long long *b, int ldb, int reverse) {
for (int i = 0; i < h; i++) {
long long *d = dst + (size_t)i*h; const long long *y=b+(size_t)i*ldb;
for (int j=0; j<h; j+=4) {
__m256i vd=LDU(d+j),vy=LDU(y+j);
STU(d+j, reverse ? _mm256_sub_epi64(vy,vd) : _mm256_sub_epi64(vd,vy));
}
}
}
static void wg_combine(int h, long long *C, int ldc, const long long *M1,
const long long *M2, const long long *M3, const long long *M4,
const long long *M5, const long long *M6, const long long *M7) {
for (int i=0; i<h; i++) {
const long long *a=M1+(size_t)i*h,*b=M2+(size_t)i*h,*c=M3+(size_t)i*h;
const long long *d=M4+(size_t)i*h,*e=M5+(size_t)i*h,*f=M6+(size_t)i*h,*g=M7+(size_t)i*h;
long long *c11=C+(size_t)i*ldc,*c12=c11+h,*c21=c11+(size_t)h*ldc,*c22=c21+h;
/* 7 条 M 读流各"每行踢一脚"：行内 64 B 步长由硬件跟踪，行间 4 KB 的跳越超出
L2 流预取器的能力 ⇒ 每行补一条非临时预取（M 只被读一次 ⇒ nta 不污染 cache）。
代价 7 uop / 行（行体 512 元素）≈ 0.014 uop/元素。 */
if (i + 4 < h) {
long long *pa=(long long*)a+(size_t)4*h, *pb=(long long*)b+(size_t)4*h;
long long *pc=(long long*)c+(size_t)4*h, *pd=(long long*)d+(size_t)4*h;
long long *pe=(long long*)e+(size_t)4*h, *pf=(long long*)f+(size_t)4*h;
long long *pg=(long long*)g+(size_t)4*h;
__builtin_prefetch(pa,0,0); __builtin_prefetch(pb,0,0); __builtin_prefetch(pc,0,0);
__builtin_prefetch(pd,0,0); __builtin_prefetch(pe,0,0); __builtin_prefetch(pf,0,0);
__builtin_prefetch(pg,0,0);
}
for (int j=0; j<h; j+=4) {
__m256i p1=LDU(a+j),p2=LDU(b+j),p3=LDU(c+j),p4=LDU(d+j);
__m256i p5=LDU(e+j),p6=LDU(f+j),p7=LDU(g+j);
__m256i u2=_mm256_add_epi64(p1,p6),u3=_mm256_add_epi64(u2,p7);
__m256i u4=_mm256_add_epi64(u2,p5);
STNT(c11+j,_mm256_add_epi64(p1,p2));
STNT(c12+j,_mm256_add_epi64(u4,p3));
STNT(c21+j,_mm256_sub_epi64(u3,p4));
STNT(c22+j,_mm256_add_epi64(u3,p5));
}
}
}
/* ---- sm_mid: materialised Winograd at the 512 level (a further layer) ---- */
static long long wg_mid[9 * (256 * 256 + 3072 + 16)] __attribute__((aligned(1 << 21)));
__attribute__((constructor)) static void sm_thp_mid_init(void) {
uintptr_t a = (uintptr_t)wg_mid, e = a + (uintptr_t)sizeof(wg_mid);
uintptr_t ha = (a + ((uintptr_t)1 << 21) - 1) & ~(((uintptr_t)1 << 21) - 1);
uintptr_t he = e & ~(((uintptr_t)1 << 21) - 1);
if (he > ha) madvise((void *)ha, (size_t)(he - ha), MADV_HUGEPAGE);
}
static void winograd_mid(int n, const long long *X, int ldx, const long long *Y, int ldy,
long long *Z, int ldz) {
const int g = n >> 1; const size_t z2 = (size_t)g * g;
long long *m1=wg_mid,*m2=m1+z2+3072,*m3=m2+z2+3072,*m4=m3+z2+3072;
long long *m5=m4+z2+3072,*m6=m5+z2+3072,*m7=m6+z2+3072,*pp=m7+z2+3072,*qq=pp+z2+3072;
const long long *X11=X,*X12=X+g,*X21=X+(size_t)g*ldx,*X22=X21+g;
const long long *Y11=Y,*Y12=Y+g,*Y21=Y+(size_t)g*ldy,*Y22=Y21+g;
#define MPROD(A_,LNA_,B_,LNB_,Z_) do { Opd ao=o1((A_),(LNA_)),bo=o1((B_),(LNB_)); mmxR(g,&ao,&bo,(Z_),g,0,g_scratch); } while(0)
MPROD(X11,ldx,Y11,ldy,m1);
MPROD(X12,ldx,Y21,ldy,m2);
wg_add(g,pp,X21,ldx,X22,ldx,+1); wg_add(g,qq,Y12,ldy,Y11,ldy,-1); MPROD(pp,g,qq,g,m5);
wg_mutate(g,pp,X11,ldx,0); wg_mutate(g,qq,Y22,ldy,1); MPROD(pp,g,qq,g,m6);
wg_mutate(g,pp,X12,ldx,1); MPROD(pp,g,Y22,ldy,m3);
wg_mutate(g,qq,Y21,ldy,0); MPROD(X22,ldx,qq,g,m4);
wg_add(g,pp,X11,ldx,X21,ldx,-1); wg_add(g,qq,Y22,ldy,Y12,ldy,-1); MPROD(pp,g,qq,g,m7);
#undef MPROD
wg_combine(g,Z,ldz,m1,m2,m3,m4,m5,m6,m7);
}
static void winograd_top(int n, const long long *A, const long long *B, long long *C) {
const int h=n/2; const size_t sz=(size_t)h*h;
long long *M1=wg_buf,*M2=M1+sz+3072,*M3=M2+sz+3072,*M4=M3+sz+3072;
long long *M5=M4+sz+3072,*M6=M5+sz+3072,*M7=M6+sz+3072,*P=M7+sz+3072,*Q=P+sz+3072;
const long long *A11=A,*A12=A+h,*A21=A+(size_t)h*n,*A22=A21+h;
const long long *B11=B,*B12=B+h,*B21=B+(size_t)h*n,*B22=B21+h;
#define PROD(X,Y,Z,LX,LY) do { winograd_mid(h,(X),(LX),(Y),(LY),(Z),h); } while(0)
PROD(A11,B11,M1,n,n);
PROD(A12,B21,M2,n,n);
wg_add(h,P,A21,n,A22,n,+1); /* S1 */
wg_add(h,Q,B12,n,B11,n,-1); /* T1 */
PROD(P,Q,M5,h,h);
wg_mutate(h,P,A11,n,0); /* S2 = S1 - A11 */
wg_mutate(h,Q,B22,n,1); /* T2 = B22 - T1 */
PROD(P,Q,M6,h,h);
wg_mutate(h,P,A12,n,1); /* S4 = A12 - S2 */
PROD(P,B22,M3,h,n);
wg_mutate(h,Q,B21,n,0); /* T4 = T2 - B21 */
PROD(A22,Q,M4,n,h);
wg_add(h,P,A11,n,A21,n,-1); /* S3 */
wg_add(h,Q,B22,n,B12,n,-1); /* T3 */
PROD(P,Q,M7,h,h);
#undef PROD
wg_combine(h,C,n,M1,M2,M3,M4,M5,M6,M7);
_mm_sfence();
}
void matrix_multiply(int n, const long long *A, const long long *B, long long *C) {
Opd Ao = o1(A, n), Bo = o1(B, n);
if (n <= CUTL || (n & 1)) leaf(n, &Ao, &Bo, C, n, 0, 0, 0, 0);
else if (n == 1024) winograd_top(n, A, B, C);
else mmxR(n, &Ao, &Bo, C, n, 0, g_scratch);
}
#pragma GCC pop_options
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 111.019 ms | 30 MB + 668 KB | Accepted | Score: 100 | 显示更多 |