/* References:
- saffah_cc_v41_agg1, https://duck.ac/submission/123578: copied current indexed modular64 engine, retaining its complete inherited references and source notices below.
No separate software license notice appears in the public source.
Idea: The six-row hot loop already broadcasts next iteration's A0 into ymm14 at its end, after the previous product consumes that register. The next iteration redundantly broadcasts the same address again. Prime ymm14 once before the loop and use the carried broadcast thereafter; all six actual A products and the existing early load remain unchanged.
Purpose: Eliminate one redundant broadcast per contraction step while preserving the existing load latency overlap and exact arithmetic.
*/
// ===== [ml1g_] mmml1k 席 ===== **刀 K-A：mikro 热环的「入口对齐 + 不执行填充」**
//
// ===== REFERENCES（依 RULES.md §3；引用他人代码列原作者账号 + 原提交 URL + 所用内容）=====
// [1] duck.ac 用户 **saffah_cc_v41_260924**，提交 **#122376**
// <https://duck.ac/submission/122376>（Accepted，106.718334 ms = 本题榜面 T）——
// 本件正文的**最底层基座**（逐字节复制；Strassen 递归、mikro 汇编环、面板布局、
// buildA_blk/buildB/cmb2/cmb3/cmb*_blk 全部保留）。
// [2] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#123056**（106.205045 ms =
// 我方现役 mine）<https://duck.ac/submission/123056> —— 本件正文的**直接基座**（逐字节
// 复制）。相对 [1] 它只加了本账号 ml1d_ 席的「冷写 RFO」刀（`matop` n>=512 与 `cmb3`
// dst>=256^2 用 `vmovntdq`，见 ml1d_ 席 notes 节）。**本件不动那一刀的任何一个字节。**
// [3] 仪器（本账号自有，非引用）：`tools/exact.py`（判题机 ns 级真值/T 与严支）、
// `tools/probe.py`（立即体验通道 = 判题机 in-process 整进程时间）、`ref/gcc9/g9.sh`
// （判题同款 gcc 9.3，`-O2 -std=c++17 -static -U_FORTIFY_SOURCE`）。
//
// ===== 新思路（ml1g_ 席）=====
// 【真靶】开棒现读 `tools/exact.py mmml1k`：`mine #123056 = 106.205045` · `T #122376 =
// 106.718334` · 严支 **0.99·T+1µs = 105.652151** ⇒ 缺 **0.552894 ms（0.521%）· RED**。
//
// 【本刀的唯一改动】`mikro` 的汇编里，在「已剥出的第一次迭代」与热环标号 `1:` 之间原本
// 有 **56 字节 `nop` 填充**（作者在原始引擎里用来把热环对齐）。它在**落穿路径上**：
// 每次调用（n=1024 路径上 422 576 次）都要发射这 7 条 nop。
// 本刀把它换成 **`jmp 1f` + `.p2align 6`** 两个指令级的改动：
// ① `jmp 1f` → 落穿不再经过填充（7 条 nop uop → 1 条 jmp uop）；
// ② `.p2align 6` → 由汇编器把热环起点对齐到 **64 字节边界**，不再依赖「函数起始
// 地址 mod 64」这个**会被任何 .text 变更平移的量**（本席实测：现役件里 mikro 起始
// mod64 = 0、热环起点 mod64 = **3**；本刀后热环起点 mod64 = **0**，且与函数落点无关）。
// **语义零变化**：`jmp 1f`/`.p2align` 不触碰任何寄存器、内存或标志位；热环体、A/B 面板
// 读图案、所有存指令一字未动。**逐位等价闸门**：判题同款 gcc 9.3，n = 3,4,5,7,8,9,11,
// 15,16,17,23,31,32,33,47,63,64,65,96,127,128,129,192,255,256,257,384,512,700,1024
// 共 30 档，FNV 与 SUM **全档逐位同基座** ✓（n=1024 的 FNV = 10686382741757082234 与
// 在册 ck 一致）。
//
// 【判据 / 定价】判题机 in-process 三臂篮（`work/ml1g__tmp/ml1g_arms3s.cpp`，驱动钉在引擎
// 之前；同二进制内三臂逐位同 ck）读到本刀族 −0.41 ~ −0.54 ms；但本席另用**对齐受控**的
// 三臂篮（`ml1g_arms6s.cpp`，所有臂的 mikro 都 `aligned(64)` + 环前 `.p2align 6`，实测
// 四份 mikro 的热环**全部落 mod64 = 0**、环长 201 B）复测，**同一批内 A0(现役对齐+填充)
// = 106.887 ms** ⇒ 本刀（=A0 的代码 + jmp + .p2align 6）的板面价**必须用板面判**（探针
// 装置对「改 .text 体积」的候选有 ~0.9% 的布局敏感性，本席两篮读数符号相反 ⇒ 不可用）。
// ⇒ **本件是「板面定价件」**：配合同批锚 `ml1g_anchor.cpp`（= #123056 逐字节 + `MR 6`
// → `(6)` 的词法 no-op，机码与 #123056 逐字节同）分离刀价与窗口。
//
// 【实验性提交声明】本件为**实验性提交**：目的是在**判题机板面**（ns 级、已证窗口 ±0.06 ms）
// 上给「热环入口对齐 + 去掉执行填充」这一个单点定号。若板面 ≤ 现 mine 则本件即新 mine；
// 若 > 现 mine 则本件不进榜面最快件，`mine` 不受影响。
// ===== [ml1d_] mmml1k ===== **刀：冷写的 RFO —— 按作用域律只铺 n≥512 级（非临时存 `vmovntdq`）** =====
//
// ===== REFERENCES（依 RULES.md §3；引用他人代码列原作者账号 + 原提交 URL + 所用内容）=====
// [1] duck.ac 用户 **saffah_cc_v41_260924**，提交 **#122376**
// <https://duck.ac/submission/122376>（Accepted，106.718334 ms = 本题榜面 T）——
// **本件正文的基座**（逐字节复制；其正文里的一切注释、宏、面板布局、Strassen 递归、
// mikro 汇编环、buildA_blk/buildB/cmb2/cmb3/cmb*_blk 全部保留）。其自述刀：
// `l1k6_rv swp2`（B 面板存预交换形态、环内两条相隔 4 B 的载入替掉 `vpshufd`）、
// `l1k15_rv fus3`（叶相位三目标融合）、`l1k7_rv rba`（A 面板按 MR 行分块重建）。
// [2] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#122780**（106.625010 ms = 我方
// 现役 mine；机器码与 [1] 逐字节同）—— 本件的**对照件 / 同批锚**（见 §闸门）。
// [3] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），**mmml2k** 提交 **#123023**
// <https://duck.ac/submission/123023>（749.966163 ms，−4.94 ms）—— **刀法来源**（本账号
// 自有的另一题，非逐字复制）：该件给出「**冷写的 RFO 是本族唯一被板面兑现的写侧杠杆**」
// 与**作用域判据**（该级平面在消费者读到它之前是否已被逐出；实测 8 MB/2 MB 级付钱、
// 512 KB 级 `vmovntdq` 反而 +6.5 ms）。**本件把该判据按 1k 的尺寸逐级重算后只铺在付钱的那一级。**
// [4] 仪器（本账号自有，非引用）：`tools/exact.py`（判题机 ns 级真值）、`tools/probe.py`
// （立即体验通道 = 判题机 in-process 整进程时间）、本席生成器
// `work/ml1d__tmp/{ml1d_gen.py,ml1d_gen2.py,ml1d_mksub.py}`；闸门 = 判题同款 gcc 9.3
// `-O2 -std=c++17`。**本件的每一支臂的实测读数出自本席本轮的探针（见 §定价）。**
// ===== 新思路 =====（ml1d_ 席的工作）
// 【真靶】开棒现读 `tools/exact.py mmml1k` + `tools/subhist.py mmml1k` + 榜面逐行复核：
// `mine #122780 = 106.625010` · `T #122376 = 106.718334` · 严支 **0.99·T+1µs = 105.652151**
// ⇒ 缺 **0.972859 ms（0.912%）**。（派单里的 108.884787/107.483462/1.401325 是 **ml1y7_ 时期
// 的两代前口径**，与判题机现读不符；以现读为准。）
// 【判据：谁的钱】判题机 = i3-8100，L1d 32 KB / L2 256 KB / **L3 6 MB**。NT 只对
// 「**写后不再原样读回、且写出前已被逐出**」的块付钱（本族在册判词：`#104960` NT 写 C = −0.27% ✓、
// `ntadd` NT 写 M 目的地 = +0.53% ✗，差别就在消费者还读不读这份数据）。
// 本席用计数器把本题**全部写侧账量清**（`work/ml1d__tmp/ml1d_cnt4.cpp`，n=1024 实跑）：
// `matop` 80 呼 / **55.000 MB 纯存**（n=512 十呼 = 20 MB；n=256 七十呼 = 35 MB）
// `cmb3` 2736 呼 / 16.531 MB 纯存，其中 dst ≥ 256×256 只有 **8.875 MB**（其余 2695 呼是叶级热瓦片）
// `cmb2` 162 呼 / 63.000 MB、`cmb2_blk` 10780 呼 / 61.250 MB —— **全是 RMW（先读后写，
// 行已独占、本来就无 RFO）⇒ NT 在结构上不可用，不在本刀域内**
// 【作用域律的 1k 版（逐级数出来的）】
// * n=**512** 的 `matop` dst = **2 MB**：写完还要先写另一侧 2 MB 才轮到消费者，之后 7 个 256-子节点
// 各自再物化 1 MB + 各跑 7 个叶（buildB 每叶 229 KB 写）⇒ 写—读间隔 ≈ 6~7 MB > L3 6 MB
// ⇒ **逐出成立 = 冷** ⇒ 每写白付一次 RFO（20 MB）。
// * n=**256** 的 `matop` dst = **512 KB**：其 7 个 128-孙节点的 L3 足迹 ≈ 512 KB + 512 KB
// + 230 KB（B 面板回收复用）≈ **1.3 MB ≪ 6 MB** ⇒ **驻留 = 热**（消费者 122 MB 读全在 L3）
// ⇒ **NT 会把那 122 MB 打到 DRAM**。**这是本刀唯一的域边界，本席用两支探针臂把它量成了数。**
// 【本件改动】（只改**存指令形态**，值/次序/字节数/面板布局一字未动）
// ① `matop2`：`n >= 512` 时用 `vmovntdq`（消 20 MB 写-分配 RFO），`n == 256` 保持 `vmovdqu`；
// ② `cmb3`：dst 行数 ≥ 256 且 dst 32 B 对齐时用 `vmovntdq`（再消 8.875 MB）。
// `vmovntdq` 需 32 B 对齐：`matop` 的 dst = `matbuf[0/1]`（64 B 对齐，行距 n*8 = 2048/4096 B）✓；
// `cmb3` 在**入口**做一次性对齐 + `ldd` 检查，不满足的调用原样走可缓存存 ✓。单线程同核
// 「同地址 load-after-store」由 x86 内存模型保证（2k 的 #123023 已在 n=2048 全元素验证过同一条）。
// 【定价（判题机 in-process，本席本轮两支探针，A0=基座同二进制同机码）】
// 探针#1：A0 107.302 · A1(`matop` n≥512 NT) **106.826 → −0.476** · A2(再加 n=256 NT) 108.925
// → **+1.623** · A3(A1+cmb3≥256 NT) **106.629 → −0.673**
// 探针#2：A0 106.985 · A1 **106.482 → −0.503** · A2(+cmb3≥256) 106.486 → −0.499 ·
// A3(A2+cmb3≥128+目的地 `prefetchw`) 106.713 → −0.272
// ⇒ ★ **两发同号：`matop` n≥512 的 NT = −0.48/−0.50 ms（本刀的全部确定性收益）**；
// `cmb3` ≥256 = −0.197/+0.004（噪声内，取了低风险的一半）；**n=256 级 NT = +2.10 ms（判负，
// 与 mmml2k 的 512 KB 级 +6.5 ms 同形）**；目的地 `prefetchw` = +0.23（判负，未采纳）。
// 四臂 FNV 校验和逐位相同（`10686382741757082234`）✓。
// 【诚实交底】本刀只买到 ≈ **−0.50 ms**，而缺口是 0.972859 ms ⇒ **本件发射后预期仍 RED（约缺 0.47 ms）**，
// 发射的目的是**把已兑现的收益落袋 + 给下一手留下「哪一级付钱」的判题机实测边界**，
// 不是为了冲线。**未采纳项（n=256 NT / prefetchw）已判负，别再重走。**
// 【闸门】判题同款 gcc 9.3、`-O2 -std=c++17`：n = 128/256/512/1024 四档 FNV 与 [1]/[2] 逐位相同 ✓
// ===== [ml1b_] mmml1k ===== **同批基座锚（诊断件，非候选）** =====
//
// ===== REFERENCES（依 RULES.md §3；引用他人代码列原作者账号 + 原提交 URL + 所用内容）=====
// [1] duck.ac 用户 **saffah_cc_v41_260924**，提交 **#122376**
// <https://duck.ac/submission/122376>（Accepted，106.718334 ms = 本题榜面 T）——
// **本件正文的基座**（逐字节复制 + 本头部与 §思路 所述改动；其正文里的一切注释、宏、
// 面板布局、Strassen 递归、mikro 汇编环、buildA_blk/buildB/cmb* 全部保留）。
// 其自述刀：`l1k6_rv swp2`（B 面板存预交换形态、环内两条相隔 4 B 的载入替掉 `vpshufd`）、
// `l1k15_rv fus3`（叶相位三目标融合）、`l1k7_rv rba`（A 面板按 MR 行分块重建）。
// [2] duck.ac 本账号 **saffah_cc_v41_agg1**，提交 **#122723**（106.759485 ms = 我方现役 mine，
// 机器码与 [1] 逐字节同，`.text` md5 `a420b4f863263b3d`）—— 本件的**对照件**（闸门 + 探针标定）。
// [3] 仪器（本账号自有，非引用）：`tools/exact.py`（判题机 ns 级真值）、`tools/probe.py`
// （立即体验通道 = 判题机 in-process 整进程 Time (ms)）、本席生成器
// `work/ml1b__tmp/gen1.py`（从 [1] 的正文机器化生成本件）、闸门 = 判题同款 gcc 9.3 `-O2`。
// ===== 思路 =====（ml1b_ 席的新工作）
// 【真靶】T = #122376 = 106.718334 ⇒ **真靶 = 0.99·T+1µs = 105.652151**；mine = 106.759485
// ⇒ 缺 **1.107334 ms（1.037%）**。派单方向 =「64 级叶路径专用化 / 去通用代码」（mmml4k 的
// 对手链实测真价：定长打包器 + store-only 专用驱动器买 −29.9/−23.4 ms）。
// 【判题机相位表（ml1a_ 装置，基座口径）】引擎 106.612 = mikro 81.37（76.3%）· matop 8.41 ·
// buildB 5.09 · cmb2 4.85 · cmb3 1.46 · cmb2_blk 1.05 · cmb3_blk 0.30 · 其余 4.1。
// 【本刀的单点性】mikro 是 `static inline`，其**汇编环逐字节未动**（34 指令 / 39 uop / 24 MAC，
// 环体字节数与指令次序完全一致），改的只有**每次调用的调用方/收尾侧**：
// ① 形参 10 个 → 5 个：原形参表 `(cnt,ap,bp,c1,ldc1,op1,c2,ldc2,op2,rows)` 里 **后 4 个走栈**
// （第 7~10 参）⇒ 每次调用都有 4 处栈写 + 栈读（`movl 16(%rbp),%r10d` 等）；
// ② 尾声的**逐行守卫链**：`if (rows > 0) … if (rows > 5)` = 6 组 `cmpl/je`；
// ③ 尾声的**逐行 c2 判空**：`testq %r9,%r9 / je` × 6；
// ④ 序言的 `cmp $4, %[rows] / je` 形状分派。
// 现改为 `template<int ROWS,int OP> mikro64s(cnt,ap,bp,c1,ldc1)`：**ROWS∈{6,4}、OP∈{0,1,2}
// 全部编译期**，尾声是 6（或 4）段直落代码（无分支、无守卫、无栈参）。
// **调用方 = leaflevel3 → leafP_blk，其 7 个调用点全部满足 `C2==0 ∧ n==KC==64 ∧ rows∈{6,4}
// ∧ op1∈{0,1,2}`**（L=0,1,2 存储；L=5,6 加；L=3,4 先落 T 再 cmb2_blk，也是 C2=0）；不满足该
// 条件的调用（其它 n、双目标）**原样走旧通用 mikro**，行为逐位不变。
// 【为什么这不是「环内减 uop」】环体一字未动 ⇒ 不触碰 ml4j_/ml4k_ 判死的「环体字节数不增」判据，
// 也不动端口配比；改掉的是**每 64 k 步一次**的固定开销（422 576 次/全程）。
// 【定价（判题机 in-process，同批；probe 双胞胎臂噪声 0.0086%）】基座 **112.475784 / 112.522063**
// ⇒ 本件 **111.699004 / 111.805461** = **−0.7768 / −0.7166 ms**（两发同号，均值 ≈ **−0.75 ms**）；
// 两件 stdout 校验和逐位同（`ck=10686382741757082234`）✓。按在册「减指令族 1:1 传导」
// （a14ah_：探针 −1.63% ↔ 板面 −1.54%）预期板面同量级。
// 【闸门】判题同款 gcc 9.3、`-O2 -std=c++17`：**30 档** n = 3,5,6,7,9,10,11,13,15,17,20,24,27,
// 33,39,40,64,100,128,129,130,200,256,260,400,511,512,700,1023,1024 的 FNV/SUM 与 [1]/[2]
// **逐位全同** ✓（含全部 n≥2·CUTL 的叶路径档）。`.text`：86 246 B → 87 974 B（+1.7 KB）。
// 【未动】mikro 汇编环、面板布局与尺寸、`matop`/`buildB`/`buildA_blk`/`cmb2`/`cmb3`/`cmb*_blk`、
// 全部宏常量（MR/KC/MC/NC/UU/CUTL/KCP）、Strassen 递归结构、C 输出寻址 ✓
// ============================================================================================
#pragma GCC optimize("O3","unroll-loops")
/* 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)
l1k6_rv `swp2`: the B panel is STORED PRE-SWAPPED and the mikro reads BOTH operands out
of it with two 32-byte loads 4 bytes apart -- the p5 `vpshufd` is deleted, one L1 load
replaces it, no second panel and no extra memory.
plus Strassen with fused operand combination and dual-destination stores
(no extra memory). MR=6 KC=512 MC=48 NC=512 CUTL=512 U=1 BSW=0 EXP=0 */
#include <immintrin.h>
/* Single-instruction 256-bit unaligned access, NO alignment precondition.
g++-9 -O2 with no -march expands the unaligned 256-bit intrinsics into 2-3 uops, one
of them a p5-only vinsertf128/vextractf128; these emit the single instruction the
intrinsic is named for. AT&T order is `vmovdqu src,dst`: the STORE is %0,%1 (ymm
first, memory second). A REVERSED store operand silently emits a LOAD instead --
that is how this was first written and it produced wrong output with no diagnostic.
volatile + "memory" is REQUIRED: a non-volatile asm with only an "m" input is pure
and gets eliminated. */
#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")
/* SWQ == `vpshufd $0xB1`: exchange the two 32-bit halves inside every 64-bit lane.
The B panel is STORED pre-swapped so the mikro can get both of the operands it needs
out of the SAME panel with two loads 4 bytes apart (see the mikro). */
#define SWQ(v) _mm256_shuffle_epi32((v), 0xB1)
#include <stddef.h>
#pragma GCC push_options
#pragma GCC target("avx2")
#define STQNT(p, v) __asm__ volatile("vmovntdq %0,%1" :: "x"(v), "m"(*(__m256i *)(p)) : "memory")
static void cmb3(long long *d, int ldd, const long long *a, const long long *b, const long long *c, int ld, int rows, int ncol);
typedef long long TE;
#define MR (6)
#define KC 64
#define MC ( 12 )
#define NC 512
#define UU 1
#define CUTL 64
/* A-panel row stride PADDED by 8 uint64 (64 B). With the natural stride KC*8 = 4096 B
every one of the 6 rows of a micro-tile maps to the SAME L1 set (4096/64 = 64 sets
apart, 64 sets in the cache) -- a textbook conflict. +8 uint64 makes the stride
4160 B = 65 lines, so consecutive rows land in consecutive sets. */
#define KCP (KC + 12)
static TE Apanel[(size_t)MC * KCP + 64] __attribute__((aligned(64)));
static TE Bpanel[(size_t)2 * (NC / 4 + 1) * KC * 8 + 64] __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) {
__m256i AC[12];
long idx = -(cnt << 3);
long long *apE = (long long *)ap + cnt;
long long *bpE = (long long *)bp + (cnt << 2);
__asm__ volatile(
"test %[idx], %[idx]\n\t"
"jz 8f\n\t"
"vmovdqu (%[bp],%[idx],4), %%ymm12\n\t"
"vmovdqu 4(%[bp],%[idx],4), %%ymm13\n\t"
"cmp $4, %[rows]\n\t"
"je 6f\n\t"
"vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %[c0]\n\t"
"vpmulld %%ymm12, %%ymm14, %[c1]\n\t"
"vbroadcastsd 608(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %[c2]\n\t"
"vpmulld %%ymm12, %%ymm14, %[c3]\n\t"
"vbroadcastsd 1216(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %[c4]\n\t"
"vpmulld %%ymm12, %%ymm14, %[c5]\n\t"
"vbroadcastsd 1824(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %[c6]\n\t"
"vpmulld %%ymm12, %%ymm14, %[c7]\n\t"
"vbroadcastsd 2432(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %[c8]\n\t"
"vpmulld %%ymm12, %%ymm14, %[c9]\n\t"
"vbroadcastsd 3040(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %[c10]\n\t"
"vpmulld %%ymm12, %%ymm14, %[c11]\n\t"
"vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
"vmovdqu 36(%[bp],%[idx],4), %%ymm13\n\t"
"add $8, %[idx]\n\t"
"jz 2f\n\t"
".byte 0x0f,0x1f,0x84,0x00,0x00,0x00,0x00,0x00, 0x0f,0x1f,0x84,0x00,0x00,0x00,0x00,0x00, 0x0f,0x1f,0x84,0x00,0x00,0x00,0x00,0x00, 0x0f,0x1f,0x84,0x00,0x00,0x00,0x00,0x00, 0x0f,0x1f,0x84,0x00,0x00,0x00,0x00,0x00, 0x0f,0x1f,0x84,0x00,0x00,0x00,0x00,0x00, 0x0f,0x1f,0x84,0x00,0x00,0x00,0x00,0x00\n\t"
"vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
"1:\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %[c1], %[c1]\n\t"
"vbroadcastsd 608(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %[c0], %[c0]\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %[c3], %[c3]\n\t"
"vbroadcastsd 1216(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %[c2], %[c2]\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %[c5], %[c5]\n\t"
"vbroadcastsd 1824(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %[c4], %[c4]\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %[c7], %[c7]\n\t"
"vbroadcastsd 2432(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %[c6], %[c6]\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %[c9], %[c9]\n\t"
"vbroadcastsd 3040(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %[c8], %[c8]\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %[c11], %[c11]\n\t"
"vbroadcastsd 8(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %[c10], %[c10]\n\t"
"vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
"vmovdqu 36(%[bp],%[idx],4), %%ymm13\n\t"
"add $8, %[idx]\n\t"
"jnz 1b\n\t"
"jmp 2f\n\t"
"6:\n\t"
"vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %[c0]\n\t"
"vpmulld %%ymm12, %%ymm14, %[c1]\n\t"
"vbroadcastsd 608(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %[c2]\n\t"
"vpmulld %%ymm12, %%ymm14, %[c3]\n\t"
"vbroadcastsd 1216(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %[c4]\n\t"
"vpmulld %%ymm12, %%ymm14, %[c5]\n\t"
"vbroadcastsd 1824(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %[c6]\n\t"
"vpmulld %%ymm12, %%ymm14, %[c7]\n\t"
"vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
"vmovdqu 36(%[bp],%[idx],4), %%ymm13\n\t"
"add $8, %[idx]\n\t"
"jz 2f\n\t"
"5:\n\t"
"vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %[c1], %[c1]\n\t"
"vbroadcastsd 608(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %[c0], %[c0]\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %[c3], %[c3]\n\t"
"vbroadcastsd 1216(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %[c2], %[c2]\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %[c5], %[c5]\n\t"
"vbroadcastsd 1824(%[ap],%[idx],1), %%ymm14\n\t"
"vpaddq %%ymm15, %[c4], %[c4]\n\t"
"vpmuludq %%ymm14, %%ymm13, %%ymm15\n\t"
"vpmulld %%ymm12, %%ymm14, %%ymm14\n\t"
"vpaddd %%ymm14, %[c7], %[c7]\n\t"
"vpaddq %%ymm15, %[c6], %[c6]\n\t"
"vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
"vmovdqu 36(%[bp],%[idx],4), %%ymm13\n\t"
"add $8, %[idx]\n\t"
"jnz 5b\n\t"
"jmp 2f\n\t"
"8:\n\t"
"vpxor %[c0], %[c0], %[c0]\n\t"
"vmovdqa %[c0], %[c1]\n\t"
"vmovdqa %[c0], %[c2]\n\t"
"vmovdqa %[c0], %[c3]\n\t"
"vmovdqa %[c0], %[c4]\n\t"
"vmovdqa %[c0], %[c5]\n\t"
"vmovdqa %[c0], %[c6]\n\t"
"vmovdqa %[c0], %[c7]\n\t"
"vmovdqa %[c0], %[c8]\n\t"
"vmovdqa %[c0], %[c9]\n\t"
"vmovdqa %[c0], %[c10]\n\t"
"vmovdqa %[c0], %[c11]\n\t"
"2:\n\t"
""
: [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx), [c0] "=x"(AC[0]), [c1] "=x"(AC[1]), [c2] "=x"(AC[2]), [c3] "=x"(AC[3]), [c4] "=x"(AC[4]), [c5] "=x"(AC[5]), [c6] "=x"(AC[6]), [c7] "=x"(AC[7]), [c8] "=x"(AC[8]), [c9] "=x"(AC[9]), [c10] "=x"(AC[10]), [c11] "=x"(AC[11])
: [rows] "r"((long)rows)
: "ymm12", "ymm13", "ymm14", "ymm15", "cc", "memory");
if (rows > 0) {
__m256i lo = AC[0], hi = AC[1];
__m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
__m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
long long *cp = c1 + (size_t)0 * ldc1;
if (op1 == 0) STQ((__m256i *)cp, res);
else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
if (c2) { cp = c2 + (size_t)0 * ldc2;
if (op2 == 0) STQ((__m256i *)cp, res);
else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
if (rows > 1) {
__m256i lo = AC[2], hi = AC[3];
__m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
__m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
long long *cp = c1 + (size_t)1 * ldc1;
if (op1 == 0) STQ((__m256i *)cp, res);
else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
if (c2) { cp = c2 + (size_t)1 * ldc2;
if (op2 == 0) STQ((__m256i *)cp, res);
else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
if (rows > 2) {
__m256i lo = AC[4], hi = AC[5];
__m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
__m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
long long *cp = c1 + (size_t)2 * ldc1;
if (op1 == 0) STQ((__m256i *)cp, res);
else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
if (c2) { cp = c2 + (size_t)2 * ldc2;
if (op2 == 0) STQ((__m256i *)cp, res);
else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
if (rows > 3) {
__m256i lo = AC[6], hi = AC[7];
__m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
__m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
long long *cp = c1 + (size_t)3 * ldc1;
if (op1 == 0) STQ((__m256i *)cp, res);
else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
if (c2) { cp = c2 + (size_t)3 * ldc2;
if (op2 == 0) STQ((__m256i *)cp, res);
else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
if (rows > 4) {
__m256i lo = AC[8], hi = AC[9];
__m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
__m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
long long *cp = c1 + (size_t)4 * ldc1;
if (op1 == 0) STQ((__m256i *)cp, res);
else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
if (c2) { cp = c2 + (size_t)4 * ldc2;
if (op2 == 0) STQ((__m256i *)cp, res);
else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
if (rows > 5) {
__m256i lo = AC[10], hi = AC[11];
__m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
__m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
long long *cp = c1 + (size_t)5 * ldc1;
if (op1 == 0) STQ((__m256i *)cp, res);
else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
if (c2) { cp = c2 + (size_t)5 * ldc2;
if (op2 == 0) STQ((__m256i *)cp, res);
else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
}
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[8], c[8], s[8];
};
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) {
const int nt = o->nt;
const long long *P = o->p; const int ld = o->ld;
for (int r = 0; r < MC; r++) {
long long *p = Apanel + (size_t)r * KCP;
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 = LDQ((const __m256i *)(s0 + k));
if (nt == 2) {
__m256i w = LDQ((const __m256i *)(s1 + k));
v = (sgn1 > 0) ? _mm256_add_epi64(v, w) : _mm256_sub_epi64(v, w);
}
STQ((__m256i *)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 {
for (; k + 16 <= kc; k += 16) {
__m256i v0 = _mm256_setzero_si256(), v1 = v0, v2 = v0, v3 = v0;
for (int t = 0; t < nt; t++) {
const long long *sp = P + (size_t)(o->r[t] + ic + r) * ld + o->c[t] + pc + k;
__m256i w0 = LDQ((const __m256i *)sp);
__m256i w1 = LDQ((const __m256i *)(sp + 4));
__m256i w2 = LDQ((const __m256i *)(sp + 8));
__m256i w3 = LDQ((const __m256i *)(sp + 12));
if (o->s[t] < 0) {
v0 = _mm256_sub_epi64(v0, w0); v1 = _mm256_sub_epi64(v1, w1);
v2 = _mm256_sub_epi64(v2, w2); v3 = _mm256_sub_epi64(v3, w3); }
else {
v0 = _mm256_add_epi64(v0, w0); v1 = _mm256_add_epi64(v1, w1);
v2 = _mm256_add_epi64(v2, w2); v3 = _mm256_add_epi64(v3, w3); }
}
STQ((__m256i *)p, v0);
STQ((__m256i *)(p + 4), v1);
STQ((__m256i *)(p + 8), v2);
STQ((__m256i *)(p + 12), v3);
p += 16;
}
for (; k + 4 <= kc; k += 4) {
__m256i v = _mm256_setzero_si256();
for (int t = 0; t < nt; t++) {
__m256i w = LDQ((const __m256i *)(P + (size_t)(o->r[t] + ic + r) * ld + o->c[t] + pc + k));
v = (o->s[t] < 0) ? _mm256_sub_epi64(v, w) : _mm256_add_epi64(v, w);
}
STQ((__m256i *)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);
}
}
}
p += KC - k; /* k >= kc is never read: mikro runs kc steps */
}
}
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;
const int jbn = (nc + 3) & ~3; /* mikro reads only jg < nc, step 4 */
int jb = 0;
/* PACK8: 8 columns per k. The 4-column form loads 32 B at a stride of
ld*8 (4 KB or more), so every load is a distinct, half-used line and the
hardware prefetcher cannot see the stream. Reading jb and jb+4 together
makes each k one FULL 64 B line, halving the line count and the DRAM
traffic. All 8 columns come from one line because jc+jb is a multiple of
8 (jb steps by 8) and P is 4096-byte aligned. Byte-identical: the two
stores write exactly the values the 4-column loop wrote to the same
addresses. */
if (nt <= 2) {
for (; jb + 8 <= nc; jb += 8) {
long long *q0 = Bpanel + (size_t)(jb / 4) * KC * 4;
long long *q1 = q0 + (size_t)KC * 4;
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;
int k = 0;
for (; k + 2 <= kc; k += 2) {
const long long *a0 = sq + (size_t)(pc + k) * ld;
const long long *a1 = a0 + ld;
if (k + 3 < kc) { _mm_prefetch((const char *)(a0 + (size_t)3 * ld), _MM_HINT_T0);
_mm_prefetch((const char *)(a1 + (size_t)3 * ld), _MM_HINT_T0); }
__m256i v0 = LDQ((const __m256i *)a0), u0 = LDQ((const __m256i *)(a0 + 4));
__m256i v1 = LDQ((const __m256i *)a1), u1 = LDQ((const __m256i *)(a1 + 4));
if (nt == 2) {
const long long *b0 = sq1 + (size_t)(pc + k) * ld;
const long long *b1 = b0 + ld;
__m256i w0 = LDQ((const __m256i *)b0), x0 = LDQ((const __m256i *)(b0 + 4));
__m256i w1 = LDQ((const __m256i *)b1), x1 = LDQ((const __m256i *)(b1 + 4));
if (o->s[1] > 0) { v0 = _mm256_add_epi64(v0, w0); u0 = _mm256_add_epi64(u0, x0);
v1 = _mm256_add_epi64(v1, w1); u1 = _mm256_add_epi64(u1, x1); }
else { v0 = _mm256_sub_epi64(v0, w0); u0 = _mm256_sub_epi64(u0, x0);
v1 = _mm256_sub_epi64(v1, w1); u1 = _mm256_sub_epi64(u1, x1); }
}
_mm256_store_si256((__m256i *)(q0 + (size_t)k * 4), SWQ(v0));
_mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), SWQ(u0));
_mm256_store_si256((__m256i *)(q0 + (size_t)(k + 1) * 4), SWQ(v1));
_mm256_store_si256((__m256i *)(q1 + (size_t)(k + 1) * 4), SWQ(u1));
}
for (; k < kc; k++) {
const long long *a = sq + (size_t)(pc + k) * ld;
__m256i v = LDQ((const __m256i *)a);
__m256i u = LDQ((const __m256i *)(a + 4));
if (nt == 2) {
const long long *b = sq1 + (size_t)(pc + k) * ld;
__m256i w = LDQ((const __m256i *)b);
__m256i x = LDQ((const __m256i *)(b + 4));
if (o->s[1] > 0) { v = _mm256_add_epi64(v, w); u = _mm256_add_epi64(u, x); }
else { v = _mm256_sub_epi64(v, w); u = _mm256_sub_epi64(u, x); }
}
_mm256_store_si256((__m256i *)(q0 + (size_t)k * 4), SWQ(v));
_mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), SWQ(u));
}
}
}
for (; jb < jbn; 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 = LDQ((const __m256i *)(sq + (size_t)(pc + k) * ld));
if (nt == 2) {
__m256i w = LDQ((const __m256i *)(sq1 + (size_t)(pc + k) * ld));
v = (o->s[1] > 0) ? _mm256_add_epi64(v, w) : _mm256_sub_epi64(v, w);
}
_mm256_store_si256((__m256i *)(q + (size_t)k * 4), SWQ(v));
}
} else if (full) {
int k = 0;
for (; k + 8 <= kc; k += 8) {
__m256i v0 = _mm256_setzero_si256(), v1 = v0, v2 = v0, v3 = v0;
__m256i v4 = v0, v5 = v0, v6 = v0, v7 = v0;
for (int tt = 0; tt < nt; tt++) {
const long long *sp = P + (size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb;
__m256i w0 = LDQ((const __m256i *)sp);
__m256i w1 = LDQ((const __m256i *)(sp + ld));
__m256i w2 = LDQ((const __m256i *)(sp + 2 * ld));
__m256i w3 = LDQ((const __m256i *)(sp + 3 * ld));
__m256i w4 = LDQ((const __m256i *)(sp + 4 * ld));
__m256i w5 = LDQ((const __m256i *)(sp + 5 * ld));
__m256i w6 = LDQ((const __m256i *)(sp + 6 * ld));
__m256i w7 = LDQ((const __m256i *)(sp + 7 * ld));
if (o->s[tt] < 0) {
v0 = _mm256_sub_epi64(v0, w0); v1 = _mm256_sub_epi64(v1, w1);
v2 = _mm256_sub_epi64(v2, w2); v3 = _mm256_sub_epi64(v3, w3);
v4 = _mm256_sub_epi64(v4, w4); v5 = _mm256_sub_epi64(v5, w5);
v6 = _mm256_sub_epi64(v6, w6); v7 = _mm256_sub_epi64(v7, w7);
} else {
v0 = _mm256_add_epi64(v0, w0); v1 = _mm256_add_epi64(v1, w1);
v2 = _mm256_add_epi64(v2, w2); v3 = _mm256_add_epi64(v3, w3);
v4 = _mm256_add_epi64(v4, w4); v5 = _mm256_add_epi64(v5, w5);
v6 = _mm256_add_epi64(v6, w6); v7 = _mm256_add_epi64(v7, w7);
}
}
_mm256_store_si256((__m256i *)(q + (size_t)k * 4), SWQ(v0));
_mm256_store_si256((__m256i *)(q + (size_t)(k + 1) * 4), SWQ(v1));
_mm256_store_si256((__m256i *)(q + (size_t)(k + 2) * 4), SWQ(v2));
_mm256_store_si256((__m256i *)(q + (size_t)(k + 3) * 4), SWQ(v3));
_mm256_store_si256((__m256i *)(q + (size_t)(k + 4) * 4), SWQ(v4));
_mm256_store_si256((__m256i *)(q + (size_t)(k + 5) * 4), SWQ(v5));
_mm256_store_si256((__m256i *)(q + (size_t)(k + 6) * 4), SWQ(v6));
_mm256_store_si256((__m256i *)(q + (size_t)(k + 7) * 4), SWQ(v7));
}
for (; k < kc; k++) {
__m256i v = _mm256_setzero_si256();
for (int tt = 0; tt < nt; tt++) {
__m256i w = LDQ((const __m256i *)(P + (size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb));
v = (o->s[tt] < 0) ? _mm256_sub_epi64(v, w) : _mm256_add_epi64(v, w);
}
_mm256_store_si256((__m256i *)(q + (size_t)k * 4), SWQ(v));
}
} 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 = LDQ((const __m256i *)t);
_mm256_store_si256((__m256i *)(q + (size_t)k * 4), SWQ(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 * KCP * 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;
}
}
}
/* recursive Strassen: C1 op1 op= Ao*Bo and C2 op2 op= Ao*Bo */
/* Scratch: one h*h product block per recursion level, carved from the tail of the
buffer. Only the levels that actually recurse consume space, so a walk of the
recursion computes exactly the bytes matrix_multiply() hands us below. */
/* ===================== FUSED LEAF-OPERAND PANEL BUILDER =====================
Replaces, at a node of size 2*CUTL (whose 7 children are leaves), the pair
[ matop(node operand) ] + [ packA x7 / packB x7 inside the leaves ]
with ONE pass per side that reads the node operand's four quadrants once and
writes all 7 leaf-operand panels in the mikro's own layout. */
#define QN 64 /* quadrant size == CUTL */
#define PN 7
typedef char pnl_cutl_is_64[(CUTL == 64) ? 1 : -1];
typedef char pnl_mr_is_6[(MR == 6) ? 1 : -1];
typedef char pnl_kc_is_kcp_off[(KCP == KC + 12) ? 1 : -1];
static long long PAN_A_RAW[PN][(size_t)MR*KCP + 32 + 8*1] __attribute__((aligned(64)));
static long long PAN_B_RAW[PN][(size_t)16*KC*4 + 32 + 8*0] __attribute__((aligned(64)));
static long long *const PAN_A[PN] = {
PAN_A_RAW[0]+8*1,
PAN_A_RAW[1]+8*1,
PAN_A_RAW[2]+8*1,
PAN_A_RAW[3]+8*1,
PAN_A_RAW[4]+8*1,
PAN_A_RAW[5]+8*1,
PAN_A_RAW[6]+8*1 };
static long long *const PAN_B[PN] = {
PAN_B_RAW[0]+8*0,
PAN_B_RAW[1]+8*0,
PAN_B_RAW[2]+8*0,
PAN_B_RAW[3]+8*0,
PAN_B_RAW[4]+8*0,
PAN_B_RAW[5]+8*0,
PAN_B_RAW[6]+8*0 };
/* A side: rows i in [0,64), k in [0,64). Each row is 64 elements = 8 whole
64B lines, so the read pattern is the same as matop2's plane read. */
/* ---- HOISTED BUILDERS ------------------------------------------------------------------
`LDQ`/`STQ` carry a "memory" clobber, so gcc may NOT keep any memory-derived value
in a register across them. The shipped builders therefore re-load `o->r[t]`,
`o->c[t]`, `o->s[t]` and re-multiply `o->r[t]*ld` INSIDE the k-loop (verified in the
g++-9 .s: `movslq 20(%rax),%rsi / imulq %r14,%rsi / movl 84(%rax),%edi / js` appear
in the innermost body). Hoisting them into locals is the same loads, the same
stores and the same addresses -- only loop-invariant arithmetic moves out.
nt==2 (this row's dominant arity) is specialised; a hoisted generic path covers
every other arity and sign pattern. Measured on the judge: buildB -6.09 %,
buildA null; both byte-exact at all four sign combinations. */
static void buildA(const Opd *o) {
const int nt = o->nt; const long long *P = o->p; const int ld = o->ld;
if (nt == 2) {
const size_t qld = (size_t)QN * ld;
const long long *B0 = P + (size_t)o->r[0] * ld + o->c[0];
const long long *B1 = P + (size_t)o->r[1] * ld + o->c[1];
const int n0 = (o->s[0] < 0), n1 = (o->s[1] < 0);
const long long *pA = n0 ? B1 : B0, *pB = n0 ? B0 : B1;
const int sub = (n0 != n1), negall = (n0 && n1);
for (int i = 0; i < QN; i++) {
const size_t pi = (size_t)(i % MR);
long long *p0 = PAN_A[0] + pi * KCP, *p1 = PAN_A[1] + pi * KCP,
*p2 = PAN_A[2] + pi * KCP, *p3 = PAN_A[3] + pi * KCP,
*p4 = PAN_A[4] + pi * KCP, *p5 = PAN_A[5] + pi * KCP,
*p6 = PAN_A[6] + pi * KCP;
const long long *f0 = pA + (size_t)i * ld, *f1 = pB + (size_t)i * ld;
__m256i a[2][4], w, z = _mm256_setzero_si256();
for (int k = 0; k < QN; k += 8) {
/* BOTH 4-blocks' loads issue before either block's stores: the shipped loop
cannot do this (the "memory" clobber forbids gcc the reorder). */
for (int b = 0; b < 2; b++) {
const int kk = k + b * 4;
a[b][0] = LDQ((const __m256i *)(f0 + kk));
w = LDQ((const __m256i *)(f1 + kk));
a[b][0] = sub ? _mm256_sub_epi64(a[b][0], w) : _mm256_add_epi64(a[b][0], w);
a[b][1] = LDQ((const __m256i *)(f0 + QN + kk));
w = LDQ((const __m256i *)(f1 + QN + kk));
a[b][1] = sub ? _mm256_sub_epi64(a[b][1], w) : _mm256_add_epi64(a[b][1], w);
a[b][2] = LDQ((const __m256i *)(f0 + qld + kk));
w = LDQ((const __m256i *)(f1 + qld + kk));
a[b][2] = sub ? _mm256_sub_epi64(a[b][2], w) : _mm256_add_epi64(a[b][2], w);
a[b][3] = LDQ((const __m256i *)(f0 + qld + QN + kk));
w = LDQ((const __m256i *)(f1 + qld + QN + kk));
a[b][3] = sub ? _mm256_sub_epi64(a[b][3], w) : _mm256_add_epi64(a[b][3], w);
if (negall) { for (int q = 0; q < 4; q++) a[b][q] = _mm256_sub_epi64(z, a[b][q]); }
}
for (int b = 0; b < 2; b++) {
const int kk = k + b * 4;
STQ((__m256i *)(p0 + kk), _mm256_add_epi64(a[b][0], a[b][3])); /* M1 A11+A22 */
STQ((__m256i *)(p1 + kk), _mm256_add_epi64(a[b][2], a[b][3])); /* M2 A21+A22 */
STQ((__m256i *)(p2 + kk), a[b][0]); /* M3 A11 */
STQ((__m256i *)(p3 + kk), a[b][3]); /* M4 A22 */
STQ((__m256i *)(p4 + kk), _mm256_add_epi64(a[b][0], a[b][1])); /* M5 A11+A12 */
STQ((__m256i *)(p5 + kk), _mm256_sub_epi64(a[b][2], a[b][0])); /* M6 A21-A11 */
STQ((__m256i *)(p6 + kk), _mm256_sub_epi64(a[b][1], a[b][3])); /* M7 A12-A22 */
}
}
}
return;
}
{ /* generic hoisted path: identical arithmetic, pointers and signs hoisted */
const long long *Bp[8]; int sg[8];
for (int t = 0; t < nt; t++) { Bp[t] = P + (size_t)o->r[t] * ld + o->c[t]; sg[t] = (o->s[t] < 0); }
const size_t qld = (size_t)QN * ld;
for (int i = 0; i < QN; i++) {
const size_t pi = (size_t)(i % MR);
long long *p0 = PAN_A[0] + pi * KCP, *p1 = PAN_A[1] + pi * KCP,
*p2 = PAN_A[2] + pi * KCP, *p3 = PAN_A[3] + pi * KCP,
*p4 = PAN_A[4] + pi * KCP, *p5 = PAN_A[5] + pi * KCP,
*p6 = PAN_A[6] + pi * KCP;
for (int k = 0; k < QN; k += 4) {
__m256i a11, a12, a21, a22;
if (nt == 1 && !sg[0]) {
a11 = LDQ((const __m256i *)(Bp[0] + (size_t)i * ld + k));
a12 = LDQ((const __m256i *)(Bp[0] + (size_t)i * ld + QN + k));
a21 = LDQ((const __m256i *)(Bp[0] + (size_t)(i + QN) * ld + k));
a22 = LDQ((const __m256i *)(Bp[0] + (size_t)(i + QN) * ld + QN + k));
} else {
a11 = _mm256_setzero_si256(); a12 = a11; a21 = a11; a22 = a11;
for (int t = 0; t < nt; t++) {
const long long *b = Bp[t] + (size_t)i * ld;
const int neg = sg[t]; __m256i w;
w = LDQ((const __m256i *)(b + k));
a11 = neg ? _mm256_sub_epi64(a11, w) : _mm256_add_epi64(a11, w);
w = LDQ((const __m256i *)(b + QN + k));
a12 = neg ? _mm256_sub_epi64(a12, w) : _mm256_add_epi64(a12, w);
w = LDQ((const __m256i *)(b + (size_t)QN * ld + k));
a21 = neg ? _mm256_sub_epi64(a21, w) : _mm256_add_epi64(a21, w);
w = LDQ((const __m256i *)(b + (size_t)QN * ld + QN + k));
a22 = neg ? _mm256_sub_epi64(a22, w) : _mm256_add_epi64(a22, w);
}
}
STQ((__m256i *)(p0 + k), _mm256_add_epi64(a11, a22)); /* M1 A11+A22 */
STQ((__m256i *)(p1 + k), _mm256_add_epi64(a21, a22)); /* M2 A21+A22 */
STQ((__m256i *)(p2 + k), a11); /* M3 A11 */
STQ((__m256i *)(p3 + k), a22); /* M4 A22 */
STQ((__m256i *)(p4 + k), _mm256_add_epi64(a11, a12)); /* M5 A11+A12 */
STQ((__m256i *)(p5 + k), _mm256_sub_epi64(a21, a11)); /* M6 A21-A11 */
STQ((__m256i *)(p6 + k), _mm256_sub_epi64(a12, a22)); /* M7 A12-A22 */
}
}
}
}
static void buildB(const Opd *o) {
const int nt = o->nt; const long long *P = o->p; const int ld = o->ld;
const size_t qld = (size_t)QN * ld;
const long long *Bp[8]; int sg[8]; int allpos = 1;
for (int t = 0; t < nt; t++) {
Bp[t] = P + (size_t)o->r[t] * ld + o->c[t]; sg[t] = (o->s[t] < 0);
if (sg[t]) allpos = 0;
}
for (int jb = 0; jb < QN; jb += 8) {
const size_t g = (size_t)(jb >> 2) * KC * 4;
long long *p0 = PAN_B[0] + g, *p1 = PAN_B[1] + g, *p2 = PAN_B[2] + g,
*p3 = PAN_B[3] + g, *p4 = PAN_B[4] + g, *p5 = PAN_B[5] + g,
*p6 = PAN_B[6] + g;
for (int k = 0; k < QN; k++) {
__m256i b11a, b11b, b12a, b12b, b21a, b21b, b22a, b22b;
if (nt == 1 && allpos) {
const long long *r0p = Bp[0] + (size_t)k * ld;
const long long *r1p = r0p + qld;
b11a = LDQ((const __m256i *)(r0p + jb)); b11b = LDQ((const __m256i *)(r0p + jb + 4));
b12a = LDQ((const __m256i *)(r0p + QN + jb)); b12b = LDQ((const __m256i *)(r0p + QN + jb + 4));
b21a = LDQ((const __m256i *)(r1p + jb)); b21b = LDQ((const __m256i *)(r1p + jb + 4));
b22a = LDQ((const __m256i *)(r1p + QN + jb)); b22b = LDQ((const __m256i *)(r1p + QN + jb + 4));
} else if (nt == 2 && allpos) {
/* pairwise combine straight off the two base pointers: the shipped generic
loop re-derives Bp[0]/Bp[1] and re-tests the signs every (k,t). */
const long long *r0p = Bp[0] + (size_t)k * ld, *q0p = Bp[1] + (size_t)k * ld;
const long long *r1p = r0p + qld, *q1p = q0p + qld;
__m256i w, x;
b11a = LDQ((const __m256i *)(r0p + jb)); x = LDQ((const __m256i *)(q0p + jb));
b11b = LDQ((const __m256i *)(r0p + jb + 4)); w = LDQ((const __m256i *)(q0p + jb + 4));
b11a = _mm256_add_epi64(b11a, x); b11b = _mm256_add_epi64(b11b, w);
b12a = LDQ((const __m256i *)(r0p + QN + jb)); x = LDQ((const __m256i *)(q0p + QN + jb));
b12b = LDQ((const __m256i *)(r0p + QN + jb + 4)); w = LDQ((const __m256i *)(q0p + QN + jb + 4));
b12a = _mm256_add_epi64(b12a, x); b12b = _mm256_add_epi64(b12b, w);
b21a = LDQ((const __m256i *)(r1p + jb)); x = LDQ((const __m256i *)(q1p + jb));
b21b = LDQ((const __m256i *)(r1p + jb + 4)); w = LDQ((const __m256i *)(q1p + jb + 4));
b21a = _mm256_add_epi64(b21a, x); b21b = _mm256_add_epi64(b21b, w);
b22a = LDQ((const __m256i *)(r1p + QN + jb)); x = LDQ((const __m256i *)(q1p + QN + jb));
b22b = LDQ((const __m256i *)(r1p + QN + jb + 4)); w = LDQ((const __m256i *)(q1p + QN + jb + 4));
b22a = _mm256_add_epi64(b22a, x); b22b = _mm256_add_epi64(b22b, w);
} else {
b11a = b11b = b12a = b12b = b21a = b21b = b22a = b22b = _mm256_setzero_si256();
for (int t = 0; t < nt; t++) {
const long long *r0p = Bp[t] + (size_t)k * ld;
const long long *r1p = r0p + qld;
const int neg = sg[t];
__m256i w;
w = LDQ((const __m256i *)(r0p + jb)); b11a = neg ? _mm256_sub_epi64(b11a, w) : _mm256_add_epi64(b11a, w);
w = LDQ((const __m256i *)(r0p + jb + 4)); b11b = neg ? _mm256_sub_epi64(b11b, w) : _mm256_add_epi64(b11b, w);
w = LDQ((const __m256i *)(r0p + QN + jb)); b12a = neg ? _mm256_sub_epi64(b12a, w) : _mm256_add_epi64(b12a, w);
w = LDQ((const __m256i *)(r0p + QN + jb + 4));b12b = neg ? _mm256_sub_epi64(b12b, w) : _mm256_add_epi64(b12b, w);
w = LDQ((const __m256i *)(r1p + jb)); b21a = neg ? _mm256_sub_epi64(b21a, w) : _mm256_add_epi64(b21a, w);
w = LDQ((const __m256i *)(r1p + jb + 4)); b21b = neg ? _mm256_sub_epi64(b21b, w) : _mm256_add_epi64(b21b, w);
w = LDQ((const __m256i *)(r1p + QN + jb)); b22a = neg ? _mm256_sub_epi64(b22a, w) : _mm256_add_epi64(b22a, w);
w = LDQ((const __m256i *)(r1p + QN + jb + 4));b22b = neg ? _mm256_sub_epi64(b22b, w) : _mm256_add_epi64(b22b, w);
}
}
const size_t ka = (size_t)k * 4, kb = ka + (size_t)KC * 4;
STQ((__m256i *)(p0 + ka), SWQ(_mm256_add_epi64(b11a, b22a)));
STQ((__m256i *)(p0 + kb), SWQ(_mm256_add_epi64(b11b, b22b)));
STQ((__m256i *)(p1 + ka), SWQ(b11a));
STQ((__m256i *)(p1 + kb), SWQ(b11b));
STQ((__m256i *)(p2 + ka), SWQ(_mm256_sub_epi64(b12a, b22a)));
STQ((__m256i *)(p2 + kb), SWQ(_mm256_sub_epi64(b12b, b22b)));
STQ((__m256i *)(p3 + ka), SWQ(_mm256_sub_epi64(b21a, b11a)));
STQ((__m256i *)(p3 + kb), SWQ(_mm256_sub_epi64(b21b, b11b)));
STQ((__m256i *)(p4 + ka), SWQ(b22a));
STQ((__m256i *)(p4 + kb), SWQ(b22b));
STQ((__m256i *)(p5 + ka), SWQ(_mm256_add_epi64(b11a, b12a)));
STQ((__m256i *)(p5 + kb), SWQ(_mm256_add_epi64(b11b, b12b)));
STQ((__m256i *)(p6 + ka), SWQ(_mm256_add_epi64(b21a, b22a)));
STQ((__m256i *)(p6 + kb), SWQ(_mm256_add_epi64(b21b, b22b)));
}
}
}
/* the leaf, driven from prebuilt panels. Loop nest, order, mikro call
sequence and destination addressing are the shipped leaf's; packA / packB
are replaced by the panel pointers. */
static void leafP(int n, int L, long long *C1, int ldc1, int op1,
long long *C2, int ldc2, int op2) {
const long long *pa = PAN_A[L];
const long long *pb = PAN_B[L];
for (int ic = 0; ic < n; ic += MC) {
int mc = (n - ic < MC) ? (n - ic) : MC;
for (int ir = 0; ir < mc; ir += MR) {
int rows = (mc - ir < MR) ? (mc - ir) : MR;
const long long *ap = pa + (size_t)(ic + ir) * KCP;
long long *c1p = C1 + (size_t)(ic + ir) * ldc1;
long long *c2p = C2 ? C2 + (size_t)(ic + ir) * ldc2 : 0;
for (int jg = 0; jg < n; jg += 4) {
mikro(n, ap, pb + (size_t)(jg >> 2) * KC * 4, c1p + jg, ldc1, op1,
c2p ? c2p + jg : 0, ldc2, op2, rows);
}
}
}
}
static void leafP_blk(int n, int L, int ic0, int ic1, long long *C1, int ldc1, int op1,
long long *C2, int ldc2, int op2) {
const long long *pa = PAN_A[L];
const long long *pb = PAN_B[L];
for (int ic = ic0; ic < ic1; ic += MC) {
int mc = (n - ic < MC) ? (n - ic) : MC;
if (mc > ic1 - ic) mc = ic1 - ic;
for (int ir = 0; ir < mc; ir += MR) {
int rows = (mc - ir < MR) ? (mc - ir) : MR;
const long long *ap = pa + (size_t)(ic + ir - ic0) * KCP;
long long *c1p = C1 + (size_t)(ic + ir) * ldc1;
long long *c2p = C2 ? C2 + (size_t)(ic + ir) * ldc2 : 0;
for (int jg = 0; jg < n; jg += 4) {
mikro(n, ap, pb + (size_t)(jg >> 2) * KC * 4, c1p + jg, ldc1, op1,
c2p ? c2p + jg : 0, ldc2, op2, rows);
}
}
}
}
static long long g_scratch[6u * 1024u * 1024u] __attribute__((aligned(64)));
/* dst op= sign * T (op: 0 store, 1 add, 2 sub) */
static void cmb(long long *dst, int ldd, const long long *T, int ldt, int h, int sign, int op) {
__m256i sgn = _mm256_set1_epi64x(sign < 0 ? -1 : 0);
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;
/* ARM cmbu16: the j-loop's inner step raised 4 -> 16 elements so FOUR loads
are issued before the first store. `LDQ` carries a "memory" clobber, so
the shipped 4-element step cannot be software-pipelined by gcc; here the
pipelining is explicit. Same loads, same stores, same addresses. */
if (op == 0) {
for (; j + 16 <= h; j += 16) {
__m256i v0 = LDQ((const __m256i *)(t + j)), v1 = LDQ((const __m256i *)(t + j + 4));
__m256i v2 = LDQ((const __m256i *)(t + j + 8)), v3 = LDQ((const __m256i *)(t + j + 12));
if (sign < 0) { __m256i z = _mm256_setzero_si256();
v0 = _mm256_sub_epi64(z, v0); v1 = _mm256_sub_epi64(z, v1);
v2 = _mm256_sub_epi64(z, v2); v3 = _mm256_sub_epi64(z, v3); }
STQ((__m256i *)(c + j), v0); STQ((__m256i *)(c + j + 4), v1);
STQ((__m256i *)(c + j + 8), v2); STQ((__m256i *)(c + j + 12), v3);
}
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 + 16 <= h; j += 16) {
__m256i v0 = LDQ((const __m256i *)(t + j)), v1 = LDQ((const __m256i *)(t + j + 4));
__m256i v2 = LDQ((const __m256i *)(t + j + 8)), v3 = LDQ((const __m256i *)(t + j + 12));
__m256i w0 = LDQ((const __m256i *)(c + j)), w1 = LDQ((const __m256i *)(c + j + 4));
__m256i w2 = LDQ((const __m256i *)(c + j + 8)), w3 = LDQ((const __m256i *)(c + j + 12));
if ((op == 1) != (sign < 0)) {
STQ((__m256i *)(c + j), _mm256_add_epi64(w0, v0));
STQ((__m256i *)(c + j + 4), _mm256_add_epi64(w1, v1));
STQ((__m256i *)(c + j + 8), _mm256_add_epi64(w2, v2));
STQ((__m256i *)(c + j + 12), _mm256_add_epi64(w3, v3));
} else {
STQ((__m256i *)(c + j), _mm256_sub_epi64(w0, v0));
STQ((__m256i *)(c + j + 4), _mm256_sub_epi64(w1, v1));
STQ((__m256i *)(c + j + 8), _mm256_sub_epi64(w2, v2));
STQ((__m256i *)(c + j + 12), _mm256_sub_epi64(w3, v3));
}
}
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]; }
}
}
}
/* fused two-destination combine: reads T ONCE for both destinations.
Semantics identical to cmb(d1,s1,o1) then cmb(d2,s2,o2) on the same T. */
static void cmb2(long long *d1,int ld1,int op1,int s1,
long long *d2,int ld2,int op2,int s2,
const long long *T,int ldt,int h) {
for (int i = 0; i < h; i++) {
const long long *t = T + (size_t)i * ldt;
long long *c1 = d1 + (size_t)i * ld1, *c2 = d2 + (size_t)i * ld2;
int j = 0;
for (; j + 16 <= h; j += 16) {
__m256i v0=LDQ((const __m256i*)(t+j)), v1=LDQ((const __m256i*)(t+j+4));
__m256i v2=LDQ((const __m256i*)(t+j+8)), v3=LDQ((const __m256i*)(t+j+12));
__m256i a0=LDQ((const __m256i*)(c1+j)), a1=LDQ((const __m256i*)(c1+j+4));
__m256i a2=LDQ((const __m256i*)(c1+j+8)), a3=LDQ((const __m256i*)(c1+j+12));
__m256i b0=LDQ((const __m256i*)(c2+j)), b1=LDQ((const __m256i*)(c2+j+4));
__m256i b2=LDQ((const __m256i*)(c2+j+8)), b3=LDQ((const __m256i*)(c2+j+12));
if ((op1 == 1) != (s1 < 0)) { STQ((__m256i*)(c1+j),_mm256_add_epi64(a0,v0));
STQ((__m256i*)(c1+j+4),_mm256_add_epi64(a1,v1)); STQ((__m256i*)(c1+j+8),_mm256_add_epi64(a2,v2));
STQ((__m256i*)(c1+j+12),_mm256_add_epi64(a3,v3)); }
else { STQ((__m256i*)(c1+j),_mm256_sub_epi64(a0,v0));
STQ((__m256i*)(c1+j+4),_mm256_sub_epi64(a1,v1)); STQ((__m256i*)(c1+j+8),_mm256_sub_epi64(a2,v2));
STQ((__m256i*)(c1+j+12),_mm256_sub_epi64(a3,v3)); }
if ((op2 == 1) != (s2 < 0)) { STQ((__m256i*)(c2+j),_mm256_add_epi64(b0,v0));
STQ((__m256i*)(c2+j+4),_mm256_add_epi64(b1,v1)); STQ((__m256i*)(c2+j+8),_mm256_add_epi64(b2,v2));
STQ((__m256i*)(c2+j+12),_mm256_add_epi64(b3,v3)); }
else { STQ((__m256i*)(c2+j),_mm256_sub_epi64(b0,v0));
STQ((__m256i*)(c2+j+4),_mm256_sub_epi64(b1,v1)); STQ((__m256i*)(c2+j+8),_mm256_sub_epi64(b2,v2));
STQ((__m256i*)(c2+j+12),_mm256_sub_epi64(b3,v3)); }
}
for (; j + 4 <= h; j += 4) {
__m256i v = LDQ((const __m256i *)(t + j));
__m256i a = LDQ((const __m256i *)(c1 + j));
STQ((__m256i *)(c1 + j), ((op1 == 1) != (s1 < 0)) ? _mm256_add_epi64(a,v)
: _mm256_sub_epi64(a,v));
__m256i b = LDQ((const __m256i *)(c2 + j));
STQ((__m256i *)(c2 + j), ((op2 == 1) != (s2 < 0)) ? _mm256_add_epi64(b,v)
: _mm256_sub_epi64(b,v));
}
for (; j < h; j++) {
long long v = t[j];
if ((op1 == 1) != (s1 < 0)) c1[j] += v; else c1[j] -= v;
if ((op2 == 1) != (s2 < 0)) c2[j] += v; else c2[j] -= v;
}
}
}
/* l1k15_rv `fus3`: the three products whose SECOND destination is C22 (M1, M2, M3)
land in three quadrants that are ALREADY LIVE (C11, C21, C12) -- so the three
sequential combine passes over C22 collapse into ONE pass. Source bytes are
identical; the destination is written ONCE instead of three times.
Semantics == cmb(C22,C11,+1); cmb(C22,C21,-1); cmb(C22,C12,+1) [op==0] */
static void cmb3_p(long long *d, int ldd,
const long long *a, const long long *b, const long long *c, int ld,
int rows, int ncol) {
for (int i = 0; i < rows; i++) {
const long long *pa = a + (size_t)i * ld, *pb = b + (size_t)i * ld, *pc = c + (size_t)i * ld;
long long *pd = d + (size_t)i * ldd;
int j = 0;
for (; j + 16 <= ncol; j += 16) {
__m256i A0=LDQ((const __m256i*)(pa+j)), A1=LDQ((const __m256i*)(pa+j+4));
__m256i A2=LDQ((const __m256i*)(pa+j+8)), A3=LDQ((const __m256i*)(pa+j+12));
__m256i B0=LDQ((const __m256i*)(pb+j)), B1=LDQ((const __m256i*)(pb+j+4));
__m256i B2=LDQ((const __m256i*)(pb+j+8)), B3=LDQ((const __m256i*)(pb+j+12));
__m256i C0=LDQ((const __m256i*)(pc+j)), C1=LDQ((const __m256i*)(pc+j+4));
__m256i C2=LDQ((const __m256i*)(pc+j+8)), C3=LDQ((const __m256i*)(pc+j+12));
STQ((__m256i*)(pd+j), _mm256_sub_epi64(_mm256_add_epi64(A0,C0),B0));
STQ((__m256i*)(pd+j+4), _mm256_sub_epi64(_mm256_add_epi64(A1,C1),B1));
STQ((__m256i*)(pd+j+8), _mm256_sub_epi64(_mm256_add_epi64(A2,C2),B2));
STQ((__m256i*)(pd+j+12), _mm256_sub_epi64(_mm256_add_epi64(A3,C3),B3));
}
for (; j + 4 <= ncol; j += 4)
STQ((__m256i*)(pd+j), _mm256_sub_epi64(
_mm256_add_epi64(LDQ((const __m256i*)(pa+j)), LDQ((const __m256i*)(pc+j))),
LDQ((const __m256i*)(pb+j))));
for (; j < ncol; j++) pd[j] = pa[j] - pb[j] + pc[j];
}
}
static void cmb3_ntx(long long *d, int ldd,
const long long *a, const long long *b, const long long *c, int ld,
int rows, int ncol) {
for (int i = 0; i < rows; i++) {
const long long *pa = a + (size_t)i * ld, *pb = b + (size_t)i * ld, *pc = c + (size_t)i * ld;
long long *pd = d + (size_t)i * ldd;
int j = 0;
for (; j + 16 <= ncol; j += 16) {
__m256i A0=LDQ((const __m256i*)(pa+j)), A1=LDQ((const __m256i*)(pa+j+4));
__m256i A2=LDQ((const __m256i*)(pa+j+8)), A3=LDQ((const __m256i*)(pa+j+12));
__m256i B0=LDQ((const __m256i*)(pb+j)), B1=LDQ((const __m256i*)(pb+j+4));
__m256i B2=LDQ((const __m256i*)(pb+j+8)), B3=LDQ((const __m256i*)(pb+j+12));
__m256i C0=LDQ((const __m256i*)(pc+j)), C1=LDQ((const __m256i*)(pc+j+4));
__m256i C2=LDQ((const __m256i*)(pc+j+8)), C3=LDQ((const __m256i*)(pc+j+12));
STQNT((__m256i*)(pd+j), _mm256_sub_epi64(_mm256_add_epi64(A0,C0),B0));
STQNT((__m256i*)(pd+j+4), _mm256_sub_epi64(_mm256_add_epi64(A1,C1),B1));
STQNT((__m256i*)(pd+j+8), _mm256_sub_epi64(_mm256_add_epi64(A2,C2),B2));
STQNT((__m256i*)(pd+j+12), _mm256_sub_epi64(_mm256_add_epi64(A3,C3),B3));
}
for (; j + 4 <= ncol; j += 4)
STQNT((__m256i*)(pd+j), _mm256_sub_epi64(
_mm256_add_epi64(LDQ((const __m256i*)(pa+j)), LDQ((const __m256i*)(pc+j))),
LDQ((const __m256i*)(pb+j))));
for (; j < ncol; j++) pd[j] = pa[j] - pb[j] + pc[j];
}
}
static void cmb3_blk(long long *d, int ldd, const long long *a, const long long *b,
const long long *c, int ld, int h, int i0, int i1) {
if (i1 <= i0) return;
cmb3(d + (size_t)i0 * ldd, ldd, a + (size_t)i0 * ld, b + (size_t)i0 * ld,
c + (size_t)i0 * ld, ld, i1 - i0, h);
}
/* single-destination recursive Strassen: d op= Ao*Bo, with T a scratch block of h*h */
/* materialise one level's operands ONCE per parent call, so the LEAF's packing reads 2
contiguous sources per element instead of 4 strided ones. Only applied at the level whose
children are leaves (n == 2*CUTL): the cost is 2*n^2 element-writes per parent, and it is
amortised over that parent's 7 leaf calls, which read 4x more elements than it writes. */
static long long matbufA[3][262144 + 64], matbufB[3][262144 + 64];
/* lane m142: THE FUSION'S MAXIMUM FORM, built as a measurement.
`matop`'s t-loop over nt is a RUNTIME loop with a sign branch per term; the dominant
arity on this row is nt==2 (m101's firing table: "by arity: all nt=2"). A FULLY
UNROLLED nt==2 path removes exactly the work a "fuse the materialisation into the
combine" construction could remove -- the per-element CONTROL -- while leaving every
byte of traffic identical. If THIS does not pay, no pass-fusion can: the fusion
changes no traffic either (2 dense reads + 1 dense write, either way). Correct by
construction == the generic path's arithmetic. */
static void matop2_p(const Opd *o, int n, long long *dst) {
const long long *P = o->p; const int ld = o->ld;
const long long *p0 = P + (size_t)o->r[0] * ld + o->c[0];
const long long *p1 = P + (size_t)o->r[1] * ld + o->c[1];
const int s1neg = (o->s[1] < 0);
int i = 0;
/* TEST ARM: the i-loop unrolled by 2. Two destination rows are in flight,
so the row-jump latency that serialises the shipped loop is exposed once
per TWO rows instead of once per row. Every load, store and address is
the same as the shipped function's; only the ORDER changes.
The unrolled path is taken only when n % 16 == 0, because the shipped
tail loop writes up to 3 elements PAST the row end and rows are written
in order there; the real call sizes are 128/256/512 and are all % 16 == 0. */
if ((n & 15) == 0) for (; i + 2 <= n; i += 2) {
long long *dpa = dst + (size_t)i * n;
long long *dpb = dpa + n;
const long long *a = p0 + (size_t)i * ld, *b = p1 + (size_t)i * ld;
const long long *c = a + ld, *e = b + ld;
_mm_prefetch((const char *)(p0 + (size_t)(i + 4) * ld), _MM_HINT_T0);
_mm_prefetch((const char *)(p1 + (size_t)(i + 4) * ld), _MM_HINT_T0);
int j = 0;
if (!s1neg) {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)), w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
__m256i x0=LDQ((const __m256i *)(c+j)), x1=LDQ((const __m256i *)(c+j+4));
__m256i x2=LDQ((const __m256i *)(c+j+8)), x3=LDQ((const __m256i *)(c+j+12));
__m256i y0=LDQ((const __m256i *)(b+j)), y1=LDQ((const __m256i *)(b+j+4));
__m256i y2=LDQ((const __m256i *)(b+j+8)), y3=LDQ((const __m256i *)(b+j+12));
__m256i z0=LDQ((const __m256i *)(e+j)), z1=LDQ((const __m256i *)(e+j+4));
__m256i z2=LDQ((const __m256i *)(e+j+8)), z3=LDQ((const __m256i *)(e+j+12));
STQ((__m256i *)(dpa+j), _mm256_add_epi64(w0,y0));
STQ((__m256i *)(dpa+j+4), _mm256_add_epi64(w1,y1));
STQ((__m256i *)(dpa+j+8), _mm256_add_epi64(w2,y2));
STQ((__m256i *)(dpa+j+12), _mm256_add_epi64(w3,y3));
STQ((__m256i *)(dpb+j), _mm256_add_epi64(x0,z0));
STQ((__m256i *)(dpb+j+4), _mm256_add_epi64(x1,z1));
STQ((__m256i *)(dpb+j+8), _mm256_add_epi64(x2,z2));
STQ((__m256i *)(dpb+j+12), _mm256_add_epi64(x3,z3));
}
for (; j < n; j += 4) {
STQ((__m256i *)(dpa+j), _mm256_add_epi64(LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
STQ((__m256i *)(dpb+j), _mm256_add_epi64(LDQ((const __m256i *)(c+j)), LDQ((const __m256i *)(e+j))));
}
} else {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)), w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
__m256i x0=LDQ((const __m256i *)(c+j)), x1=LDQ((const __m256i *)(c+j+4));
__m256i x2=LDQ((const __m256i *)(c+j+8)), x3=LDQ((const __m256i *)(c+j+12));
__m256i y0=LDQ((const __m256i *)(b+j)), y1=LDQ((const __m256i *)(b+j+4));
__m256i y2=LDQ((const __m256i *)(b+j+8)), y3=LDQ((const __m256i *)(b+j+12));
__m256i z0=LDQ((const __m256i *)(e+j)), z1=LDQ((const __m256i *)(e+j+4));
__m256i z2=LDQ((const __m256i *)(e+j+8)), z3=LDQ((const __m256i *)(e+j+12));
STQ((__m256i *)(dpa+j), _mm256_sub_epi64(w0,y0));
STQ((__m256i *)(dpa+j+4), _mm256_sub_epi64(w1,y1));
STQ((__m256i *)(dpa+j+8), _mm256_sub_epi64(w2,y2));
STQ((__m256i *)(dpa+j+12), _mm256_sub_epi64(w3,y3));
STQ((__m256i *)(dpb+j), _mm256_sub_epi64(x0,z0));
STQ((__m256i *)(dpb+j+4), _mm256_sub_epi64(x1,z1));
STQ((__m256i *)(dpb+j+8), _mm256_sub_epi64(x2,z2));
STQ((__m256i *)(dpb+j+12), _mm256_sub_epi64(x3,z3));
}
for (; j < n; j += 4) {
STQ((__m256i *)(dpa+j), _mm256_sub_epi64(LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
STQ((__m256i *)(dpb+j), _mm256_sub_epi64(LDQ((const __m256i *)(c+j)), LDQ((const __m256i *)(e+j))));
}
}
}
for (; i < n; i++) {
long long *dp = dst + (size_t)i * n;
const long long *a = p0 + (size_t)i * ld, *b = p1 + (size_t)i * ld;
_mm_prefetch((const char *)(p0 + (size_t)(i + 2) * ld), _MM_HINT_T0);
_mm_prefetch((const char *)(p1 + (size_t)(i + 2) * ld), _MM_HINT_T0);
int j = 0;
if (!s1neg) {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)), w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
w0=_mm256_add_epi64(w0,LDQ((const __m256i *)(b+j)));
w1=_mm256_add_epi64(w1,LDQ((const __m256i *)(b+j+4)));
w2=_mm256_add_epi64(w2,LDQ((const __m256i *)(b+j+8)));
w3=_mm256_add_epi64(w3,LDQ((const __m256i *)(b+j+12)));
STQ((__m256i *)(dp+j), w0); STQ((__m256i *)(dp+j+4), w1);
STQ((__m256i *)(dp+j+8), w2); STQ((__m256i *)(dp+j+12), w3);
}
for (; j < n; j += 4) STQ((__m256i *)(dp+j), _mm256_add_epi64(
LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
} else {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)), w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
w0=_mm256_sub_epi64(w0,LDQ((const __m256i *)(b+j)));
w1=_mm256_sub_epi64(w1,LDQ((const __m256i *)(b+j+4)));
w2=_mm256_sub_epi64(w2,LDQ((const __m256i *)(b+j+8)));
w3=_mm256_sub_epi64(w3,LDQ((const __m256i *)(b+j+12)));
STQ((__m256i *)(dp+j), w0); STQ((__m256i *)(dp+j+4), w1);
STQ((__m256i *)(dp+j+8), w2); STQ((__m256i *)(dp+j+12), w3);
}
for (; j < n; j += 4) STQ((__m256i *)(dp+j), _mm256_sub_epi64(
LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
}
}
}
static void matop2_ntx(const Opd *o, int n, long long *dst) {
const long long *P = o->p; const int ld = o->ld;
const long long *p0 = P + (size_t)o->r[0] * ld + o->c[0];
const long long *p1 = P + (size_t)o->r[1] * ld + o->c[1];
const int s1neg = (o->s[1] < 0);
int i = 0;
/* TEST ARM: the i-loop unrolled by 2. Two destination rows are in flight,
so the row-jump latency that serialises the shipped loop is exposed once
per TWO rows instead of once per row. Every load, store and address is
the same as the shipped function's; only the ORDER changes.
The unrolled path is taken only when n % 16 == 0, because the shipped
tail loop writes up to 3 elements PAST the row end and rows are written
in order there; the real call sizes are 128/256/512 and are all % 16 == 0. */
if ((n & 15) == 0) for (; i + 2 <= n; i += 2) {
long long *dpa = dst + (size_t)i * n;
long long *dpb = dpa + n;
const long long *a = p0 + (size_t)i * ld, *b = p1 + (size_t)i * ld;
const long long *c = a + ld, *e = b + ld;
_mm_prefetch((const char *)(p0 + (size_t)(i + 4) * ld), _MM_HINT_T0);
_mm_prefetch((const char *)(p1 + (size_t)(i + 4) * ld), _MM_HINT_T0);
int j = 0;
if (!s1neg) {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)), w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
__m256i x0=LDQ((const __m256i *)(c+j)), x1=LDQ((const __m256i *)(c+j+4));
__m256i x2=LDQ((const __m256i *)(c+j+8)), x3=LDQ((const __m256i *)(c+j+12));
__m256i y0=LDQ((const __m256i *)(b+j)), y1=LDQ((const __m256i *)(b+j+4));
__m256i y2=LDQ((const __m256i *)(b+j+8)), y3=LDQ((const __m256i *)(b+j+12));
__m256i z0=LDQ((const __m256i *)(e+j)), z1=LDQ((const __m256i *)(e+j+4));
__m256i z2=LDQ((const __m256i *)(e+j+8)), z3=LDQ((const __m256i *)(e+j+12));
STQNT((__m256i *)(dpa+j), _mm256_add_epi64(w0,y0));
STQNT((__m256i *)(dpa+j+4), _mm256_add_epi64(w1,y1));
STQNT((__m256i *)(dpa+j+8), _mm256_add_epi64(w2,y2));
STQNT((__m256i *)(dpa+j+12), _mm256_add_epi64(w3,y3));
STQNT((__m256i *)(dpb+j), _mm256_add_epi64(x0,z0));
STQNT((__m256i *)(dpb+j+4), _mm256_add_epi64(x1,z1));
STQNT((__m256i *)(dpb+j+8), _mm256_add_epi64(x2,z2));
STQNT((__m256i *)(dpb+j+12), _mm256_add_epi64(x3,z3));
}
for (; j < n; j += 4) {
STQNT((__m256i *)(dpa+j), _mm256_add_epi64(LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
STQNT((__m256i *)(dpb+j), _mm256_add_epi64(LDQ((const __m256i *)(c+j)), LDQ((const __m256i *)(e+j))));
}
} else {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)), w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
__m256i x0=LDQ((const __m256i *)(c+j)), x1=LDQ((const __m256i *)(c+j+4));
__m256i x2=LDQ((const __m256i *)(c+j+8)), x3=LDQ((const __m256i *)(c+j+12));
__m256i y0=LDQ((const __m256i *)(b+j)), y1=LDQ((const __m256i *)(b+j+4));
__m256i y2=LDQ((const __m256i *)(b+j+8)), y3=LDQ((const __m256i *)(b+j+12));
__m256i z0=LDQ((const __m256i *)(e+j)), z1=LDQ((const __m256i *)(e+j+4));
__m256i z2=LDQ((const __m256i *)(e+j+8)), z3=LDQ((const __m256i *)(e+j+12));
STQNT((__m256i *)(dpa+j), _mm256_sub_epi64(w0,y0));
STQNT((__m256i *)(dpa+j+4), _mm256_sub_epi64(w1,y1));
STQNT((__m256i *)(dpa+j+8), _mm256_sub_epi64(w2,y2));
STQNT((__m256i *)(dpa+j+12), _mm256_sub_epi64(w3,y3));
STQNT((__m256i *)(dpb+j), _mm256_sub_epi64(x0,z0));
STQNT((__m256i *)(dpb+j+4), _mm256_sub_epi64(x1,z1));
STQNT((__m256i *)(dpb+j+8), _mm256_sub_epi64(x2,z2));
STQNT((__m256i *)(dpb+j+12), _mm256_sub_epi64(x3,z3));
}
for (; j < n; j += 4) {
STQNT((__m256i *)(dpa+j), _mm256_sub_epi64(LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
STQNT((__m256i *)(dpb+j), _mm256_sub_epi64(LDQ((const __m256i *)(c+j)), LDQ((const __m256i *)(e+j))));
}
}
}
for (; i < n; i++) {
long long *dp = dst + (size_t)i * n;
const long long *a = p0 + (size_t)i * ld, *b = p1 + (size_t)i * ld;
_mm_prefetch((const char *)(p0 + (size_t)(i + 2) * ld), _MM_HINT_T0);
_mm_prefetch((const char *)(p1 + (size_t)(i + 2) * ld), _MM_HINT_T0);
int j = 0;
if (!s1neg) {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)), w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
w0=_mm256_add_epi64(w0,LDQ((const __m256i *)(b+j)));
w1=_mm256_add_epi64(w1,LDQ((const __m256i *)(b+j+4)));
w2=_mm256_add_epi64(w2,LDQ((const __m256i *)(b+j+8)));
w3=_mm256_add_epi64(w3,LDQ((const __m256i *)(b+j+12)));
STQNT((__m256i *)(dp+j), w0); STQNT((__m256i *)(dp+j+4), w1);
STQNT((__m256i *)(dp+j+8), w2); STQNT((__m256i *)(dp+j+12), w3);
}
for (; j < n; j += 4) STQNT((__m256i *)(dp+j), _mm256_add_epi64(
LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
} else {
for (; j + 16 <= n; j += 16) {
__m256i w0=LDQ((const __m256i *)(a+j)), w1=LDQ((const __m256i *)(a+j+4));
__m256i w2=LDQ((const __m256i *)(a+j+8)), w3=LDQ((const __m256i *)(a+j+12));
w0=_mm256_sub_epi64(w0,LDQ((const __m256i *)(b+j)));
w1=_mm256_sub_epi64(w1,LDQ((const __m256i *)(b+j+4)));
w2=_mm256_sub_epi64(w2,LDQ((const __m256i *)(b+j+8)));
w3=_mm256_sub_epi64(w3,LDQ((const __m256i *)(b+j+12)));
STQNT((__m256i *)(dp+j), w0); STQNT((__m256i *)(dp+j+4), w1);
STQNT((__m256i *)(dp+j+8), w2); STQNT((__m256i *)(dp+j+12), w3);
}
for (; j < n; j += 4) STQNT((__m256i *)(dp+j), _mm256_sub_epi64(
LDQ((const __m256i *)(a+j)), LDQ((const __m256i *)(b+j))));
}
}
}
/* [ml1d_] NT cold-write knife. matop's destination is a *fresh* n*n plane written
wholesale and consumed only after >L3 of intervening traffic (see head). For n>=512
that plane is 2 MB and provably evicted before its consumer -> the write-allocate
(RFO) read is pure waste -> vmovntdq. For n==256 the plane is 512 KB and stays
L3-resident -> NT would only destroy the consumer's hits (MEASURED +2.10 ms) -> keep
the cacheable store. cmb3's destination is never read back at all -> NT for the
large planes; the 2695 sub-64 leaf tiles stay cacheable. */
static void matop2(const Opd *o, int n, long long *dst) {
if (n >= 512) { matop2_ntx(o, n, dst); return; }
matop2_p(o, n, dst);
}
static void cmb3(long long *d, int ldd, const long long *a, const long long *b,
const long long *c, int ld, int rows, int ncol) {
if (rows >= 256 && ((unsigned long)d & 31u) == 0 && (ldd & 3) == 0) {
cmb3_ntx(d, ldd, a, b, c, ld, rows, ncol); return;
}
cmb3_p(d, ldd, a, b, c, ld, rows, ncol);
}
static void matop(const Opd *o, int n, long long *dst) {
const int nt = o->nt; const long long *P = o->p; const int ld = o->ld;
if (nt == 2) { matop2(o, n, dst); return; }
for (int i = 0; i < n; i++) {
long long *dp = dst + (size_t)i * n;
int j = 0;
for (; j + 16 <= n; j += 16) {
__m256i v0=_mm256_setzero_si256(),v1=v0,v2=v0,v3=v0;
for (int t = 0; t < nt; t++) {
const long long *sp = P + (size_t)(o->r[t] + i) * ld + o->c[t] + j;
__m256i w0=LDQ((const __m256i *)sp), w1=LDQ((const __m256i *)(sp+4));
__m256i w2=LDQ((const __m256i *)(sp+8)), w3=LDQ((const __m256i *)(sp+12));
if (o->s[t] < 0) { v0=_mm256_sub_epi64(v0,w0); v1=_mm256_sub_epi64(v1,w1);
v2=_mm256_sub_epi64(v2,w2); v3=_mm256_sub_epi64(v3,w3); }
else { v0=_mm256_add_epi64(v0,w0); v1=_mm256_add_epi64(v1,w1);
v2=_mm256_add_epi64(v2,w2); v3=_mm256_add_epi64(v3,w3); }
}
STQ((__m256i *)(dp+j), v0); STQ((__m256i *)(dp+j+4), v1);
STQ((__m256i *)(dp+j+8), v2); STQ((__m256i *)(dp+j+12), v3);
}
for (; j < n; j += 4) {
__m256i v = _mm256_setzero_si256();
for (int t = 0; t < nt; t++) {
__m256i w = LDQ((const __m256i *)(P + (size_t)(o->r[t] + i) * ld + o->c[t] + j));
v = (o->s[t] < 0) ? _mm256_sub_epi64(v, w) : _mm256_add_epi64(v, w);
}
STQ((__m256i *)(dp + j), v);
}
}
}
/* the strassen node whose 7 children are leaves, driven from one fused pass */
static void leaflevel(int n, const Opd *Ao, const Opd *Bo, long long *d, int ld, int op,
long long *T) {
const int h = n >> 1;
long long *Q[4];
Q[0] = d; Q[1] = d + h; Q[2] = d + (size_t)h * ld; Q[3] = Q[2] + h;
int wr[4] = {0, 0, 0, 0};
buildA(Ao); buildB(Bo);
#define MKL(qi, si, qj, sj, L) 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 (oa == 0) { \
leafP(h, (L), Q[qi], ld, 0, 0, 0, 0); \
cmb(Q[qj], ld, Q[qi], ld, h, (sj), ob); \
} else { \
leafP(h, (L), T, h, 0, 0, 0, 0); \
cmb2(Q[qi], ld, oa, (si), Q[qj], ld, ob, (sj), T, h, h); \
} \
} while (0)
#define MKL1(qi, si, L) do { \
int oa = wr[qi] ? effop(1, (si)) : effop(op, (si)); \
wr[qi] = 1; \
leafP(h, (L), Q[qi], ld, oa, 0, 0, 0); \
} while (0)
{ MKL(0, 1, 3, 1, 0); } /* M1 = (A11+A22)(B11+B22) -> +C11, +C22 */
{ MKL(2, 1, 3, 2, 1); } /* M2 = (A21+A22)B11 -> +C21, -C22 */
{ MKL(1, 1, 3, 1, 2); } /* M3 = A11(B12-B22) -> +C12, +C22 */
{ MKL(0, 1, 2, 1, 3); } /* M4 = A22(B21-B11) -> +C11, +C21 */
{ MKL(0, 2, 1, 1, 4); } /* M5 = (A11+A12)B22 -> -C11, +C12 */
{ MKL1(3, 1, 5); } /* M6 = (A21-A11)(B11+B12) -> +C22 */
{ MKL1(0, 1, 6); } /* M7 = (A12-A22)(B21+B22) -> +C11 */
#undef MKL
#undef MKL1
}
static void buildA_blk(const Opd *o, int i0, int i1) {
const int nt = o->nt; const long long *P = o->p; const int ld = o->ld;
if (nt == 2) {
const size_t qld = (size_t)QN * ld;
const long long *B0 = P + (size_t)o->r[0] * ld + o->c[0];
const long long *B1 = P + (size_t)o->r[1] * ld + o->c[1];
const int n0 = (o->s[0] < 0), n1 = (o->s[1] < 0);
const long long *pA = n0 ? B1 : B0, *pB = n0 ? B0 : B1;
const int sub = (n0 != n1), negall = (n0 && n1);
for (int i = i0; i < i1; i++) {
const size_t pi = (size_t)(i - i0);
long long *p0 = PAN_A[0] + pi * KCP, *p1 = PAN_A[1] + pi * KCP,
*p2 = PAN_A[2] + pi * KCP, *p3 = PAN_A[3] + pi * KCP,
*p4 = PAN_A[4] + pi * KCP, *p5 = PAN_A[5] + pi * KCP,
*p6 = PAN_A[6] + pi * KCP;
const long long *f0 = pA + (size_t)i * ld, *f1 = pB + (size_t)i * ld;
__m256i a[2][4], w, z = _mm256_setzero_si256();
for (int k = 0; k < QN; k += 8) {
/* BOTH 4-blocks' loads issue before either block's stores: the shipped loop
cannot do this (the "memory" clobber forbids gcc the reorder). */
for (int b = 0; b < 2; b++) {
const int kk = k + b * 4;
a[b][0] = LDQ((const __m256i *)(f0 + kk));
w = LDQ((const __m256i *)(f1 + kk));
a[b][0] = sub ? _mm256_sub_epi64(a[b][0], w) : _mm256_add_epi64(a[b][0], w);
a[b][1] = LDQ((const __m256i *)(f0 + QN + kk));
w = LDQ((const __m256i *)(f1 + QN + kk));
a[b][1] = sub ? _mm256_sub_epi64(a[b][1], w) : _mm256_add_epi64(a[b][1], w);
a[b][2] = LDQ((const __m256i *)(f0 + qld + kk));
w = LDQ((const __m256i *)(f1 + qld + kk));
a[b][2] = sub ? _mm256_sub_epi64(a[b][2], w) : _mm256_add_epi64(a[b][2], w);
a[b][3] = LDQ((const __m256i *)(f0 + qld + QN + kk));
w = LDQ((const __m256i *)(f1 + qld + QN + kk));
a[b][3] = sub ? _mm256_sub_epi64(a[b][3], w) : _mm256_add_epi64(a[b][3], w);
if (negall) { for (int q = 0; q < 4; q++) a[b][q] = _mm256_sub_epi64(z, a[b][q]); }
}
for (int b = 0; b < 2; b++) {
const int kk = k + b * 4;
STQ((__m256i *)(p0 + kk), _mm256_add_epi64(a[b][0], a[b][3])); /* M1 A11+A22 */
STQ((__m256i *)(p1 + kk), _mm256_add_epi64(a[b][2], a[b][3])); /* M2 A21+A22 */
STQ((__m256i *)(p2 + kk), a[b][0]); /* M3 A11 */
STQ((__m256i *)(p3 + kk), a[b][3]); /* M4 A22 */
STQ((__m256i *)(p4 + kk), _mm256_add_epi64(a[b][0], a[b][1])); /* M5 A11+A12 */
STQ((__m256i *)(p5 + kk), _mm256_sub_epi64(a[b][2], a[b][0])); /* M6 A21-A11 */
STQ((__m256i *)(p6 + kk), _mm256_sub_epi64(a[b][1], a[b][3])); /* M7 A12-A22 */
}
}
}
return;
}
{ /* generic hoisted path: identical arithmetic, pointers and signs hoisted */
const long long *Bp[8]; int sg[8];
for (int t = 0; t < nt; t++) { Bp[t] = P + (size_t)o->r[t] * ld + o->c[t]; sg[t] = (o->s[t] < 0); }
const size_t qld = (size_t)QN * ld;
for (int i = i0; i < i1; i++) {
const size_t pi = (size_t)(i - i0);
long long *p0 = PAN_A[0] + pi * KCP, *p1 = PAN_A[1] + pi * KCP,
*p2 = PAN_A[2] + pi * KCP, *p3 = PAN_A[3] + pi * KCP,
*p4 = PAN_A[4] + pi * KCP, *p5 = PAN_A[5] + pi * KCP,
*p6 = PAN_A[6] + pi * KCP;
for (int k = 0; k < QN; k += 4) {
__m256i a11, a12, a21, a22;
if (nt == 1 && !sg[0]) {
a11 = LDQ((const __m256i *)(Bp[0] + (size_t)i * ld + k));
a12 = LDQ((const __m256i *)(Bp[0] + (size_t)i * ld + QN + k));
a21 = LDQ((const __m256i *)(Bp[0] + (size_t)(i + QN) * ld + k));
a22 = LDQ((const __m256i *)(Bp[0] + (size_t)(i + QN) * ld + QN + k));
} else {
a11 = _mm256_setzero_si256(); a12 = a11; a21 = a11; a22 = a11;
for (int t = 0; t < nt; t++) {
const long long *b = Bp[t] + (size_t)i * ld;
const int neg = sg[t]; __m256i w;
w = LDQ((const __m256i *)(b + k));
a11 = neg ? _mm256_sub_epi64(a11, w) : _mm256_add_epi64(a11, w);
w = LDQ((const __m256i *)(b + QN + k));
a12 = neg ? _mm256_sub_epi64(a12, w) : _mm256_add_epi64(a12, w);
w = LDQ((const __m256i *)(b + (size_t)QN * ld + k));
a21 = neg ? _mm256_sub_epi64(a21, w) : _mm256_add_epi64(a21, w);
w = LDQ((const __m256i *)(b + (size_t)QN * ld + QN + k));
a22 = neg ? _mm256_sub_epi64(a22, w) : _mm256_add_epi64(a22, w);
}
}
STQ((__m256i *)(p0 + k), _mm256_add_epi64(a11, a22)); /* M1 A11+A22 */
STQ((__m256i *)(p1 + k), _mm256_add_epi64(a21, a22)); /* M2 A21+A22 */
STQ((__m256i *)(p2 + k), a11); /* M3 A11 */
STQ((__m256i *)(p3 + k), a22); /* M4 A22 */
STQ((__m256i *)(p4 + k), _mm256_add_epi64(a11, a12)); /* M5 A11+A12 */
STQ((__m256i *)(p5 + k), _mm256_sub_epi64(a21, a11)); /* M6 A21-A11 */
STQ((__m256i *)(p6 + k), _mm256_sub_epi64(a12, a22)); /* M7 A12-A22 */
}
}
}
}
static void cmb_blk(long long *dst, int ldd, const long long *T, int ldt, int h, int sign, int op, int i0, int i1) {
__m256i sgn = _mm256_set1_epi64x(sign < 0 ? -1 : 0);
for (int i = i0; i < i1; i++) {
const long long *t = T + (size_t)i * ldt;
long long *c = dst + (size_t)i * ldd;
int j = 0;
/* ARM cmbu16: the j-loop's inner step raised 4 -> 16 elements so FOUR loads
are issued before the first store. `LDQ` carries a "memory" clobber, so
the shipped 4-element step cannot be software-pipelined by gcc; here the
pipelining is explicit. Same loads, same stores, same addresses. */
if (op == 0) {
for (; j + 16 <= h; j += 16) {
__m256i v0 = LDQ((const __m256i *)(t + j)), v1 = LDQ((const __m256i *)(t + j + 4));
__m256i v2 = LDQ((const __m256i *)(t + j + 8)), v3 = LDQ((const __m256i *)(t + j + 12));
if (sign < 0) { __m256i z = _mm256_setzero_si256();
v0 = _mm256_sub_epi64(z, v0); v1 = _mm256_sub_epi64(z, v1);
v2 = _mm256_sub_epi64(z, v2); v3 = _mm256_sub_epi64(z, v3); }
STQ((__m256i *)(c + j), v0); STQ((__m256i *)(c + j + 4), v1);
STQ((__m256i *)(c + j + 8), v2); STQ((__m256i *)(c + j + 12), v3);
}
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 + 16 <= h; j += 16) {
__m256i v0 = LDQ((const __m256i *)(t + j)), v1 = LDQ((const __m256i *)(t + j + 4));
__m256i v2 = LDQ((const __m256i *)(t + j + 8)), v3 = LDQ((const __m256i *)(t + j + 12));
__m256i w0 = LDQ((const __m256i *)(c + j)), w1 = LDQ((const __m256i *)(c + j + 4));
__m256i w2 = LDQ((const __m256i *)(c + j + 8)), w3 = LDQ((const __m256i *)(c + j + 12));
if ((op == 1) != (sign < 0)) {
STQ((__m256i *)(c + j), _mm256_add_epi64(w0, v0));
STQ((__m256i *)(c + j + 4), _mm256_add_epi64(w1, v1));
STQ((__m256i *)(c + j + 8), _mm256_add_epi64(w2, v2));
STQ((__m256i *)(c + j + 12), _mm256_add_epi64(w3, v3));
} else {
STQ((__m256i *)(c + j), _mm256_sub_epi64(w0, v0));
STQ((__m256i *)(c + j + 4), _mm256_sub_epi64(w1, v1));
STQ((__m256i *)(c + j + 8), _mm256_sub_epi64(w2, v2));
STQ((__m256i *)(c + j + 12), _mm256_sub_epi64(w3, v3));
}
}
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]; }
}
}
}
static void cmb2_blk(long long *d1,int ld1,int op1,int s1,
long long *d2,int ld2,int op2,int s2,
const long long *T,int ldt,int h, int i0, int i1) {
for (int i = i0; i < i1; i++) {
const long long *t = T + (size_t)i * ldt;
long long *c1 = d1 + (size_t)i * ld1, *c2 = d2 + (size_t)i * ld2;
int j = 0;
for (; j + 16 <= h; j += 16) {
__m256i v0=LDQ((const __m256i*)(t+j)), v1=LDQ((const __m256i*)(t+j+4));
__m256i v2=LDQ((const __m256i*)(t+j+8)), v3=LDQ((const __m256i*)(t+j+12));
__m256i a0=LDQ((const __m256i*)(c1+j)), a1=LDQ((const __m256i*)(c1+j+4));
__m256i a2=LDQ((const __m256i*)(c1+j+8)), a3=LDQ((const __m256i*)(c1+j+12));
__m256i b0=LDQ((const __m256i*)(c2+j)), b1=LDQ((const __m256i*)(c2+j+4));
__m256i b2=LDQ((const __m256i*)(c2+j+8)), b3=LDQ((const __m256i*)(c2+j+12));
if ((op1 == 1) != (s1 < 0)) { STQ((__m256i*)(c1+j),_mm256_add_epi64(a0,v0));
STQ((__m256i*)(c1+j+4),_mm256_add_epi64(a1,v1)); STQ((__m256i*)(c1+j+8),_mm256_add_epi64(a2,v2));
STQ((__m256i*)(c1+j+12),_mm256_add_epi64(a3,v3)); }
else { STQ((__m256i*)(c1+j),_mm256_sub_epi64(a0,v0));
STQ((__m256i*)(c1+j+4),_mm256_sub_epi64(a1,v1)); STQ((__m256i*)(c1+j+8),_mm256_sub_epi64(a2,v2));
STQ((__m256i*)(c1+j+12),_mm256_sub_epi64(a3,v3)); }
if ((op2 == 1) != (s2 < 0)) { STQ((__m256i*)(c2+j),_mm256_add_epi64(b0,v0));
STQ((__m256i*)(c2+j+4),_mm256_add_epi64(b1,v1)); STQ((__m256i*)(c2+j+8),_mm256_add_epi64(b2,v2));
STQ((__m256i*)(c2+j+12),_mm256_add_epi64(b3,v3)); }
else { STQ((__m256i*)(c2+j),_mm256_sub_epi64(b0,v0));
STQ((__m256i*)(c2+j+4),_mm256_sub_epi64(b1,v1)); STQ((__m256i*)(c2+j+8),_mm256_sub_epi64(b2,v2));
STQ((__m256i*)(c2+j+12),_mm256_sub_epi64(b3,v3)); }
}
for (; j + 4 <= h; j += 4) {
__m256i v = LDQ((const __m256i *)(t + j));
__m256i a = LDQ((const __m256i *)(c1 + j));
STQ((__m256i *)(c1 + j), ((op1 == 1) != (s1 < 0)) ? _mm256_add_epi64(a,v)
: _mm256_sub_epi64(a,v));
__m256i b = LDQ((const __m256i *)(c2 + j));
STQ((__m256i *)(c2 + j), ((op2 == 1) != (s2 < 0)) ? _mm256_add_epi64(b,v)
: _mm256_sub_epi64(b,v));
}
for (; j < h; j++) {
long long v = t[j];
if ((op1 == 1) != (s1 < 0)) c1[j] += v; else c1[j] -= v;
if ((op2 == 1) != (s2 < 0)) c2[j] += v; else c2[j] -= v;
}
}
}
/* row-blocked leaflevel: the A-side build for rows [ic0,ic1) of all 7 panels is
done immediately before those rows are consumed. Every operand row is still
read exactly once, the instruction stream of each phase is unchanged.
The per-product destination operators (oa/ob) are functions of the ORIGINAL
sequential wr[] evolution, so they are computed ONCE here, before the loop. */
static void leaflevel3(int n, const Opd *Ao, const Opd *Bo, long long *d, int ld, int op,
long long *T) {
const int h = n >> 1;
long long *Q[4];
Q[0] = d; Q[1] = d + h; Q[2] = d + (size_t)h * ld; Q[3] = Q[2] + h;
int OA[7], OB[7];
{
int wr[4] = {0, 0, 0, 0};
const int qi[7] = {0, 2, 1, 0, 0, 3, 0};
const int si[7] = {1, 1, 1, 1, 2, 1, 1};
const int qj[7] = {3, 3, 3, 2, 1, -1, -1};
const int sj[7] = {1, 2, 1, 1, 1, 0, 0};
for (int L = 0; L < 7; L++) {
OA[L] = wr[qi[L]] ? effop(1, si[L]) : effop(op, si[L]);
wr[qi[L]] = 1;
if (qj[L] >= 0) { OB[L] = wr[qj[L]] ? effop(1, sj[L]) : effop(op, sj[L]); wr[qj[L]] = 1; }
}
}
buildB(Bo);
/* ROW-BLOCKED A BUILD (l1k7_rv `rba`): the A panel depends on the ROW index and
NOT on the column group, and buildA_blk reads every operand row exactly ONCE
across the blocks, so blocking the build by MR rows costs no extra read -- while
the live A window drops from 7 x 64 x 608 B = 268 KB to 7 x MR x 512 B = 21.5 KB.
B stays a single up-front pass (its k is the contraction index and every mikro
call reads all KC of it). Only the ORDER of the mikro calls changes. */
{
static const int DQ1[7] = {0, 2, 1, 0, 0, 3, 0};
static const int DQ2[7] = {3, 3, 3, 2, 1, -1, -1};
static const int SIx[7] = {1, 1, 1, 1, 2, 1, 1};
static const int SJx[7] = {1, 2, 1, 1, 1, 0, 0};
for (int ic0 = 0; ic0 < h; ic0 += MR) {
int ic1 = ic0 + MR; if (ic1 > h) ic1 = h;
buildA_blk(Ao, ic0, ic1);
if (op == 0) {
/* FUS3 at the leaf: M1->Q[0], M2->Q[2], M3->Q[1], then ONE C22 pass. */
leafP_blk(h, 0, ic0, ic1, Q[0], ld, 0, 0, 0, 0);
leafP_blk(h, 1, ic0, ic1, Q[2], ld, 0, 0, 0, 0);
leafP_blk(h, 2, ic0, ic1, Q[1], ld, 0, 0, 0, 0);
cmb3_blk(Q[3], ld, Q[0], Q[2], Q[1], ld, h, ic0, ic1);
for (int L = 3; L < 7; L++) {
long long *D1 = Q[DQ1[L]];
long long *D2 = (DQ2[L] >= 0) ? Q[DQ2[L]] : 0;
const int oa = OA[L], ob = OB[L];
if (D2 == 0) {
leafP_blk(h, L, ic0, ic1, D1, ld, oa, 0, 0, 0);
} else if (oa == 0) {
leafP_blk(h, L, ic0, ic1, D1, ld, 0, 0, 0, 0);
cmb_blk(D2, ld, D1, ld, h, SJx[L], ob, ic0, ic1);
} else {
leafP_blk(h, L, ic0, ic1, T, h, 0, 0, 0, 0);
cmb2_blk(D1, ld, oa, SIx[L], D2, ld, ob, SJx[L], T, h, h, ic0, ic1);
}
}
} else {
for (int L = 0; L < 7; L++) {
long long *D1 = Q[DQ1[L]];
long long *D2 = (DQ2[L] >= 0) ? Q[DQ2[L]] : 0;
const int oa = OA[L], ob = OB[L];
if (D2 == 0) {
leafP_blk(h, L, ic0, ic1, D1, ld, oa, 0, 0, 0);
} else if (oa == 0) {
leafP_blk(h, L, ic0, ic1, D1, ld, 0, 0, 0, 0);
cmb_blk(D2, ld, D1, ld, h, SJx[L], ob, ic0, ic1);
} else {
leafP_blk(h, L, ic0, ic1, T, h, 0, 0, 0, 0);
cmb2_blk(D1, ld, oa, SIx[L], D2, ld, ob, SJx[L], T, h, h, ic0, ic1);
}
}
}
}
}
}
static void mmxR(int n, const Opd *Ao0, const Opd *Bo0, long long *d, int ld, int op,
long long *T) {
const Opd *Ao = Ao0, *Bo = Bo0; Opd A1, B1;
if (n <= CUTL || (n & 1) || Ao->nt > 4 || Bo->nt > 4) {
leaf(n, Ao, Bo, d, ld, op, 0, 0, 0); return; }
if (n == 2 * CUTL) {
{ leaflevel3(n, Ao, Bo, d, ld, op, T); return; }
leaflevel(n, Ao, Bo, d, ld, op, T); return; }
if (n <= 8 * CUTL && (Ao->nt > 1 || Bo->nt > 1)) {
int Lm = (n > 4 * CUTL) ? 0 : (n > 2 * CUTL) ? 1 : 2;
if (Ao->nt > 1) { matop(Ao, n, matbufA[Lm]); A1 = o1(matbufA[Lm], n); Ao = &A1; }
if (Bo->nt > 1) { matop(Bo, n, matbufB[Lm]); B1 = o1(matbufB[Lm], n); Bo = &B1; }
}
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 (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); \
cmb2(Q[qi], ld, oa, (si), Q[qj], ld, ob, (sj), T, h, h); \
} \
} 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)
if (op == 0) {
/* FUS3: M1, M2, M3 all write C22. Their three sources (C11, C21, C12) are
live simultaneously, so one pass replaces three. Same mmxR call order. */
wr[0] = 1; wr[3] = 1;
{ Opd Ma = oadd(&A11, &A22, 1), Mb = oadd(&B11, &B22, 1);
mmxR(h, &Ma, &Mb, Q[0], ld, 0, T); }
wr[2] = 1;
{ Opd Ma = oadd(&A21, &A22, 1), Mb = B11;
mmxR(h, &Ma, &Mb, Q[2], ld, 0, T); }
wr[1] = 1;
{ Opd Ma = A11, Mb = oadd(&B12, &B22, -1);
mmxR(h, &Ma, &Mb, Q[1], ld, 0, T); }
cmb3(Q[3], ld, Q[0], Q[2], Q[1], ld, h, h);
} else {
{ /* 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
}
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); return; }
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 | 108.39 ms | 15 MB + 956 KB | Accepted | Score: 100 | 显示更多 |