// ===== d4bc_ 第二发（跨叶携带）：新增改动 =====
// 本发在本文件上方 [R1] 基座（#123995，= 本文件上一条 REFERENCES 所述）之上，**再加一处定点改动**：
// **把每个 512 叶"第一个 ic 块"的 A 面板打包（下称 P0）从"叶内同步拷 123 KB"改成"由上一叶最后一个
// ic 块的核调用携带出来"**（跨叶携带；两缓冲奇偶交替：叶 r 起始缓冲 = r&1，其末块把下一叶首块面板
// 写进 cur^1）。
// 依据链（全部判题机 n=2048 免费通道、同基座同驱动、取 r1；量具口径沿用 `d4ba_` 封条 (2191) 的消融 Δ）：
// 1. ★**P0 定价（本席同批定价臂）**：`p0skipall`（整体删 P0）= **−2.19e6 拍 = −0.236%**（同批 base
// 928.974e6 ↔ 926.780e6）⇒ P0 是一块"每叶工作量与尺度无关"的真池（×343 叶 ≈ **−4.3 ms@4096**）。
// 2. ★**P0 不可删、只能转移**：它是该叶首块 A 面板的唯一来源；而"携带进核"这件事**是净收益** ——
// 反面对照臂 `anoap`（把 sign==0 的叶的携带整段跳过）= **+0.36% 更慢**（同批）⇒ 把 P0 的工作搬进
// 携带核的富余端口，正是可行路线。
// 3. 本发实测（判题机 n=2048，同族同驱动同批）：**−1.0046e6 拍 = −0.109%**，C 全量 FNV
// `c7f1e86920290383` 与基座逐位相同 ⇒ **值安全**（不是定价臂）。
// 4. 逐位等价闸门：**n=4096 + 两组独立随机输入种子（12345 / 98765）**，逐 256 行分 16 块全量 FNV
// ⇒ **32/32 块与在役件 #123995 同源编译的对照件逐位相同** ✓。
// 5. ★**开发记录（诚实留档）**：本改动第一次实现时**在 `rec` 的叶分支上引入了一个 7× 冗余环**
// （把 7 个积的循环套进了原有的 `for(r=0..6)` 外层循环 ⇒ 每节点跑 49 次而非 7 次），
// 本机 21.3 s vs 3.8 s（5.6×）而**值仍逐位正确**（重复计算同一批输出）—— 这解释了初版的
// "值等价但 5.6× 慢"之谜，与携带机制本身无关。修正环结构后本机 3.86 s（与基座同档）。
// 6. 口径：本发与 [R1] 同属"每叶工作量与尺度无关"类（每叶首块面板固定 123 KB、343 叶全覆盖）
// ⇒ 按 `d4bb_` 封条 (2195) 的机理，预期兑现 ≈ Δ×叶比(7) ≈ **−1.96 ms**，用于把缺口（1.667725 ms）转绿。
// ===== [R1] 以下为上一发的原文（未改） =====
// ===== REFERENCES =====
// [R1] 本账号 saffah_cc_v41_agg1 提交 **#123993**（`d4bb_` 席，1859.660926 ms，本账号在册最好）
// <https://duck.ac/submission/123993>
// 用途：**本提交的直接基座**。代码区与该件除本席（`d4bc_`）标注的那处定点删除外逐字节相同；
// 该件保留的上游 Credit 链原样保留，未作删改。
// [R2] 同账号同题在册账与封条：`d4bb_` 封条 **(2195)**（消融真值台账 + 12.1× 尺度高估律）·
// `d4ba_` 封条 **(2191)**（「贴背内联核的 rdtscp 括弧不可加、只有消融 Δ 可加」律）·
// `d4az_` 封条 **(2184)** · `d4ax_` 封条 **(2178)**。本席入场读数：`mine 1859.660926`（#123993）·
// `T 1876.693695`（#122097，`saffah_codex_6a_agg3`）· 严支 `1857.927758` ⇒ 缺 **1.733168 ms (0.093%)**。
// [R3] 只读借用的队内工件（同账号、同题）：`problems/mmmd4k/work/d4av__tmp/d4av_base.cpp`
// （= #123975 的逐字节副本，md5 fcd27a62ab44c5866b25daba93ec9466，作消融/差分基座）·
// `work/d4ba__tmp/d4ba_mkab.py`（「一旗一消融」生成器形态，本席照此自建 `d4bc_mkarm.py`）·
// `work/d4bb__tmp/d4bb_mk4096.py`（n=4096 逐 256 行 16 块 FNV 闸门形态）。
//
// ===== 思路 =====
// **本发只做一件事：删掉 `kernel_tail16` 内层循环里的 8 条 `prefetcht0`（不动循环结构、不动任何
// 数据通路、不动一条 FMA）。** 全部读数口径沿用 `d4ba_` 封条 (2191)：判题机 n=2048 免费通道消融 Δ、
// 同基座同驱动、取轮次 r1。
// 1. 尺子先对齐：base 928.786e6（前席 928.793e6）· packB 消融 −29.99e6（册 −30.66e6）·
// `kernel<SG>` 消融 −165.55e6（册 −164.86e6）· PEND 族消融（入口 `active=false` 版）−30.71e6
// （册 −30.31e6）⇒ 四数全部对上。
// 2. 三个点名池逐一复核后判死（细节入 notes）：packB 读侧 20.17e6 / 写侧 23.41e6（两臂非可加，
// 合价 30.66e6），本机汇编证明该环已被 GCC 展开到 64 列/迭代、环开销≈0，且融合表在现有 `lperm`
// 下已达「两个纯象限积 ⇒ 最多 2 次融合」的理论上限 ⇒ 余下全是"必须搬的字节"；
// PEND 30.31e6 里 14.76e6 是 n=2048 测试台独有（4096 走 `root_last`/`pend2_partial` 重叠路径）；
// SG 携带核的 18.5e6 富余就是 A 面板打包本身的必要流量（预取已被前席判负）。
// 3. ★本席新价（此前从未定价）：`kernel_tail16`（每个 512 叶最后一个 ic 块的 2 行 × 16 列输出尾）
// = **−4.26e6 拍 = 0.459%**；2732 拍/调用 vs FMA 地板 2048 拍 = **1.33×**。
// 同批还判出 **`kernel2` 是死代码**（掏空后 C 的 FNV 与基座逐位相同、计数恒 0 ⇒ 该函数在
// n=512 叶上永不触发）。
// 4. ★机理：该环每 k 步 = 2 条广播 + 4 条 b 载入 + **2 条预取** = 8 个载入口 uop（2 个载入口
// ⇒ 4 拍），与 8 条 FMA（2 个 FMA 口 ⇒ 4 拍）**正好打平** ⇒ 任一点抖动都直接堆到关键路径。
// 删掉预取后载入腿降到 3 拍，FMA 成为唯一瓶颈。
// 5. ★反面对照（同批实测、且都是值安全臂）：删 `kernel<SG>` 的 8 条预取 = **+2.44e6（更慢）**、
// 删 `kernel_plain` 的 8 条预取 = **+9.43e6（大幅更慢）** ⇒ 那两条是承重预取（掩盖 DRAM 延迟），
// 只有 tail16 因载入口正好饱和才是"预取反而害事"的特例。这一对照本身就是机理的证据。
// 6. 本发实测 Δ（判题机 n=2048）：**−0.2546e6 拍（−0.027%）**，C 全量 FNV `c7f1e86920290383`
// 与基座逐位相同 ⇒ 口径上**值安全**，不是定价臂。
// 7. 逐位等价闸门（提交硬闸）：**n=4096 + 两组独立随机输入种子（12345 / 98765）**，逐 256 行分
// 16 块全量 FNV ⇒ **32/32 块与在役件 #123993 同源编译的对照件逐位相同** ✓（本机只用于正确性，
// 量级判读一律走判题机）。
// 8. ★口径与预期（承接 `d4bb_` 封条 (2195) 的 12.1× 尺度高估律）：本刀是**计算/端口类**（改的是
// 每 k 步的端口配平，与尺度无关；每叶尾工作量固定、343 叶全覆盖、n=2048 与 n=4096 命中同一
// 代码路径）⇒ 按封条 (2103) 的分档律属"可信尺子"档，预期兑现 ≈ Δ × 叶比(7) ≈ **−0.49 ms**，
// 用于把缺口从 1.733 ms 收窄；同时本发也在册面上**实测"计算类免费通道读数在 mmmd4k 4096 上
// 的兑现比"**（本题此前只有访存/流类的一个校准点，12.1×，无计算类校准点）。
// ===== REFERENCES =====
// [R1] 本账号 saffah_cc_v41_agg1 提交 **#123975**（`d4ak_`/`d4av_` 席，1860.379542 ms，本账号在册最好）
// <https://duck.ac/submission/123975>
// 用途：**本提交的直接基座**。代码区与该件除下面 `d4bb_` 标注的那处定点改动外逐字节相同；
// 该件自身保留的 Credits 链原样保留在紧随其后的注释区里，未作删改（本席只删了正文代码区的空行）。
// [R2] 同账号同题在册账与封条：`d4ba_` 封条 **(2191)**（消融真值台账 + 「贴背内联核的 rdtscp 括弧不可加、
// 只有消融 Δ 可加」律）· `d4az_` 封条 **(2184)**（皮带距离剂量曲线 + mikro 地板律）·
// `d4ax_` 封条 **(2178)**。本席（`d4bb_`）入场读数：`mine 1860.379542`（#123975）·
// `T 1876.693695`（#122097，`saffah_codex_6a_agg3`）· 严支 `1857.927758` ⇒ 缺 **2.451784 ms (0.132%)**。
// [R3] 只读借用的队内工件（同账号、同题）：`problems/mmmd4k/work/d4av__tmp/d4av_base.cpp`
// （= #123975 的逐字节副本，md5 fcd27a62ab44c5866b25daba93ec9466，作消融/差分基座）·
// `work/d4ba__tmp/d4ba_mkab.py`（「一旗一消融」生成器形态，本席照此自建 `d4bb_mkarm.py`）。
//
// ===== 思路 =====
// **本发只做一件事：把最外层传送带节点的两个 s3 暂存面板，从「写进 C 的两个象限、按 ldc 跨步、
// 每行只写 h 列」改成「写进 C 的上半区、按 h 连续存放」（s3a -> c，s3b -> c+h*h），并让它唯一
// 的 s3 消费者（第 7 个子节点 `rec<h>(h,{c,ldc},{c+h,ldc},out[6],...)`）改读 `{c,h}` / `{c+h*h,h}`。**
//
// 依据链（全部判题机 n=2048 免费通道，同基座同驱动，取 r1；量具口径沿用 `d4ba_` 封条 (2191)）：
// 1. 先复现真值台账：plain −619.65e6 · 皮带 −63.22e6 · SG −165.78e6 拍（前席 −619.65 / −63.11 /
// −164.86e6 逐项对上）⇒ 尺子对齐。
// 2. 逐项差分定位皮带成本：皮带**源指针**全指向一块热面板 ⇒ **−10.78e6 拍**；皮带**汇指针**同样
// 处理 ⇒ **+0.50e6（更慢）**；传送带槽改调 `kernel_plain`（面板不再生成）⇒ −10.99e6。
// ⇒ **成本全在"读侧"**，写侧本来就是自由的；而 s3 面板恰是那条"每行只写 h 列、跨 2h 列"的
// 稀疏流（读侧与写侧都被它放大一倍footprint）。
// 3. 本改动实测（判题机 n=2048，两次独立发件）：**−4.34e6 / −4.17e6 拍 = −0.47%**，且 **C 的
// 33.5 MB 全量 FNV 与基座逐位相同**（`c7f1e86920290383`）——即口径上它是值安全的，不是定价臂。
// 4. 逐位等价闸门（本机跑 n=4096；本机对"逐位等价"有效）：两组独立随机输入 A/B，逐 256 行分 16 块
// 的全量 FNV，**16/16 块与基座逐位相同**；n=2048 一档由判题机免费通道逐位复核 ✓。
// 5. 为什么限在最外层：内层传送带节点（`rec<2048>`，ldc=4096）若同样连续存放，其 2h*h 连续区会
// 溢出自身 C 象限、压到兄弟节点的输出上 —— 本席实测该变体在 n=4096 上两组输入都逐位不同
// （判负），故用运行期判据 `ldc==2*h` 把改动限定在最外层那一个节点（也正因如此，n=2048 与
// n=4096 两档命中的是同一个节点，免费通道的读数可直接外推）。
// 6. 记账：−0.47% × 1860.379542 ms ≈ 8.7 ms（占比律，判题机 n=2048 直读），缺口 2.451784 ms。
//
// ===== REFERENCES =====
// [1] 本账号 saffah_cc_v41_agg1, 提交 #123974（`d4aj_` 席，1860.397775 ms，本账号在册最好）
// <https://duck.ac/submission/123974>
// 用途：**本提交的直接基座**。代码区与该件除本席下面那处定点改动外逐字节相同；
// 该件自身保留的 Credits 链原样保留在下方正文注释里，未作删改。
// [2] 在册账（同账号，本题）与封条：`d4ai_` 封条 **(2072)**（融合传送带 + 泵位归属）·
// `d4aj_` 封条 **(2076)**（爆发铺开判死 + 泵位时间表）。本席（`d4ak_`）入场读数：
// `mine 1860.397775`（#123974）· `T 1876.693695`（#122097，`saffah_codex_6a_agg3`）·
// 严支 `1857.927758` ⇒ 缺 **2.470017 ms (0.133%)**。
// [3] 在册封条中与本刀直接相关的两条：封条 **(2074)/(2073)** 所在册页的「本机结构性
// 失明」律（第 249 条：判题机 L3 = 6 MB ↔ 本机 39 MB，凡"工作集/争用/爆发"类改动
// 本机读数连符号都不可信）；以及第 248 条（相位内的拍数节省只兑现 ≈50% —— 本刀
// 针对的是**暴露在关键路径上的串行相位**，故按 1.0× 记账并在下方给出实测依据）。
//
// ===== 思路 =====
// **本发只做一件事：把本题唯一一个「完全暴露、无任何相位可重叠」的串行相位 ——
// 最后一个 rec<2048> 节点（root_last）的 h=1024 延迟 combine —— 按其真实数据依赖
// 劈成两段：M3 无关的三路输出（r11/r21/r22 + 各自的 root-u4 加法项）搬进**它自己
// 那最后一个子节点的叶相窗口**（复用该节点当时空着的 `pend2` 槽），端局只留
// `r12 = fl(fl(u4+m3)+p12)` 的一读三写。**
//
// 事实链（全部可用本地真引擎同源计数器复核，`work/d4ak__tmp/d4ak_gen_cnt.py`）：
// 1. ★ 本题「mikro 以外」的相位里，**唯一不被叶相掩蔽的一坨**就是它：在役件 #123974
// 的 `pending`/`pend2`/`root` 三个延迟 combine 逐 step 计数（n=4096 全跑）：
// pending 200704 step（12.85M 元素）: 叶内 198912 · 兜底 **1792**（0.9%）
// pend2 458752 step（ 7.34M 元素）: 叶内 393216 · 兜底 **65536**（14.3%）← 就是它
// root 131072 step（ 4.19M 元素）: 叶内 131072 · 兜底 **0**
// `combine()`（渴望路径）调用数 = **0**（n=4096 上它是死代码 ⇒ 本项定价 = 0）
// ⇒ `pend2_finish` 的 65536 step = 1.05M 元素 × (7 载入 + 4 u4 载入 + 4 存)/4 元素
// = **120 B/元素 = 126 MB**，本机 rdtsc 括号实测 **15.85e6 拍 = 全跑 0.278%**，
// 且**结构上不可重叠**（它是全程序最后一段，其后没有任何叶相调用）。
// 2. ★ 但它的四路输出里**只有一路**等最后那个子节点（wino 形）：
// r11 = m1+m2 (+p11) m1..m4 = 该节点的 C 四象限 = M1,M2,M3,M4
// r12 = (m1+m6+m5)+m3 (+p12) m3 = **最后那个子节点**的输出
// r21 = (m1+m6+m7)-m4 (+p21) m5,m6,m7 = 该节点的三块面板 M5,M6,M7
// r22 = (m1+m6+m7)+m5 (+p22) p* = root 的 u4 面板（root_step 就地写出的）
// 而 m1,m2,m4,m5,m6,m7 来自**先跑的那六个子节点**、`p*` 来自 root —— 后者按结构
// 在最后那个子节点开始前就已跑完（root 共 131072 step，其窗口 = 该节点的 5 个
// rec<1024> 子节点共 35 叶 ≈ 191k 次 kernel 调用，root 在第 131072 次调用即
// 第 4 个子节点的叶相里跑完）。
// 3. 于是 **M3 无关的 12 个访存/4 元素 + u4 的落地** 全搬进最后子节点的叶相窗口：
// 复用该节点当时**已排空**的 `pend2` 槽（它在 child B 之后的 `pend2_finish()` 已被
// 清空）与**现有**叶相钩子 ⇒ 热循环里**不新增分支**，只在既有钩子里把
// `pend2_step()` 换成分派（`part=false` 时逐字节等价于原路径）。
// u4 = (m1+m6)+m5 就地停在该节点**自己的第一块面板**里（该槽此刻已死：本发之后
// 该节点的三块面板再无消费者；root 的面板在另一块 `workspace`，不冲突）。
// 4. 端局只剩 `r12 = fl(fl(u4 + m3) + p12)`：**3 载入 + 1 存 / 4 元素 = 32 B/元素
// = 33.5 MB**（原 126 MB 的 27%）。
//
// 逐位等价（本件的最强闸门）：`chk = 50f4e2a76848b319` **与本席在役基座逐位相同**
// （生成器 `work/d4ak__tmp/d4ak_gen_knife.py` 逐行只做定点字符串替换；所有求和**保持原
// 操作数次序**：`r11=fl(p11+fl(m1+m2))`、`r21=fl(p21+fl(u3-m4))`、`r22=fl(p22+fl(u3+m5))`、
// `r12=fl(p12+fl(u4+m3))`、`u2=fl(m1+m6)`、`u3=fl(u2+m7)`、`u4=fl(u2+m5)` 一字未改）。
// 正确性护栏（三条，全部为保守检查）：partial 步只在 `!pending.active`（child 0 的
// h=512 延迟 combine 尚未把 c11/c21/c22 写完时不许读）且 `!root_pending.active`
// （root 的 u4 面板未完成时不许读）时推进；若叶相窗口不够（任何调度下），
// `pend2_partial_finish()` 先把残余 step 排空再进端局 ⇒ **值不依赖调度**。
//
// 本机实测（同进程相位括号，n=4096 全跑，`chk` 逐位不变）：
// 在役 #123974 端局/兜底 = 15.85e6 拍 · 本件 = partial 叶内 3.78e6 + 端局 3.42e6
// = 7.20e6 拍 ⇒ **暴露相位 −54%**；partial 的 16384 step 全部在叶相窗口内跑完
// （兜底 0 step），端局元素数 262144×4 = 1.05M 不变。
// ★ 诚实记账：本机跨批墙钟漂移 ±3~5%，**本机总分号不可作判据**（在册第 249 条：
// 判题机 L3 仅为本机 1/6.5，本相位在判题机上是纯 DRAM 面）；本件按「126 MB →
// 33.5 MB 的暴露字节」记账，判题机口径下 92.5 MB ÷ 5.9 B/核拍 ≈ 15.7e6 核拍 ≈ 4.3 ms。
//
// 合规：引用已列原作者 + 提交 URL + 所用内容；未改动任何算式/操作数次序/地址映射
// （只把同一个 combine 的同一批加法**换个相位做**）；未加特例；未触碰任何判据；
// 无 /tmp 产物；本件为真拼写差异件（对基座 `cmp` 二进制不同、源码为定点插入）。
// =====================================================================
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/121415: retained frontier-gated FP64 matrix engine and all inherited citations below.
* Idea: Express the two carried-packing source addresses as constant displacements from the packing destination. Advance only the destination pointer inside the FMA loop, replacing two pointer additions per eight-K step with x86 indexed loads. Values, FMA order and packing schedule are unchanged.
* Purpose: Reduce integer updates and live output operands in the carried-packing microkernel.
*/
#include <cstdint>
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/121243,
reused the matrix engine and inherited source credits below. No independent
license declared.
Idea: Gate progressive reconstruction once at each kernel call. A row frontier
is a multiple of 128 elements, and paired 64-element steps also start at that
alignment, so a permitted pair cannot cross the frontier. Use unchecked steps
inside that gate; stalled calls skip both step checks and the prefetch check.
Purpose: Reduce repeated bounds checks while reconstruction waits for more rows.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/121181,
reused the matrix engine and all inherited source references below.
No independent license declared.
Idea: Specialize recursion for the known 4096, 2048, 1024 and 512 sizes. Each
level can discard unreachable size branches and use constant transform strides
and extents, preserving the same product order and reconstruction schedule.
Purpose: Remove dynamic size arithmetic and dispatch above the matrix leaves.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/121167,
reused the matrix engine and all inherited references below. No independent
license declared.
Idea: Every recursive multiplication in this fixed 4096-square task terminates
at a 512-square leaf. Compile that leaf extent as a constant to simplify packing
bounds, panel strides, tail selection and streaming-store conditions.
Purpose: Remove general-size arithmetic from the remaining leaf hot loops.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/121142,
reused the matrix engine and all inherited references below. No independent
license declared.
Idea: Specialize leaf reconstruction independently for each of the three
levels active at entry. Deferred tasks can retire but cannot start within a
leaf; omit checks for the levels that were absent when this leaf began.
Purpose: Extend the measured benefit of inactive-task specialization to leaves
with only one or two active reconstruction levels.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/120888,
reused the current pipelined matrix engine and all inherited source credits below.
No independent license declared.
Idea: Select a leaf implementation once based on whether deferred reconstruction
is active. No reconstruction task can begin inside a leaf, so a leaf entering
with none can omit every per-microkernel task check and row-frontier update.
The other implementation retains all existing task progression and ordering.
Purpose: Remove repeated inactive-task dispatch from most microkernel calls.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/120615,
reused the current matrix engine and plain-kernel software B-load pipeline.
All earlier source attributions remain below; no independent license declared.
Idea: Apply the one-step B-load pipeline to the A-packing microkernel too.
Its spare register is only needed for packing after the last inner FMA, so
alternating that register with the current B vector preserves all values.
Purpose: Hide B load latency in the remaining carried-packing kernel calls.
*/
/* References: pdoom https://duck.ac/submission/118506: reused FP64
recursive matrix engine. saffah_cc_v41_agg1 https://duck.ac/submission/119315:
reused plain-kernel B-load software pipeline and deferred prefetch scheduling.
saffah_codex_6a_agg3 https://duck.ac/submission/120547: B-panel reuse and merge64.
All inherited source credits retained below; no independent licenses declared.
Idea: Apply the one-step-ahead B-load scheduling to the larger matrix kernel,
using the fixed-size output path and preserving each FMA accumulation order.
Purpose: Hide B load-to-use latency across more than100000 microkernel calls.
*/
/* References:
[1] saffah_codex_6a_agg3, https://duck.ac/submission/120521.
Reused the complete FP64 matrix engine; all inherited source citations remain
below. No separate license is declared on the public source.
Idea: Reconstruct64 positions per bottom merge step instead of32. Both divide
the512-wide rows, preserving frontiers and operation grouping while amortizing
step checks and cursor updates across twice as much contiguous memory traffic.
Purpose: Measure lower reconstruction bookkeeping in the fused compute schedule.
*/
/* References:
[1] pdoom, https://duck.ac/submission/118506.
Reused the4k FP64 engine and its deferred root reconstruction, preserving its
complete inherited attribution list below. No separate license is declared.
[2] saffah_cc_v41_agg1, https://duck.ac/submission/119315.
Copied the1k fused packed-B helpers and two-reuse leaf ordering.
[3] saffah_codex_6a_agg3, https://duck.ac/submission/120514.
Reused our validated2k integration of those helpers; no separate license.
Idea: Apply two packed-B reuse positions to every bottom group of the4k
recursion, removing repeated source traversal when a pure panel is still packed.
Purpose: Measure the accumulated traffic reduction at the largest FP64 size.
*/
/*
Credits for inherited implementation and ideas:
- saffah_cc_v41_agg1, https://duck.ac/submission/105696: immediate engine
base; safely extend the upper deferred-combine window across siblings.
- saffah_codex_6s_agg2, https://duck.ac/submission/104821: preceding engine,
compiler model and upper workspace layout, including its inherited work.
- saffah_cc_v41_agg1, https://duck.ac/submission/101193: two-level deferred
reconstruction, product order and panel-layout improvements.
- saffah_cc_v41_agg1, https://duck.ac/submission/103010: produced-row
frontier and early reconstruction; saffah_codex_6s_agg2,
https://duck.ac/submission/103224: its4k integration and leaf ordering.
- saffah_codex_6s_agg2, https://duck.ac/submission/103194: eight-K plain
microkernel; https://duck.ac/submission/103275: power-of-two addressing;
https://duck.ac/submission/104809: leaf scratch-plane layout.
- saffah_cc_v41_agg1, https://duck.ac/submission/100776: dedicated two-row
tail; https://duck.ac/submission/100511 and https://duck.ac/submission/100475:
carry next-A packing in the FMA loop, double-buffer A and grouped B packing.
- saffah_cc_v41_260924, https://duck.ac/submission/99846: fused operand sums
and packing, one-pass sums and reconstruction.
- saffah_cc_v41_agg1, https://duck.ac/submission/99875 and
https://duck.ac/submission/99950: non-temporal intermediate stores.
- saffah_cc_v41_260924, https://duck.ac/submission/99141, and
saffah_cc_v41_agg1, https://duck.ac/submission/99167:6x8 AVX2 blocking,
explicit vector memory instructions, prefetch and packing-layout baseline.
- Earlier pdoom work: hybrid Winograd/Strassen recurrence, output-backed
scratch placement, delayed reconstruction, merged32-row tail and2x16
kernel, and page-spacing experiments.
New changes in this version:
1. Specialize the full6x8 kernels and hot upper reconstruction geometry for
the fixed4096 problem. Two-row output has its separate2x16 kernel.
2. The root's last product M3 only contributes to C12. Once the preceding
M1 is ready, reconstruct C11,C21,C22 and the C12 base (M1+M6)+M5 during
M3's remaining five1024 subproblems. Retired M5 storage holds that base.
3. Add the base while reconstructing M3's four quadrants. The root task
processes32 elements per step only when bottom reconstruction is inactive.
This moves320MiB of logical root traffic into the last multiplication
subtree; total nominal traffic and arithmetic count do not decrease.
All arithmetic remains binary64, with the original expression grouping.
No input quantization. Specialized to the problem's n=4096 contract.
*/
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing")
#include <sys/mman.h>
#pragma GCC target("avx2,fma,arch=skylake")
#include <immintrin.h>
#ifndef B2M_NO_BARA
#define B2M_NO_BARA 1
#endif
#define MD9_T2BASE 1
#include <stddef.h>
#ifndef BASE
#define BASE 512
#endif
#ifndef MC
#define MC 30
#endif
#ifndef NTMIN
#define NTMIN 512
#endif
#ifndef UPPER_PAGE_PAD
#define UPPER_PAGE_PAD 15872
#endif
#ifndef LEAF_PAD
#define LEAF_PAD 1152
#endif
#ifndef SUMNT
#define SUMNT 1024
#endif
static constexpr int MR=6,NR=8;
alignas(4096) static double workspace[48*1024*1024];
alignas(64) static double ap[(MC+MR)*(BASE+24)],ap2[(MC+MR)*(BASE+24)],bp[BASE*(BASE+24)];
/* m2s/md5 fused-pack support: a zero row (for pad rows and for pure copies) and the
two sign masks (+0.0 / -0.0) used to apply Strassen's operand sign in-register. */
alignas(64) static double pk_zero[BASE+8];
alignas(64) static double pk_mp[4]={0.0,0.0,0.0,0.0}, pk_mn[4]={-0.0,-0.0,-0.0,-0.0};
struct Mat { const double *x; int ld; const double *y=nullptr; int ly=0; int sign=0; };
static inline __m256d ld(Mat a,int i,int j) {
__m256d x=_mm256_load_pd(a.x+(size_t)i*a.ld+j);
if(a.sign) {
__m256d y=_mm256_load_pd(a.y+(size_t)i*a.ly+j);
x=a.sign>0?_mm256_add_pd(x,y):_mm256_sub_pd(x,y);
}
return x;
}
static inline double scalar(Mat a,int i,int k) {
double x=a.x[(size_t)i*a.ld+k];
if(a.sign)x+=a.sign>0?a.y[(size_t)i*a.ly+k]:-a.y[(size_t)i*a.ly+k];
return x;
}
template<int SG> static __attribute__((always_inline)) inline void kernel(int n,const double *a,const double *b,double *c,int ldc,int rows, const double *pks,const double *pky,double *pkd) {
double *out=c;
ptrdiff_t stride=(ptrdiff_t)ldc*8;
int k=n;
pks=(const double*)((uintptr_t)pks-(uintptr_t)pkd);
pky=(const double*)((uintptr_t)pky-(uintptr_t)pkd);
asm volatile (
"vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
"vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
"vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
"vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
"vxorpd %%ymm4, %%ymm4, %%ymm4\n\t"
"vxorpd %%ymm5, %%ymm5, %%ymm5\n\t"
"vxorpd %%ymm6, %%ymm6, %%ymm6\n\t"
"vxorpd %%ymm7, %%ymm7, %%ymm7\n\t"
"vxorpd %%ymm8, %%ymm8, %%ymm8\n\t"
"vxorpd %%ymm9, %%ymm9, %%ymm9\n\t"
"vxorpd %%ymm10, %%ymm10, %%ymm10\n\t"
"vxorpd %%ymm11, %%ymm11, %%ymm11\n\t"
"vmovapd 0(%[b]), %%ymm12\n\t"
".p2align 6\n\t"
"1:\n\t"
"vmovapd 64(%[b]), %%ymm15\n\t"
"vmovapd 32(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 128(%[b]), %%ymm12\n\t"
"vmovapd 96(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 512(%[b])\n\t"
"prefetcht0 576(%[b])\n\t"
"prefetcht0 640(%[b])\n\t"
"prefetcht0 704(%[b])\n\t"
"vmovapd 192(%[b]), %%ymm15\n\t"
"vmovapd 160(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 256(%[b]), %%ymm12\n\t"
"vmovapd 224(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 768(%[b])\n\t"
"prefetcht0 832(%[b])\n\t"
"prefetcht0 896(%[b])\n\t"
"prefetcht0 960(%[b])\n\t"
"vmovapd 320(%[b]), %%ymm15\n\t"
"vmovapd 288(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 384(%[b]), %%ymm12\n\t"
"vmovapd 352(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 448(%[b]), %%ymm15\n\t"
"vmovapd 416(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 512(%[b]), %%ymm12\n\t"
"vmovapd 480(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovupd 0(%[pkd],%[pks]), %%ymm15\n\t"
".if %c[sg] == 1\n\t"
"vaddpd 0(%[pkd],%[pky]), %%ymm15, %%ymm15\n\t"
".elseif %c[sg] == -1\n\t"
"vsubpd 0(%[pkd],%[pky]), %%ymm15, %%ymm15\n\t"
".endif\n\t"
"vmovupd %%ymm15, 0(%[pkd])\n\t"
"add $32, %[pkd]\n\t"
"add $64, %[a]\n\t"
"add $512, %[b]\n\t"
"sub $8, %[k]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, (%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm2, (%[out])\n\t"
"vmovntpd %%ymm3, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm4, (%[out])\n\t"
"vmovntpd %%ymm5, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm6, (%[out])\n\t"
"vmovntpd %%ymm7, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm8, (%[out])\n\t"
"vmovntpd %%ymm9, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm10, (%[out])\n\t"
"vmovntpd %%ymm11, 32(%[out])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k), [pkd] "+&r"(pkd)
: [pks] "r"(pks), [pky] "r"(pky), [sg] "i"(SG), [stride] "r"(stride), [o0] "i"(0*(BASE+24)*8), [o1] "i"(1*(BASE+24)*8), [o2] "i"(2*(BASE+24)*8), [o3] "i"(3*(BASE+24)*8), [o4] "i"(4*(BASE+24)*8), [o5] "i"(5*(BASE+24)*8)
: "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
template<int S> static __attribute__((noinline)) void kpA(int n,const double *a,const double *b,double *c,int ldc,int rows, const double *L0,const double *L1,const double *L2,double *D0,double *D1) {
double *out=c;
ptrdiff_t stride=(ptrdiff_t)ldc*8;
int k=n;
ptrdiff_t idx=0;
asm volatile (
"vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
"vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
"vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
"vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
"vxorpd %%ymm4, %%ymm4, %%ymm4\n\t"
"vxorpd %%ymm5, %%ymm5, %%ymm5\n\t"
"vxorpd %%ymm6, %%ymm6, %%ymm6\n\t"
"vxorpd %%ymm7, %%ymm7, %%ymm7\n\t"
"vxorpd %%ymm8, %%ymm8, %%ymm8\n\t"
"vxorpd %%ymm9, %%ymm9, %%ymm9\n\t"
"vxorpd %%ymm10, %%ymm10, %%ymm10\n\t"
"vxorpd %%ymm11, %%ymm11, %%ymm11\n\t"
"vmovapd 0(%[b]), %%ymm12\n\t"
".p2align 6\n\t"
"1:\n\t"
"vmovapd 64(%[b]), %%ymm15\n\t"
"vmovapd 32(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 128(%[b]), %%ymm12\n\t"
"vmovapd 96(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 512(%[b])\n\t"
"prefetcht0 576(%[b])\n\t"
"prefetcht0 640(%[b])\n\t"
"prefetcht0 704(%[b])\n\t"
"vmovapd 192(%[b]), %%ymm15\n\t"
"vmovapd 160(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 256(%[b]), %%ymm12\n\t"
"vmovapd 224(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 768(%[b])\n\t"
"prefetcht0 832(%[b])\n\t"
"prefetcht0 896(%[b])\n\t"
"prefetcht0 960(%[b])\n\t"
"vmovapd 320(%[b]), %%ymm15\n\t"
"vmovapd 288(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 384(%[b]), %%ymm12\n\t"
"vmovapd 352(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 448(%[b]), %%ymm15\n\t"
"vmovapd 416(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 512(%[b]), %%ymm12\n\t"
"vmovapd 480(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovupd (%[L0],%[idx]), %%ymm14\n\t"
"vmovupd (%[L2],%[idx]), %%ymm15\n\t"
"vsubpd %%ymm14, %%ymm15, %%ymm15\n\t"
"vmovntpd %%ymm15, (%[D1],%[idx])\n\t"
"vmovupd (%[L1],%[idx]), %%ymm15\n\t"
".if %c[s] == 1\n\t"
"vaddpd %%ymm15, %%ymm14, %%ymm14\n\t"
".else\n\t"
"vsubpd %%ymm15, %%ymm14, %%ymm14\n\t"
".endif\n\t"
"vmovupd %%ymm14, (%[D0],%[idx])\n\t"
"add $32, %[idx]\n\t"
"add $64, %[a]\n\t"
"add $512, %[b]\n\t"
"sub $8, %[k]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, (%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm2, (%[out])\n\t"
"vmovntpd %%ymm3, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm4, (%[out])\n\t"
"vmovntpd %%ymm5, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm6, (%[out])\n\t"
"vmovntpd %%ymm7, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm8, (%[out])\n\t"
"vmovntpd %%ymm9, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm10, (%[out])\n\t"
"vmovntpd %%ymm11, 32(%[out])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k), [idx] "+&r"(idx)
: [L0] "r"(L0), [L1] "r"(L1), [L2] "r"(L2), [D0] "r"(D0), [D1] "r"(D1), [s] "i"(S), [stride] "r"(stride), [o0] "i"(0*(BASE+24)*8), [o1] "i"(1*(BASE+24)*8), [o2] "i"(2*(BASE+24)*8), [o3] "i"(3*(BASE+24)*8), [o4] "i"(4*(BASE+24)*8), [o5] "i"(5*(BASE+24)*8)
: "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
template<int S> static __attribute__((noinline)) void kpB(int n,const double *a,const double *b,double *c,int ldc,int rows, const double *L0,const double *L1,const double *L2,double *D0,double *D1) {
double *out=c;
ptrdiff_t stride=(ptrdiff_t)ldc*8;
int k=n;
ptrdiff_t idx=0;
asm volatile (
"vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
"vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
"vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
"vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
"vxorpd %%ymm4, %%ymm4, %%ymm4\n\t"
"vxorpd %%ymm5, %%ymm5, %%ymm5\n\t"
"vxorpd %%ymm6, %%ymm6, %%ymm6\n\t"
"vxorpd %%ymm7, %%ymm7, %%ymm7\n\t"
"vxorpd %%ymm8, %%ymm8, %%ymm8\n\t"
"vxorpd %%ymm9, %%ymm9, %%ymm9\n\t"
"vxorpd %%ymm10, %%ymm10, %%ymm10\n\t"
"vxorpd %%ymm11, %%ymm11, %%ymm11\n\t"
"vmovapd 0(%[b]), %%ymm12\n\t"
".p2align 6\n\t"
"1:\n\t"
"vmovapd 64(%[b]), %%ymm15\n\t"
"vmovapd 32(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 128(%[b]), %%ymm12\n\t"
"vmovapd 96(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 512(%[b])\n\t"
"prefetcht0 576(%[b])\n\t"
"prefetcht0 640(%[b])\n\t"
"prefetcht0 704(%[b])\n\t"
"vmovapd 192(%[b]), %%ymm15\n\t"
"vmovapd 160(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 256(%[b]), %%ymm12\n\t"
"vmovapd 224(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 768(%[b])\n\t"
"prefetcht0 832(%[b])\n\t"
"prefetcht0 896(%[b])\n\t"
"prefetcht0 960(%[b])\n\t"
"vmovapd 320(%[b]), %%ymm15\n\t"
"vmovapd 288(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 384(%[b]), %%ymm12\n\t"
"vmovapd 352(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 448(%[b]), %%ymm15\n\t"
"vmovapd 416(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 512(%[b]), %%ymm12\n\t"
"vmovapd 480(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovupd (%[L0],%[idx]), %%ymm14\n\t"
"vmovupd (%[L1],%[idx]), %%ymm15\n\t"
"vsubpd %%ymm15, %%ymm14, %%ymm14\n\t"
"vmovntpd %%ymm14, (%[D0],%[idx])\n\t"
"vmovupd (%[L2],%[idx]), %%ymm15\n\t"
".if %c[s] == 0\n\t"
"vsubpd %%ymm14, %%ymm15, %%ymm15\n\t"
".else\n\t"
"vsubpd %%ymm15, %%ymm14, %%ymm15\n\t"
".endif\n\t"
"vmovntpd %%ymm15, (%[D1],%[idx])\n\t"
"add $32, %[idx]\n\t"
"add $64, %[a]\n\t"
"add $512, %[b]\n\t"
"sub $8, %[k]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, (%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm2, (%[out])\n\t"
"vmovntpd %%ymm3, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm4, (%[out])\n\t"
"vmovntpd %%ymm5, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm6, (%[out])\n\t"
"vmovntpd %%ymm7, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm8, (%[out])\n\t"
"vmovntpd %%ymm9, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm10, (%[out])\n\t"
"vmovntpd %%ymm11, 32(%[out])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k), [idx] "+&r"(idx)
: [L0] "r"(L0), [L1] "r"(L1), [L2] "r"(L2), [D0] "r"(D0), [D1] "r"(D1), [s] "i"(S), [stride] "r"(stride), [o0] "i"(0*(BASE+24)*8), [o1] "i"(1*(BASE+24)*8), [o2] "i"(2*(BASE+24)*8), [o3] "i"(3*(BASE+24)*8), [o4] "i"(4*(BASE+24)*8), [o5] "i"(5*(BASE+24)*8)
: "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
static __attribute__((always_inline)) inline void kernel_plain(int n,const double *a,const double *b,double *c,int ldc,int rows) {
double *out=c;
ptrdiff_t stride=(ptrdiff_t)ldc*8;
int k=n;
asm volatile (
"vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
"vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
"vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
"vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
"vxorpd %%ymm4, %%ymm4, %%ymm4\n\t"
"vxorpd %%ymm5, %%ymm5, %%ymm5\n\t"
"vxorpd %%ymm6, %%ymm6, %%ymm6\n\t"
"vxorpd %%ymm7, %%ymm7, %%ymm7\n\t"
"vxorpd %%ymm8, %%ymm8, %%ymm8\n\t"
"vxorpd %%ymm9, %%ymm9, %%ymm9\n\t"
"vxorpd %%ymm10, %%ymm10, %%ymm10\n\t"
"vxorpd %%ymm11, %%ymm11, %%ymm11\n\t"
"vmovapd 0(%[b]), %%ymm12\n\t"
".p2align 6\n\t"
"1:\n\t"
"vmovapd 64(%[b]), %%ymm15\n\t"
"vmovapd 32(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 128(%[b]), %%ymm12\n\t"
"vmovapd 96(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 512(%[b])\n\t"
"prefetcht0 576(%[b])\n\t"
"prefetcht0 640(%[b])\n\t"
"prefetcht0 704(%[b])\n\t"
"vmovapd 192(%[b]), %%ymm15\n\t"
"vmovapd 160(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 256(%[b]), %%ymm12\n\t"
"vmovapd 224(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 768(%[b])\n\t"
"prefetcht0 832(%[b])\n\t"
"prefetcht0 896(%[b])\n\t"
"prefetcht0 960(%[b])\n\t"
"vmovapd 320(%[b]), %%ymm15\n\t"
"vmovapd 288(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 384(%[b]), %%ymm12\n\t"
"vmovapd 352(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 448(%[b]), %%ymm15\n\t"
"vmovapd 416(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 512(%[b]), %%ymm12\n\t"
"vmovapd 480(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"add $64, %[a]\n\t"
"add $512, %[b]\n\t"
"sub $8, %[k]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, (%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm2, (%[out])\n\t"
"vmovntpd %%ymm3, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm4, (%[out])\n\t"
"vmovntpd %%ymm5, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm6, (%[out])\n\t"
"vmovntpd %%ymm7, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm8, (%[out])\n\t"
"vmovntpd %%ymm9, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm10, (%[out])\n\t"
"vmovntpd %%ymm11, 32(%[out])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k)
: [stride] "r"(stride), [o0] "i"(0*(BASE+24)*8), [o1] "i"(1*(BASE+24)*8), [o2] "i"(2*(BASE+24)*8), [o3] "i"(3*(BASE+24)*8), [o4] "i"(4*(BASE+24)*8), [o5] "i"(5*(BASE+24)*8)
: "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
/* ===== mdA_dp2: deferred h=1024 combine (Winograd form, exact same arithmetic) ===== */
struct Pend2 { double *c; int ldc; const double *m; int h; int pz; int pos; bool active; const double *add=nullptr; bool part=false; };
static Pend2 pend2={0,0,0,0,0,0,false};
static constexpr int STEP2=16;
template<int H=1024> static inline void pend2_prefetch(void){
if(!pend2.active)return;
if(pend2.part){ /* d4ak_: partial pass = 64 elements/step + u4 park */
int i=pend2.pos>>10, j=pend2.pos&1023;
const double *c=pend2.c+(size_t)i*pend2.ldc+j;
const double *mm=pend2.m+pend2.pos;
const double *pp=pend2.add+(size_t)i*2048+j;
for(int t=0;t<64;t+=8) {
_mm_prefetch((const char*)(c+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+1024+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)1024*pend2.ldc+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)1024*pend2.ldc+1024+t),_MM_HINT_T1);
_mm_prefetch((const char*)(mm+t),_MM_HINT_T1);
_mm_prefetch((const char*)(mm+(1024*1024+UPPER_PAGE_PAD)+t),_MM_HINT_T1);
_mm_prefetch((const char*)(mm+2*(1024*1024+UPPER_PAGE_PAD)+t),_MM_HINT_T1);
_mm_prefetch((const char*)(pp+t),_MM_HINT_T1);
_mm_prefetch((const char*)(pp+(size_t)1024*2048+t),_MM_HINT_T1);
_mm_prefetch((const char*)(pp+(size_t)1024*2048+1024+t),_MM_HINT_T1);
}
return;
}
int i=pend2.pos >> __builtin_ctz((unsigned)H), j=pend2.pos & (H-1);
const double *c=pend2.c+(size_t)i*pend2.ldc+j;
const double *m=pend2.m+pend2.pos;
for(int t=0;t<STEP2;t+=8) {
_mm_prefetch((const char*)(c+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+H+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)H*pend2.ldc+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)H*pend2.ldc+H+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+(H*H+UPPER_PAGE_PAD)+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+2*(H*H+UPPER_PAGE_PAD)+t),_MM_HINT_T1);
}
}
template<int H=1024,bool ADD=false> static inline void pend2_step(void){
int i=pend2.pos >> __builtin_ctz((unsigned)H), j=pend2.pos & (H-1);
for(int t=0;t<STEP2;t+=4) {
double *c11=pend2.c+(size_t)i*pend2.ldc+j+t,*c12=c11+H,
*c21=c11+(size_t)H*pend2.ldc,*c22=c21+H;
const double *m=pend2.m+pend2.pos+t;
__m256d m1=_mm256_load_pd(c11),m2=_mm256_load_pd(c21),
m3=_mm256_load_pd(c12),m4=_mm256_load_pd(c22),
m5=_mm256_load_pd(m),m6=_mm256_load_pd(m+(H*H+UPPER_PAGE_PAD)),m7=_mm256_load_pd(m+2*(H*H+UPPER_PAGE_PAD));
__m256d u2=_mm256_add_pd(m1,m6),u3=_mm256_add_pd(u2,m7),u4=_mm256_add_pd(u2,m5);
__m256d r11=_mm256_add_pd(m1,m2),r12=_mm256_add_pd(u4,m3),
r21=_mm256_sub_pd(u3,m4),r22=_mm256_add_pd(u3,m5);
if(ADD) {
const double *p=pend2.add+(size_t)i*2048+j+t;
r11=_mm256_add_pd(_mm256_load_pd(p),r11);
r12=_mm256_add_pd(_mm256_load_pd(p+H),r12);
r21=_mm256_add_pd(_mm256_load_pd(p+(size_t)H*2048),r21);
r22=_mm256_add_pd(_mm256_load_pd(p+(size_t)H*2048+H),r22);
}
_mm256_store_pd(c11,r11);_mm256_store_pd(c12,r12);
_mm256_store_pd(c21,r21);_mm256_store_pd(c22,r22);
}
pend2.pos+=STEP2;
if(pend2.pos>=H*H)pend2.active=false;
}
static void pend2_finish(void){
if(pend2.add) {while(pend2.active)pend2_step<1024,true>();}
else if(pend2.h==1024) {while(pend2.active)pend2_step<1024>();}
else {while(pend2.active)pend2_step<2048>();}
}
struct Pending { double *c; int ldc; const double *m; int pos=0; bool active=false; int rowlim=0; bool prog=false; };
static Pending pending;
alignas(4096) static double pending_work[3*(BASE*BASE+LEAF_PAD)+LEAF_PAD];
static unsigned pending_tick=0;
static constexpr int STEP=64, PERIOD=1;
static inline void pending_prefetch() {
if(!pending.active)return;
if(pending.prog && pending.pos/BASE>=pending.rowlim)return;
int i=pending.pos/BASE,j=pending.pos%BASE;
const double *c=pending.c+(size_t)i*pending.ldc+j;
const double *m=pending.m+pending.pos;
for(int t=0;t<STEP;t+=8) {
_mm_prefetch((const char*)(c+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+BASE+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)BASE*pending.ldc+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)BASE*pending.ldc+BASE+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+(BASE*BASE+LEAF_PAD)+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+2*(BASE*BASE+LEAF_PAD)+t),_MM_HINT_T1);
}
}
static inline void pending_step_unchecked() {
int i=pending.pos/BASE,j=pending.pos%BASE;
for(int t=0;t<STEP;t+=4) {
double *c11=pending.c+(size_t)i*pending.ldc+j+t,
*c12=c11+BASE,*c21=c11+(size_t)BASE*pending.ldc,*c22=c21+BASE;
const double *m=pending.m+pending.pos+t;
__m256d m1=_mm256_load_pd(c11),m2=_mm256_load_pd(c21),
m3=_mm256_load_pd(c12),m4=_mm256_load_pd(c22),
m5=_mm256_load_pd(m),m6=_mm256_load_pd(m+(BASE*BASE+LEAF_PAD)),m7=_mm256_load_pd(m+2*(BASE*BASE+LEAF_PAD));
__m256d r11=_mm256_add_pd(_mm256_sub_pd(_mm256_add_pd(m1,m4),m5),m7),
r12=_mm256_add_pd(m3,m5),r21=_mm256_add_pd(m2,m4),
r22=_mm256_add_pd(_mm256_add_pd(_mm256_sub_pd(m1,m2),m3),m6);
_mm256_store_pd(c11,r11);_mm256_store_pd(c12,r12);
_mm256_store_pd(c21,r21);_mm256_store_pd(c22,r22);
}
pending.pos+=STEP;
if(pending.pos==BASE*BASE)pending.active=false;
}
static inline void pending_step(){if(pending.prog && pending.pos/BASE>=pending.rowlim)return;pending_step_unchecked();}
static void pending_finish() {
pending.prog=false;
while(pending.active)pending_step();
}
static __attribute__((always_inline)) inline void kernel2(int n,const double *a,const double *b,double *c,int ldc) {
int k=n;
double *out=c;
ptrdiff_t stride=(ptrdiff_t)ldc*8;
asm volatile (
"vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
"vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
"vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
"vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
".p2align 6\n\t" "1:\n\t"
"prefetcht0 512(%[b])\n\t"
"prefetcht0 576(%[b])\n\t"
"prefetcht0 640(%[b])\n\t"
"prefetcht0 704(%[b])\n\t"
"vmovapd 0(%[b]), %%ymm12\n\t"
"vmovapd 32(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vmovapd 64(%[b]), %%ymm12\n\t"
"vmovapd 96(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vmovapd 128(%[b]), %%ymm12\n\t"
"vmovapd 160(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vmovapd 192(%[b]), %%ymm12\n\t"
"vmovapd 224(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"add $32, %[a]\n\t"
"add $256, %[b]\n\t"
"sub $4, %[k]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, (%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm2, (%[out])\n\t"
"vmovntpd %%ymm3, 32(%[out])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k)
: [stride] "r"(stride), [o0] "i"(0*(BASE+24)*8), [o1] "i"(1*(BASE+24)*8)
: "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm12","ymm13","ymm14");
}
static __attribute__((noinline)) void kernel_tail16(int n,const double *a,const double *b,double *out,int ldc) {
const double *b1=b+8*(n+1);
ptrdiff_t stride=(ptrdiff_t)ldc*8;
int kk=n;
asm volatile (
"vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
"vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
"vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
"vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
"vxorpd %%ymm4, %%ymm4, %%ymm4\n\t"
"vxorpd %%ymm5, %%ymm5, %%ymm5\n\t"
"vxorpd %%ymm6, %%ymm6, %%ymm6\n\t"
"vxorpd %%ymm7, %%ymm7, %%ymm7\n\t"
".p2align 6\n\t"
"1:\n\t"
"vbroadcastsd %c[o0]+0(%[a]), %%ymm12\n\t"
"vbroadcastsd %c[o1]+0(%[a]), %%ymm13\n\t"
"vmovapd 0(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm0\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm4\n\t"
"vmovapd 32(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm5\n\t"
"vmovapd 0(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm6\n\t"
"vmovapd 32(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"vbroadcastsd %c[o0]+8(%[a]), %%ymm12\n\t"
"vbroadcastsd %c[o1]+8(%[a]), %%ymm13\n\t"
"vmovapd 64(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm0\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm4\n\t"
"vmovapd 96(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm5\n\t"
"vmovapd 64(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm6\n\t"
"vmovapd 96(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"vbroadcastsd %c[o0]+16(%[a]), %%ymm12\n\t"
"vbroadcastsd %c[o1]+16(%[a]), %%ymm13\n\t"
"vmovapd 128(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm0\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm4\n\t"
"vmovapd 160(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm5\n\t"
"vmovapd 128(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm6\n\t"
"vmovapd 160(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"vbroadcastsd %c[o0]+24(%[a]), %%ymm12\n\t"
"vbroadcastsd %c[o1]+24(%[a]), %%ymm13\n\t"
"vmovapd 192(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm0\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm4\n\t"
"vmovapd 224(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm5\n\t"
"vmovapd 192(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm6\n\t"
"vmovapd 224(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"add $32, %[a]\n\t"
"add $256, %[b]\n\t"
"add $256, %[b1]\n\t"
"sub $4, %[kk]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, 0(%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"vmovntpd %%ymm2, 64(%[out])\n\t"
"vmovntpd %%ymm3, 96(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm4, 0(%[out])\n\t"
"vmovntpd %%ymm5, 32(%[out])\n\t"
"vmovntpd %%ymm6, 64(%[out])\n\t"
"vmovntpd %%ymm7, 96(%[out])\n\t"
: [a] "+&r"(a),[b] "+&r"(b),[b1] "+&r"(b1),[out] "+&r"(out),[kk] "+&r"(kk)
: [stride] "r"(stride),[o0] "i"(0),[o1] "i"((BASE+24)*8)
: "memory","cc","ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
struct RootPending { double*c;int ldc;double*m;int pos;bool active; };
static RootPending root_pending={nullptr,0,nullptr,0,false};
static constexpr int ROOT_STEP=32;
static inline void root_prefetch() {
if(!root_pending.active)return;
int i=root_pending.pos>>11,j=root_pending.pos&2047;
const double*c=root_pending.c+(size_t)i*root_pending.ldc+j;
const double*m=root_pending.m+root_pending.pos;
constexpr size_t PZ=2048*2048+UPPER_PAGE_PAD;
for(int t=0;t<ROOT_STEP;t+=8){
_mm_prefetch((const char*)(c+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)2048*root_pending.ldc+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)2048*root_pending.ldc+2048+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+PZ+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+2*PZ+t),_MM_HINT_T1);
}
}
static inline void root_step() {
int i=root_pending.pos>>11,j=root_pending.pos&2047;
constexpr size_t PZ=2048*2048+UPPER_PAGE_PAD;
for(int t=0;t<ROOT_STEP;t+=4){
double*c11=root_pending.c+(size_t)i*root_pending.ldc+j+t;
double*c21=c11+(size_t)2048*root_pending.ldc,*c22=c21+2048;
double*m=root_pending.m+root_pending.pos+t;
__m256d m1=_mm256_load_pd(c11),m2=_mm256_load_pd(c21),m4=_mm256_load_pd(c22);
__m256d m5=_mm256_load_pd(m),m6=_mm256_load_pd(m+PZ),m7=_mm256_load_pd(m+2*PZ);
__m256d u2=_mm256_add_pd(m1,m6),u3=_mm256_add_pd(u2,m7),u4=_mm256_add_pd(u2,m5);
_mm256_store_pd(c11,_mm256_add_pd(m1,m2));
_mm256_store_pd(c21,_mm256_sub_pd(u3,m4));
_mm256_store_pd(c22,_mm256_add_pd(u3,m5));
_mm256_store_pd(m,u4);
}
root_pending.pos+=ROOT_STEP;
if(root_pending.pos==2048*2048)root_pending.active=false;
}
static void root_finish(){while(root_pending.active)root_step();}
/* ===== d4ak_: split of the last node's h=1024 combine (see file header) ===== */
static void pend2_partial_step(void) {
if(pending.active)return; /* child-0's deferred h=512 combine still writes c11/c21/c22 */
if(root_pending.active)return; /* pend2.add (the root u4 panel) must be complete */
int i=pend2.pos>>10, j=pend2.pos&1023;
for(int t=0;t<64;t+=4) {
double *c11=pend2.c+(size_t)i*pend2.ldc+j+t,*c21=c11+(size_t)1024*pend2.ldc,*c22=c21+1024;
const double *mm=pend2.m+pend2.pos+t;
const double *pp=pend2.add+(size_t)i*2048+j+t;
__m256d m1=_mm256_load_pd(c11),m2=_mm256_load_pd(c21),m4=_mm256_load_pd(c22);
__m256d m5=_mm256_load_pd(mm),m6=_mm256_load_pd(mm+(1024*1024+UPPER_PAGE_PAD)),
m7=_mm256_load_pd(mm+2*(1024*1024+UPPER_PAGE_PAD));
__m256d p11=_mm256_load_pd(pp),p21=_mm256_load_pd(pp+(size_t)1024*2048),
p22=_mm256_load_pd(pp+(size_t)1024*2048+1024);
__m256d u2=_mm256_add_pd(m1,m6),u3=_mm256_add_pd(u2,m7),u4=_mm256_add_pd(u2,m5);
__m256d r11=_mm256_add_pd(m1,m2),r21=_mm256_sub_pd(u3,m4),r22=_mm256_add_pd(u3,m5);
r11=_mm256_add_pd(p11,r11); /* operand order == original ADD block */
r21=_mm256_add_pd(p21,r21);
r22=_mm256_add_pd(p22,r22);
_mm256_store_pd(c11,r11);_mm256_store_pd(c21,r21);_mm256_store_pd(c22,r22);
_mm256_store_pd((double*)pend2.m+pend2.pos+t,u4); /* park u4: this panel slot is dead now */
}
pend2.pos+=64;
if(pend2.pos>=1024*1024)pend2.active=false;
}
static void pend2_partial_finish(void) {
if(pending.active)pending_finish();
if(root_pending.active)root_finish();
while(pend2.active)pend2_partial_step();
/* reduced endgame: r12 = fl(fl(u4 + m3) + p12) */
for(int i=0;i<1024;i++){
double *c12=pend2.c+(size_t)i*pend2.ldc+1024;
const double *mm=pend2.m+(size_t)i*1024;
const double *pp=pend2.add+(size_t)i*2048+1024;
for(int j=0;j<1024;j+=16){
for(int t=0;t<16;t+=4){
__m256d u4=_mm256_load_pd(mm+j+t),m3=_mm256_load_pd(c12+j+t),p12=_mm256_load_pd(pp+j+t);
_mm256_store_pd(c12+j+t,_mm256_add_pd(p12,_mm256_add_pd(u4,m3)));
}
if(j+32<1024){
_mm_prefetch((const char*)(mm+j+32),_MM_HINT_T0);
_mm_prefetch((const char*)(c12+j+32),_MM_HINT_T0);
_mm_prefetch((const char*)(pp+j+32),_MM_HINT_T0);
}
if(i+1<1024)_mm_prefetch((const char*)(pend2.add+(size_t)(i+1)*2048+1024),_MM_HINT_T0);
}
}
pend2.active=false;
}
static const int lperm[7]={4,0,1,3,2,5,6};
static int g_fuse_sd[7], g_fuse_sgn[7];
static int g_fuse_sd_cur=0, g_fuse_sgn_cur=0;
static const int d1a2_bi[7]={0,0,1,2,3,0,2}, d1a2_bj[7]={3,-1,3,0,-1,1,3}, d1a2_bs[7]={1,0,-1,-1,0,1,1};
static void d1a2_fuse_tab(void){
for(int r=0;r<7;r++){g_fuse_sd[r]=0;g_fuse_sgn[r]=0;}
for(int r=1;r<7;r++){
int q=lperm[r], pq=lperm[r-1];
if(d1a2_bs[pq]!=0) continue;
if(d1a2_bj[q]==d1a2_bi[pq]){ g_fuse_sd[r]=1; g_fuse_sgn[r]=d1a2_bs[q]; }
else if(d1a2_bi[q]==d1a2_bi[pq]){ g_fuse_sd[r]=2; g_fuse_sgn[r]=d1a2_bs[q]; }
}
}
__attribute__((noinline)) static void pk_fuse(int n,Mat b,double *BP,int sgn,int side){
const double *sx=b.x; const size_t sl=b.ld;
for(int k=0;k<n;k+=2)for(int j=0;j<n;j+=8) {
double *p=BP+(size_t)j*(n+1)+k*8;
__m256d q0=_mm256_load_pd(p+0),q1=_mm256_load_pd(p+4);
__m256d q2=_mm256_load_pd(p+8),q3=_mm256_load_pd(p+12);
__m256d s0=_mm256_load_pd(sx+(size_t)(k+0)*sl+j), s1=_mm256_load_pd(sx+(size_t)(k+0)*sl+j+4);
__m256d s2=_mm256_load_pd(sx+(size_t)(k+1)*sl+j), s3=_mm256_load_pd(sx+(size_t)(k+1)*sl+j+4);
if(side==1){
if(sgn>0){q0=_mm256_add_pd(s0,q0);q1=_mm256_add_pd(s1,q1);q2=_mm256_add_pd(s2,q2);q3=_mm256_add_pd(s3,q3);}
else {q0=_mm256_sub_pd(s0,q0);q1=_mm256_sub_pd(s1,q1);q2=_mm256_sub_pd(s2,q2);q3=_mm256_sub_pd(s3,q3);}
} else {
if(sgn>0){q0=_mm256_add_pd(q0,s0);q1=_mm256_add_pd(q1,s1);q2=_mm256_add_pd(q2,s2);q3=_mm256_add_pd(q3,s3);}
else {q0=_mm256_sub_pd(q0,s0);q1=_mm256_sub_pd(q1,s1);q2=_mm256_sub_pd(q2,s2);q3=_mm256_sub_pd(q3,s3);}
}
_mm256_store_pd(p+0,q0);_mm256_store_pd(p+4,q1);
_mm256_store_pd(p+8,q2);_mm256_store_pd(p+12,q3);
}
}
static void d1a2_pk(int n,Mat b,double *BP){
if(g_fuse_sd_cur) pk_fuse(n,b,BP,g_fuse_sgn_cur,g_fuse_sd_cur);
else for(int k=0;k<n;k+=2)for(int j=0;j<n;j+=8) {
double *p=BP+(size_t)j*(n+1)+k*8;
_mm256_store_pd(p+0,ld(b,k+0,j));_mm256_store_pd(p+4,ld(b,k+0,j+4));
_mm256_store_pd(p+8,ld(b,k+1,j));_mm256_store_pd(p+12,ld(b,k+1,j+4));
}
}
/* ===== [d4ah_] pre 相传送带吸收（tile-major 版） ===== */
static void pre(int h,Mat a,double *s1,double *s2,double *s4,double *s3,int ld3,bool right);
struct D4ahPP { const double *s1; int l1; const double *s2; int l2; double *d; int l3; int sg; };
struct D4ahPT { int h; int i, j; int pass; D4ahPP p[4]; };
static D4ahPT d4ah_pt[8];
static int d4ah_ptn=0;
static int d4ah_on=1;
static long long d4ah_pchunks=0;
static const double *d4ai_L0,*d4ai_L1,*d4ai_L2; static double *d4ai_D0,*d4ai_D1;
static int d4ai_mode=0, d4ai_side=0;
/* ===== [d4aj_] root-belt duty limit: at most one root chunk per D4AJ_R leaf calls ===== */
#ifndef D4AJ_R
#define D4AJ_R 2
#endif
static int d4aj_skip=0;
static void d4ah_push(int h,const D4ahPP *q){ /* 4 个 pass = 1 个任务（tile-major） */
if(d4ah_ptn>=8) return;
D4ahPT *t=&d4ah_pt[d4ah_ptn++];
t->h=h; t->i=0; t->j=0; t->pass=0;
for(int k=0;k<4;k++) t->p[k]=q[k];
}
static void d4ah_run(D4ahPT *t){ /* 从当前 (i,j,pass) 同步做完 */
for(int i=t->i;i<t->h;i++){
for(int ps=(i==t->i?t->pass:0); ps<4; ps++){
D4ahPP *p=&t->p[ps];
const double *r1=p->s1+(size_t)i*p->l1, *r2=p->s2+(size_t)i*p->l2;
double *rr=p->d+(size_t)i*p->l3;
for(int j=0;j<t->h;j+=4){
__m256d u=_mm256_loadu_pd(r1+j), v=_mm256_loadu_pd(r2+j);
_mm256_storeu_pd(rr+j, p->sg>0?_mm256_add_pd(u,v):_mm256_sub_pd(u,v));
}
}
}
t->i=t->h;
}
static int d4ah_pre_arm(int h,Mat a,double *s1,double *s2,double *s4,double *s3,int ld3,bool right){
int n0=d4ah_ptn;
if(!d4ah_on||a.sign||h<256||(h&255)) { pre(h,a,s1,s2,s4,s3,ld3,right); return n0; }
const double *x=a.x; int ld=a.ld;
D4ahPP q[4];
if(right){ /* p=b-a ; q=d-p ; r=q-c ; t=d-b */
q[0]=(D4ahPP){x+h,ld, x,ld, s1,h, -1};
q[1]=(D4ahPP){x+(size_t)h*ld+h,ld, s1,h, s2,h, -1};
q[2]=(D4ahPP){s2,h, x+(size_t)h*ld,ld, s4,h, -1};
q[3]=(D4ahPP){x+(size_t)h*ld+h,ld, x+h,ld, s3,ld3, -1};
} else { /* p=c+d ; q=p-a ; r=b-q ; t=a-c */
q[0]=(D4ahPP){x+(size_t)h*ld,ld, x+(size_t)h*ld+h,ld, s1,h, +1};
q[1]=(D4ahPP){s1,h, x,ld, s2,h, -1};
q[2]=(D4ahPP){x+h,ld, s2,h, s4,h, -1};
q[3]=(D4ahPP){x,ld, x+(size_t)h*ld,ld, s3,ld3, -1};
}
d4ah_push(h,q);
return n0;
}
static inline int d4ah_pre_chunk(const double **q1,const double **q2,double **qd,int *qsg){
if(!d4ah_ptn) return 0;
D4ahPT *t=&d4ah_pt[d4ah_ptn-1];
if(t->h>=2048){ /* root-level task: duty-limited (leaf falls back to kernel_plain) */
if(d4aj_skip>0){ d4aj_skip--; return 0; }
d4aj_skip=D4AJ_R-1;
}
D4ahPP *p=t->p;
size_t b0=(size_t)t->i*p[0].l3+t->j, b1=(size_t)t->i*p[1].l3+t->j, b2=(size_t)t->i*p[2].l3+t->j, b3=(size_t)t->i*p[3].l3+t->j;
int side=(p[0].sg>0)?0:1;
if(t->pass==0){
d4ai_L0=p[0].s1+(size_t)t->i*p[0].l1+t->j;
d4ai_L1=p[0].s2+(size_t)t->i*p[0].l2+t->j;
d4ai_L2=(side==0? p[1].s2+(size_t)t->i*p[1].l2 : p[1].s1+(size_t)t->i*p[1].l1)+t->j;
d4ai_D0=p[0].d+b0; d4ai_D1=p[3].d+b3;
} else {
if(side==0){ d4ai_L0=p[0].d+b0; d4ai_L1=p[1].s2+(size_t)t->i*p[1].l2+t->j; }
else { d4ai_L0=p[1].s1+(size_t)t->i*p[1].l1+t->j; d4ai_L1=p[0].d+b0; }
d4ai_L2=(side==0? p[2].s1+(size_t)t->i*p[2].l1 : p[2].s2+(size_t)t->i*p[2].l2)+t->j;
d4ai_D0=p[1].d+b1; d4ai_D1=p[2].d+b2;
}
d4ai_mode=t->pass; d4ai_side=side;
*q1=d4ai_L0; *q2=d4ai_L1; *qd=(double*)d4ai_L2; *qsg=side;
d4ah_pchunks++;
if(++t->pass>=2){ t->pass=0; t->j+=256; if(t->j>=t->h){ t->j=0; if(++t->i>=t->h) d4ah_ptn--; } }
return 1;
}
static void d4ah_pre_fin_to(int k){
while(d4ah_ptn>k){ D4ahPT *t=&d4ah_pt[d4ah_ptn-1]; d4ah_run(t); d4ah_ptn--; }
_mm_sfence(); /* NT stores of the fused belt must be visible before consumers read the panels */
}
__attribute__((noinline)) static void d4bc_carry(int n,const double *a,const double *b,double *c,int ldc,int mr,
const double *pks,const double *pky,double *pkd,int sg){
if(sg>0) kernel< 1>(n,a,b,c,ldc,mr,pks,pky,pkd);
else if(sg<0) kernel<-1>(n,a,b,c,ldc,mr,pks,pky,pkd);
else kernel< 2>(n,a,b,c,ldc,mr,pks,pky,pkd);
}
template<int TASKS> static void leaf_impl(int leaf_n,Mat a,Mat b,double *c,int ldc,const Mat *anext=nullptr,bool odd=false,bool needp0=true) {
constexpr int n=BASE;
d1a2_pk(n,b,bp);
if constexpr(TASKS&1)pending_prefetch();
double *pkb[2]={ap,ap2};
int cur=(odd?1:0);
if(needp0) { /* pack the first ic-block the ordinary way */
int mc0=n<MC?n:MC;
for(int r=0;r<mc0;r+=MR) {
int rows=mc0-r<MR?mc0-r:MR;
for(int i=0;i<rows;i++)for(int k=0;k<n;k+=4)
_mm256_store_pd(pkb[cur]+(size_t)(r+i)*(BASE+24)+k,ld(a,r+i,k));
for(int i=rows;i<MR;i++)for(int k=0;k<n;k+=4)
_mm256_store_pd(pkb[cur]+(size_t)(r+i)*(BASE+24)+k,_mm256_setzero_pd());
}
}
for(int ic=0;ic<n;) {
int mc=n-ic<=MC+MR?n-ic:MC;
if(mc==MC+2) {
/* d4bc_: carry the NEXT leaf's first ic-block A panel (rows 0..29) out of these calls */
const int d4bc_nch2=anext?30*2:0;
int d4bc_ck=0, d4bc_skip=(n/NR)*((MC+MR-1)/MR)-d4bc_nch2;
const int d4bc_qd=n/2;
for(int j=0;j<n;j+=16) {
for(int jj=j;jj<j+16;jj+=NR)for(int r=0;r<MC;r+=MR) {
if(d4bc_nch2 && d4bc_ck<d4bc_nch2 && --d4bc_skip<0) {
const int row=d4bc_ck/2, off=(d4bc_ck%2)*d4bc_qd;
const double *pks=anext->x+(size_t)row*anext->ld+(size_t)off;
const double *pky=anext->sign?(anext->y+(size_t)row*anext->ly+(size_t)off):pk_zero;
double *pkd=pkb[cur^1]+(size_t)row*(BASE+24)+(size_t)off;
d4bc_ck++;
d4bc_carry(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)jj*(n+1),c+(size_t)(ic+r)*ldc+jj,ldc,MR,pks,pky,pkd,anext->sign);
} else
kernel_plain(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)jj*(n+1),c+(size_t)(ic+r)*ldc+jj,ldc,MR);
if((TASKS&1) && pending.active && (++pending_tick%PERIOD)==0 && (!pending.prog || pending.pos<BASE*pending.rowlim)) {pending_step_unchecked();if(pending.prog)pending_step_unchecked();pending_prefetch();}
if((TASKS&2) && pend2.active) { if(pend2.part)pend2_partial_step(); else pend2_step(); pend2_prefetch();}
if((TASKS&4) && root_pending.active && !pending.active) {root_step();root_prefetch();}
}
kernel_tail16(n,pkb[cur]+(size_t)MC*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+MC)*ldc+j,ldc);
}
ic+=mc;
if((TASKS&1) && pending.prog)pending.rowlim=ic;
continue;
}
/* ---- the NEXT ic-block's A panel is carried by this block's kernel calls ---- */
const int ic2=ic+mc;
const int nxt=(ic2<n)?(n-ic2<=MC+MR?n-ic2:MC):0;
const int rows2=((nxt+MR-1)/MR)*MR; /* rows incl. the zero pad rows */
const int nch=rows2*2; /* 2 chunks per row */
const int qd=n/2; /* doubles per chunk */
int ck=0; int skip=(n/NR)*((mc+MR-1)/MR)-nch;
for(int j=0;j<n;j+=NR)for(int r=0;r<mc;r+=MR) {
if(mc-r==2) { /* md9_t2: 2-row tail body (was 6-row body + tmp copy) */
kernel2(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc);
continue;
}
const double *pks=pk_zero,*pky=pk_zero,*pkm=pk_mp;
double *pkd=pk_zero;
int pkc=0; int psg=0;
if(ck<nch && --skip<0) {
const int row=ck/2, off=(ck%2)*qd;
if(row<nxt) {
pks=a.x+(size_t)(ic2+row)*a.ld+(size_t)off;
pky=a.sign?(a.y+(size_t)(ic2+row)*a.ly+(size_t)off):pk_zero;
pkm=a.sign>0?pk_mp:pk_mn;
}
pkd=pkb[cur^1]+(size_t)row*(BASE+24)+(size_t)off;
pkc=8;
ck++;
}
if(!pkc && d4ah_ptn){
if(d4ah_pre_chunk(&pks,&pky,&pkd,&psg)) pkc=9;
}
if(pkc==9){
const double *PL0=d4ai_L0,*PL1=d4ai_L1,*PL2=d4ai_L2; double *PD0=d4ai_D0,*PD1=d4ai_D1;
if(d4ai_mode==0){ if(d4ai_side==0) kpA< 1>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,PL0,PL1,PL2,PD0,PD1);
else kpA<-1>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,PL0,PL1,PL2,PD0,PD1); }
else { if(d4ai_side==0) kpB< 0>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,PL0,PL1,PL2,PD0,PD1);
else kpB< 1>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,PL0,PL1,PL2,PD0,PD1); }
}
else if(!pkc) kernel_plain(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR);
else if(psg>0) kernel<1>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,pks,pky,pkd);
else if(psg<0) kernel<-1>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,pks,pky,pkd);
else if(a.sign==0) kernel<2>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,pks,pky,pkd);
else if(a.sign>0) kernel<1>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,pks,pky,pkd);
else kernel<-1>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,pks,pky,pkd);
if((TASKS&1) && pending.active && (++pending_tick%PERIOD)==0 && (!pending.prog || pending.pos<BASE*pending.rowlim)) {
pending_step_unchecked(); if(pending.prog)pending_step_unchecked(); pending_prefetch();
}
if((TASKS&2) && pend2.active) { if(pend2.part)pend2_partial_step(); else pend2_step(); pend2_prefetch(); }
if((TASKS&4) && root_pending.active && !pending.active) {root_step();root_prefetch();}
}
/* ---- finish whatever the kernel loop did not carry ---- */
for(int cc=ck;cc<nch;cc++) {
const int row=cc/2, off=(cc%2)*qd;
if(row<nxt) {
for(int k=off;k<off+qd;k+=4)
_mm256_store_pd(pkb[cur^1]+(size_t)row*(BASE+24)+k, ld(a,ic2+row,k));
} else {
for(int k=off;k<off+qd;k+=4)
_mm256_store_pd(pkb[cur^1]+(size_t)row*(BASE+24)+k, _mm256_setzero_pd());
}
}
cur^=1;
ic+=mc;
if((TASKS&1) && pending.prog)pending.rowlim=ic;
}
if(n>=NTMIN)_mm_sfence();
}
static void leaf(int n,Mat a,Mat b,double*c,int ldc,const Mat *anext=nullptr,bool odd=false,bool needp0=true){int mode=pending.active+2*pend2.active+4*root_pending.active;switch(mode){case 0:leaf_impl<0>(n,a,b,c,ldc,anext,odd,needp0);break;case 1:leaf_impl<1>(n,a,b,c,ldc,anext,odd,needp0);break;case 2:leaf_impl<2>(n,a,b,c,ldc,anext,odd,needp0);break;case 3:leaf_impl<3>(n,a,b,c,ldc,anext,odd,needp0);break;case 4:leaf_impl<4>(n,a,b,c,ldc,anext,odd,needp0);break;case 5:leaf_impl<5>(n,a,b,c,ldc,anext,odd,needp0);break;case 6:leaf_impl<6>(n,a,b,c,ldc,anext,odd,needp0);break;case 7:leaf_impl<7>(n,a,b,c,ldc,anext,odd,needp0);break;}}
static void sum(int h,double *dst,Mat a,Mat b,int sign) {
for(int i=0;i<h;i++)for(int j=0;j<h;j+=8) {
__m256d v=sign>0?_mm256_add_pd(ld(a,i,j),ld(b,i,j)):_mm256_sub_pd(ld(a,i,j),ld(b,i,j));
__m256d w=sign>0?_mm256_add_pd(ld(a,i,j+4),ld(b,i,j+4)):_mm256_sub_pd(ld(a,i,j+4),ld(b,i,j+4));
asm volatile("" : "+x"(v), "+x"(w) : : "memory");
if(h>=SUMNT) {
_mm256_stream_pd(dst+(size_t)i*h+j,v);_mm256_stream_pd(dst+(size_t)i*h+j+4,w);
} else {
_mm256_store_pd(dst+(size_t)i*h+j,v);_mm256_store_pd(dst+(size_t)i*h+j+4,w);
}
}
if(h>=SUMNT)_mm_sfence();
}
static void combine(int h,double *c,int ldc,const double *m,bool wino=false) {
size_t z=(size_t)h*h+(wino?UPPER_PAGE_PAD:0);
for(int i=0;i<h;i++)for(int j=0;j<h;j+=4) {
size_t t=(size_t)i*h+j;
double *c11=c+(size_t)i*ldc+j,*c12=c11+h,*c21=c+(size_t)(i+h)*ldc+j,*c22=c21+h;
__m256d m1=_mm256_load_pd(c11),m2=_mm256_load_pd(c21),m3=_mm256_load_pd(c12),m4=_mm256_load_pd(c22),
m5=_mm256_load_pd(m+t),m6=_mm256_load_pd(m+z+t),m7=_mm256_load_pd(m+2*z+t);
__m256d r11,r12,r21,r22;
if(wino) {
__m256d u2=_mm256_add_pd(m1,m6),u3=_mm256_add_pd(u2,m7),u4=_mm256_add_pd(u2,m5);
r11=_mm256_add_pd(m1,m2);r12=_mm256_add_pd(u4,m3);r21=_mm256_sub_pd(u3,m4);r22=_mm256_add_pd(u3,m5);
} else {
r11=_mm256_add_pd(_mm256_sub_pd(_mm256_add_pd(m1,m4),m5),m7);r12=_mm256_add_pd(m3,m5);r21=_mm256_add_pd(m2,m4);
r22=_mm256_add_pd(_mm256_add_pd(_mm256_sub_pd(m1,m2),m3),m6);
}
_mm256_store_pd(c11,r11);_mm256_store_pd(c12,r12);_mm256_store_pd(c21,r21);_mm256_store_pd(c22,r22);
}
}
static inline __m256d raw_pre_ld(Mat a,int i,int j){return _mm256_load_pd(a.x+(size_t)i*a.ld+j);}
static void pre(int h,Mat a,double *s1,double *s2,double *s4,double *s3,int ld3,bool right) {
for(int i=0;i<h;i++)for(int j=0;j<h;j+=8) {
__m256d a0=raw_pre_ld(a,i,j),a1=raw_pre_ld(a,i,j+4),b0=raw_pre_ld(a,i,j+h),b1=raw_pre_ld(a,i,j+h+4),
c0=raw_pre_ld(a,i+h,j),c1=raw_pre_ld(a,i+h,j+4),d0=raw_pre_ld(a,i+h,j+h),d1=raw_pre_ld(a,i+h,j+h+4);
__m256d p0,p1,q0,q1,r0,r1,t0,t1;
if(!right) {
p0=_mm256_add_pd(c0,d0);p1=_mm256_add_pd(c1,d1);
q0=_mm256_sub_pd(p0,a0);q1=_mm256_sub_pd(p1,a1);
r0=_mm256_sub_pd(b0,q0);r1=_mm256_sub_pd(b1,q1);
t0=_mm256_sub_pd(a0,c0);t1=_mm256_sub_pd(a1,c1);
}else {
p0=_mm256_sub_pd(b0,a0);p1=_mm256_sub_pd(b1,a1);
q0=_mm256_sub_pd(d0,p0);q1=_mm256_sub_pd(d1,p1);
r0=_mm256_sub_pd(q0,c0);r1=_mm256_sub_pd(q1,c1);
t0=_mm256_sub_pd(d0,b0);t1=_mm256_sub_pd(d1,b1);
}
size_t off=(size_t)i*h+j,off3=(size_t)i*ld3+j;
_mm256_stream_pd(s1+off,p0);_mm256_stream_pd(s1+off+4,p1);
_mm256_stream_pd(s2+off,q0);_mm256_stream_pd(s2+off+4,q1);
_mm256_stream_pd(s4+off,r0);_mm256_stream_pd(s4+off+4,r1);
_mm256_stream_pd(s3+off3,t0);_mm256_stream_pd(s3+off3+4,t1);
}
_mm_sfence();
}
template<int N> static void rec(int rec_n,Mat a,Mat b,double *c,int ldc,double *work,bool root_last=false) {
constexpr int n=N;
if constexpr(n<=BASE){leaf(n,a,b,c,ldc);return;} else {
constexpr int h=n/2;size_t z=(size_t)h*h;
Mat aa[4]={{a.x,a.ld},{a.x+h,a.ld},{a.x+(size_t)h*a.ld,a.ld},{a.x+(size_t)h*a.ld+h,a.ld}};
Mat bb[4]={{b.x,b.ld},{b.x+h,b.ld},{b.x+(size_t)h*b.ld,b.ld},{b.x+(size_t)h*b.ld+h,b.ld}};
double *m=n==2*BASE?work+LEAF_PAD:work;
if(n==2*BASE && pending.active && pending.m==m)m=pending_work+LEAF_PAD;
size_t pz=z+(n==2*BASE?LEAF_PAD:UPPER_PAGE_PAD);
double *p=work+3*z,*q=p+z,*sub=q+z;
double *out[7]={c,c+(size_t)h*ldc,c+h,c+(size_t)h*ldc+h,m,m+pz,m+2*pz};
if(n>BASE*2) {
pending_finish();
double *pa0=m+3*pz,*pa1=pa0+pz,*pa2=pa1+pz,*pb0=pa2+pz,*pb1=pb0+pz,*pb2=pb1+pz;
sub=pb2+pz;
int d4ah_n0=d4ah_ptn;
const bool d4bb_out=(ldc==2*h); const int d4bb_l3=d4bb_out? h : ldc;
d4ah_pre_arm(h,a,pa0,pa1,pa2,c,d4bb_l3,false); /* d4bb_ 先压 a 侧（栈底） */
d4ah_pre_arm(h,b,pb0,pb1,pb2,d4bb_out? c+(size_t)h*h : c+h,d4bb_l3,true); /* d4bb_ 再压 b 侧 */
/* mdA_dp2: child order [1,3,4,5,6,0,2] -- constraint "6 before 0 and 2" is the only
hard one; this order puts the three m-panel writers (4,5,6) as late as possible,
which IS the deferral window for the previous sibling's deferred h=1024 combine. */
rec<h>(h,aa[1],bb[2],out[1],ldc,sub);
d4ah_pre_fin_to(d4ah_n0+1);
rec<h>(h,aa[3],{pb2,h},out[3],ldc,sub);
if(root_last) {
pend2_finish(); // M1 of the root is now fully reconstructed.
root_pending.active=true;
root_prefetch();
}
#if !B2M_NO_BARA
pend2_finish(); /* barrier: next call writes m, the deferred products */
#endif
d4ah_pre_fin_to(d4ah_n0);
rec<h>(h,{pa0,h},{pb0,h},out[4],h,sub);
rec<h>(h,{pa1,h},{pb1,h},out[5],h,sub);
if(d4bb_out) rec<h>(h,{c,h},{c+(size_t)h*h,h},out[6],h,sub); /* d4bb_ s3 面板连续存放 */
else rec<h>(h,{c,ldc},{c+h,ldc},out[6],h,sub);
rec<h>(h,aa[0],bb[0],out[0],ldc,sub);
if(n==4096) {
root_pending={c,ldc,m,0,false};
rec<h>(h,{pa2,h},bb[3],out[2],ldc,sub,true);
pending_finish();root_finish();pend2_finish();
return;
}
/* d4ak_: for the LAST node the deferred combine is split -- the six M3-free
outputs (+ their root-u4 ADD terms) are produced during the last child's
leaves; only r12 = fl(fl(u4+m3)+p12) is left for the endgame. */
if(root_last){ pend2={c,ldc,m,h,(int)pz,0,true}; pend2.add=root_pending.m; pend2.part=true; }
rec<h>(h,{pa2,h},bb[3],out[2],ldc,sub);
pending_finish();
if(root_last){ root_finish(); pend2_partial_finish(); }
else {
pend2_finish(); /* finish an older deferred one */
pend2={c,ldc,m,h,(int)pz,0,true}; /* defer THIS combine (wino form);
m-plane stride = the engine's pz */
}
return;
}
const int ai[7]={0,2,0,3,0,2,1},aj[7]={3,3,-1,-1,1,0,3},as[7]={1,1,0,0,1,-1,-1};
const int bi[7]={0,0,1,2,3,0,2},bj[7]={3,-1,3,0,-1,1,3},bs[7]={1,0,-1,-1,0,1,1};
if(h<=BASE) {
/* Independent product order; M7 is last so reconstruction may
advance behind the produced-row frontier of that leaf. */
Mat xa[7],ya[7];
for(int rr=0;rr<7;rr++){
const int qq=lperm[rr];
xa[rr]=aa[ai[qq]];ya[rr]=bb[bi[qq]];
if(as[qq]){xa[rr].y=aa[aj[qq]].x;xa[rr].ly=a.ld;xa[rr].sign=as[qq];}
if(bs[qq]){ya[rr].y=bb[bj[qq]].x;ya[rr].ly=b.ld;ya[rr].sign=bs[qq];}
}
for(int rr=0;rr<7;rr++){
const int qq=lperm[rr];
g_fuse_sd_cur=g_fuse_sd[rr];g_fuse_sgn_cur=g_fuse_sgn[rr];
Mat x=xa[rr],y=ya[rr];
const Mat *nx=(rr<6)?&xa[rr+1]:nullptr;
const bool oddp=((rr&1)!=0), np0=(rr==0);
if(rr==6) {
pending_finish();
pending={c,ldc,m,0,true,0,true};
leaf(h,x,y,out[qq],qq<4?ldc:h,nx,oddp,np0);
pending.prog=false;
} else leaf(h,x,y,out[qq],qq<4?ldc:h,nx,oddp,np0);
}
} else {
for(int r=0;r<7;r++) {
Mat x=aa[ai[r]],y=bb[bi[r]];
if(as[r]){sum(h,p,x,aa[aj[r]],as[r]);x={p,h};}
if(bs[r]){sum(h,q,y,bb[bj[r]],bs[r]);y={q,h};}
rec<h>(h,x,y,out[r],r<4?ldc:h,sub);
}
}
if(h==BASE) {
/* Keep any final portion deferred into the next sibling. */
} else {
pending_finish();
combine(h,c,ldc,m);
}
}
}
void matrix_multiply(int n,const double *a,const double *b,double *c) {
d1a2_fuse_tab();
madvise((void*)workspace,sizeof(workspace),MADV_HUGEPAGE);
madvise((void*)pending_work,sizeof(pending_work),MADV_HUGEPAGE);
rec<4096>(n,{a,n},{b,n},c,n,workspace);
d4ah_pre_fin_to(0);
pending_finish();
pend2_finish();
_mm_sfence();
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 1.859 s | 502 MB + 356 KB | Accepted | Score: 100 | 显示更多 |