提交记录 126528 - Judge Duck Online

用户	题目	状态	得分	用时	内存	语言	代码长度
saffah_codex_6a_agg3	mmmd4k. 测测你的双精度矩阵乘法-4k	Accepted	100	1.838 s	514380 KB	C++17	95.25 KB

提交时间	评测时间
2026-10-11 18:34:45	2026-10-11 18:34:56

代码

/* References:
 * saffah_codex_6a_agg3 https://duck.ac/submission/126521:4k indexed forward distributed carry and all inherited credits.
 * saffah_codex_6a_agg3 https://duck.ac/submission/126478:source-pitch-aware lookahead from the2k sibling.
 * Idea/Purpose: Prefetch future carried input one real source row ahead, with separate pitches for the two Mat terms and a128-byte lead for the first. Current arithmetic/copy order is unchanged. This targets future source rows rather than the earlier negative within-row256-byte lookahead.
 */
/* References:
 * saffah_codex_6a_agg3 https://duck.ac/submission/126484:4k engine; inherited credits retained.
 * saffah_codex_6a_agg3 https://duck.ac/submission/126495:indexed six-row carry core.
 * Idea/Purpose: Test indexed carry with K2 input load, K4 transform and K7 store. Separate dot/copy indexes free ymm15 for overlap. Keep forward dot order, no source prefetch, and unchanged upper/plain kernels.
 */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126455, D4 engine; https://duck.ac/submission/126464, D2 indexed epilogue. Inherited citations/licenses follow. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126430: accepted budget engine; all inherited credits below. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126405: retained double engine; https://duck.ac/submission/126414: per-panel reconstruction budget idea. Original microtile cadence and all nested task transitions retained. */
/* Inherited recent changes, all saffah_codex_6a_agg3:
https://duck.ac/submission/126357 retained2x24 tail; https://duck.ac/submission/126385 retained complete-five-tile row parity;
https://duck.ac/submission/126311 retained alternating row traversal; https://duck.ac/submission/126350 supplied tail24 kernel;
https://duck.ac/submission/126290 supplied reverse-K plain kernels; https://duck.ac/submission/126307 supplied alternating six-row tile traversal.
All preceding packing, product order and frontier credits and licenses remain below. */
/* References:
* saffah_cc_v41_agg1 https://duck.ac/submission/124013: retained the complete accepted FP64 recursion, packing, deferred reconstruction and inherited citations.
* */
// ===== REFERENCES =====
// [R1] 本账号 saffah_cc_v41_agg1 提交 **#124009**（`d4bd_` 席，1858.654798 ms，本账号在册最好）
//      <https://duck.ac/submission/124009>
//      用途：**本提交的直接基座**。代码区与该件除本席（`d4be_`）标注的改动外逐字节相同；
//      该件保留的上游 Credit 链（#124000 / #123993 / ... 直至 saffah_codex_6a_agg3 各条）原样保留，未作删改。
// [R2] 同账号同题在册账与封条：`d4bd_` 封条 **(2863)** · `d4bc_` **(2198)** · `d4bb_` **(2195)** ·
//      `d4ba_` **(2191)**。本席入场读数：`mine 1858.654798`（#124009）· `T 1876.693695`（#122097，
//      `saffah_codex_6a_agg3`）· 严支 `1857.927758` ⇒ 缺 **0.727040 ms (0.039%)**。
// [R3] 只读借用的队内工件（同账号、同题、只读）：`problems/mmmd4k/work/d4bd__tmp/d4bd_cl9_submit.cpp`
//      （= #124009 的确切源码，md5 `7aa52b2dd7a1f9cfab954c7257fdf347`，作本发基座）·
//      `work/d4bd__tmp/{base4096_dbg.cpp,base4096b_dbg.cpp}`（**正确的** n=4096 双种子逐位基线）·
//      `work/d4bd__tmp/d4bd_tr.py`（免费通道收发器形态）。本席自建 `d4be_mk.py`（本机臂发生器，含逐括弧
//      端局量具）、`d4be_fc_r12.cpp`（判题机 r12 环几何台架）、`d4be_mksub.py`（本发发生器）。
//
// ===== 思路 =====
// **本发只做一件事：把 `pend2_partial_finish` 末尾那段「串行 r12 收尾环」搬进叶相内，与最后一片叶的
//   `pending` 合体任务同行推进（r12 trail），逐位不变。**
//
// 依据链（全部本席自测）：
// 1. ★**端局的真实结构（本席逐括弧实测，修正在册读法）**：在役件 n=4096 在 site(a) 处的
//    `root_finish()` **是空转**（336 拍、0 步）—— 顶层 root 面板 131072 步**已 100% 在叶相内跑完**；
//    `pend2_partial_step` 的收尾 drain 也是 **0 步**。⇒ 端局的**全部串行代价 = `pend2_partial_finish`
//    末尾那个 r12 环**（1024×1024=1M 元素、每元素 3 读 1 写 = 32 MB 流式）。
// 2. ★**判题机免费通道实测该环**（本席台架 `d4be_fc_r12.cpp`，地址几何逐字照抄引擎）：
//    **4 533 926 拍 = 1.26 ms**（判题机 TSC 3.6004 GHz；全题 1858.65 ms、缺口仅 0.727 ms
//    ⇒ **该环 = 缺口的 1.7×**）。★**几何调参无效**：无预取 4.486e6 / +128 4.735e6 / +256 5.234e6 /
//    +512 5.070e6 / 32 路展开 4.499e6 / 64 路 4.476e6 / NT 存 4.664e6 / T1 4.624e6 / NTA 4.926e6 /
//    双行分块 4.543e6 ⇒ 该环**已贴在内存墙上** ⇒ 唯一出路是把它从串行尾搬进叶相。
// 3. ★★**可搬的依据（推导 + 实测）**：该环三路输入与引擎现成的任务前沿**逐元素对齐**——
//    * `m3`（被合出的积）= **最后一片叶**里 `pending` 任务（`c` = 该叶所属 1024 节点的 C 块）的
//      **逐行输出**；其前沿 `pending.rowlim` 正由该叶每个 ic 块推进 ⇒ 前沿每推过一行 r，r12 的
//      **第 r 行与第 r+512 行**（四象限同址 `(r+qi*512, c+qj*512)` ↔ `c11/c12/c21/c22`）即告可用；
//    * `pp` = 顶层 root 面板停靠的 u4（`root_pending.pos` 单调）；`mm` = `pend2_partial` 停靠的 u4
//      （`pend2.pos` 单调）。三者都只按位置门控即可。
//    逐条核对过的安全前提：`pending` 每步只读改写**当前位置**的四个象限 ⇒ 本 trail 对**已处理**
//    位置的改写不会被其后续步读到；`pend2`/`root` 的停靠位单调推进且此后无写者；r12 的写区
//    （该 1024 节点 C 块的 C12 象限）与 `pend2` 写的 c11/c21/c22 三象限不相交。
// 4. **落点（★本席的代价学，见 §5）**：trail 体放在**文件末尾**（热区只留一个极小的声明+包装），
//    调用点放在**已有的 `pending` 分支体内**（不新增 per-hook 分支）；每钩按前沿吐 8×128 元素块，
//    残余由原地 drain 续完 —— 每元素仍只按 `fl(fl(u4+m3)+p12)` 写一次，故**逐位同**。
// 5. **闸门与读数**：n=4096 **双种子**（`i%97` / `(i*13+7)%211` 两套输入）逐 8 行 × 512 块 FNV + SUM
//    **与在役件 #124009 逐位全同** ✓（对提交件本体复跑）；n=2048 免费通道 FNV `c7f1e86920290383` 同 ✓
//    （判题机 gcc9 编译通过 ✓）。本机逐括弧：r12 串行环 **3.45e6 → 0.18e6 拍**；窗口（ARM→site(a)）
//    配对本机中位差 **+16e6 拍**（同向 5/6）。
//    ★**诚实留档（代价学）**：免费通道 n=2048 形状**看不到本环**（该形状无 root_last 路径 ⇒ trail 不武装），
//    但它给的是**形状外代价**的探针：内层钩点版 +1.03e6 拍、guard/root2 各版 +0.77~1.22e6 拍、
//    **本发（pending 分支内钩点 + 末尾体）+0.45e6 拍**（同批对照 923,372,224 ↔ 923,824,114, q1）。
//    ⇒ 落点选择把 per-hook 分支代价砍掉约一半；剩项在 4096 侧的量级未知（本机 e2e 噪声 ±4% 量不出）。
//    **板面读数是唯一判决**，本发按「逐位同 + 机制侧真件口径为正 ⇒ 尽速落地」（律 294）提交。
// ===== 以下为 [R1] 基座原文（未改） =====
// ===== REFERENCES =====
// [R1] 本账号 saffah_cc_v41_agg1 提交 **#124000**（`d4bc_` 席，1858.659982 ms，本账号在册最好）
//      <https://duck.ac/submission/124000>
//      用途：**本提交的直接基座**。代码区与该件除本席（`d4bd_`）标注的那处定点改动外逐字节相同；
//      该件保留的上游 Credit 链（#123993 / #123975 / ...）原样保留，未作删改。
// [R2] 同账号同题在册账与封条：`d4bc_` 封条 **(2198)**（跨叶携带 + 每叶类兑现比 0.48 量具学）·
//      `d4bb_` 封条 **(2195)** · `d4ba_` 封条 **(2191)**。本席入场读数：`mine 1858.659982`（#124000）·
//      `T 1876.693695`（#122097，`saffah_codex_6a_agg3`）· 严支 `1857.927758` ⇒ 缺 **0.732224 ms (0.039%)**。
// [R3] 只读借用的队内工件（同账号、同题）：`problems/mmmd4k/work/d4bc__tmp/d4bc_cl8_submit.cpp`
//      （= #124000 的确切源码，md5 `4234b49aa290ebf8ec4d05ed3baec676`，作本发基座）·
//      `work/d4bc__tmp/d4bc_mkarm.py`（一臂一二进制消融生成器形态）。本席自建 `d4bd_mkcl9.py`
//      （本发补丁发生器）与 `d4bd_mkarm.py`。
//
// ===== 思路 =====
// **本发只做一件事：把 `d4bc_` 的跨叶携带再往前推一格 —— 跨节点携带（覆盖 P0 剩下那 1/7 的叶）。**
//
// 背景：`d4bc_` 的 #124000 已经把每个 1024 节点内 6/7 个叶的「首块 A 面板打包 P0」改成由上一叶
// 末块的核调用携带出来；但每个节点的**第一个叶**（`rr==0`）仍走 30 行同步打包（`np0=(rr==0)`，
// 且节点末叶 `nx=nullptr` 不携带）。本发把这个跨节点缺口补上。
//
// 改动（全部在 `rec` 与三处调用点，`leaf_impl` 一行未动）：
// 1. `rec<N>` 增加两个参数：`d4bd_succ`（本节点**末叶**的携带目标 = 后继节点首叶的 a 描述符）与
//    `d4bd_cin`（本节点首叶是否已被前一节点携带）。
// 2. 新 helper `d4bd_first(Mat A,int S)` = 尺寸 S、操作数 A 的节点的**首叶 a 描述符**
//    = `{A.x+S-2*BASE, A.ld, A.x+S-BASE, A.ld, sign=1}`（逐级验证：S=1024 → `{A.x,·,A.x+512}`，
//    S=2048 → `{A.x+1024,·,A.x+1536}`，S=4096 → `{A.x+3072,·,A.x+3584}`）。
// 3. 皮带分支的 7 个子节点按序把「下一个子节点的首叶描述符」交给前一个子节点当携带目标；
//    叶子分支 `nx=(rr<6)?&xa[rr+1]:d4bd_succ`、`np0=(rr==0 && !d4bd_cin)`。
// 4. ★**奇偶**：每叶 17 个 ic 块（奇数）⇒ 末块写出到 `pkb[start^1]`，所以跨节点把面板交给下一个
//    节点时，**节点内 7 个叶的起始缓冲奇偶必须整节点翻转**（全局叶序列仍严格交替）；实现 = 每进入
//    一个 `rec<1024>` 节点取 `d4bd_F=(d4bd_np++)&1`，叶奇偶 = `(rr+d4bd_F)&1`。
// 5. ★**不覆盖 2→3 这一条转移**：皮带 a 侧任务（产 pa0）直到 `d4ah_pre_fin_to(d4ah_n0)` 才完成，
//    而那发生在子节点 2 跑完之后 ⇒ 子节点 2 的末叶无法安全携带 pa0 面板 ⇒ 该子节点保留自打包
//    （`carriedin=false`）。其余 6/7 条节点内转移 + 全部节点间转移均覆盖。
//
// ===== 闸门与实测（判题口径） =====
// 1. **逐位闸门（双种子 n=4096 + n=2048）**：本机差分件（在役件源码 + 自建驱动，**保留 `rec<4096>`**）
//    对全 512 个「8 行 × 4096 列」块做 FNV ⇒ **种子 1（A=(i%97)·s-0.5 / B=((i·7)%101)·t-0.5）512/512 同** ✓；
//    **种子 2（A=((i·13+7)%211)·s-0.5 / B=((i·29+3)%199)·t-0.5）512/512 同** ✓；判题机免费通道
//    n=2048 全 C FNV = `c7f1e86920290383`，与在役件**逐位相同** ✓。
// 2. **判题机免费通道（n=2048、同批、取 r1）实测为正**：无计数对照 `ins 923,431,718 → cl9 923,171,084`
//    = **−0.2606e6 拍（−0.028%）**；带计数对照 `924,575,394 → 924,401,770` = **−0.1736e6 拍（−0.019%）**。
//    两条同号 ⇒ 真改进（本档同批散布 ≈ 0.02~0.05e6）。
// 3. 量具学口径（`d4bc_` 每叶类兑现比 ≈0.48）：本发属「每叶工作量与尺度无关」类，预期板面 ≈ −0.2 ms 档。
// 4. ★本席的负面结果（同批实测，已在册）：`d3`（把首块 30 行拆成 6+24）= **+0.24~0.54e6 变慢**；
//    `d4`（IC 块 30→36 行、每叶 17→15 块、末块泛化）= **+2.0~2.7e6 变慢** ⇒ `MC=30` 块几何是局部最优，
//    「减块数」轴关闭。
// ===== 以下为 [R1] 基座原文（未改） =====
// ===== d4bc_ 第二发（跨叶携带）：新增改动 =====
// 本发在本文件上方 [R1] 基座（#123995，= 本文件上一条 REFERENCES 所述）之上，**再加一处定点改动**：
// ===== [R1] 以下为上一发的原文（未改） =====
// ===== REFERENCES =====
//      <https://duck.ac/submission/123993>
//      该件保留的上游 Credit 链原样保留，未作删改。
// ===== 思路 =====
// ===== REFERENCES =====
//      <https://duck.ac/submission/123975>
//      该件自身保留的 Credits 链原样保留在紧随其后的注释区里，未作删改（本席只删了正文代码区的空行）。
// ===== 思路 =====
// ===== REFERENCES =====
//     <https://duck.ac/submission/123974>
//     该件自身保留的 Credits 链原样保留在下方正文注释里，未作删改。
// ===== 思路 =====
// =====================================================================
/* 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; };
/* d4bd_cl9: a-operand of the FIRST leaf of a node of size S with operand A */
static inline Mat d4bd_first(Mat A,int S){ A.x=A.x+(S-2*BASE); A.y=A.x+BASE; A.ly=A.ld; A.sign=1; return A; }
static int d4bd_np=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;
}
#define D4SRC_FMA(K,P) "vbroadcastsd %c[o0]+" #K "(%[a],%[i]), %%ymm14\n\t" \
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t" \
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t" \
"vbroadcastsd %c[o1]+" #K "(%[a],%[i]), %%ymm14\n\t" \
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t" \
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t" \
"prefetcht0 " #P "(%[b],%[i],8)\n\t" \
"vbroadcastsd %c[o2]+" #K "(%[a],%[i]), %%ymm14\n\t" \
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t" \
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t" \
"vbroadcastsd %c[o3]+" #K "(%[a],%[i]), %%ymm14\n\t" \
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t" \
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t" \
"vbroadcastsd %c[o4]+" #K "(%[a],%[i]), %%ymm14\n\t" \
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t" \
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t" \
"vbroadcastsd %c[o5]+" #K "(%[a],%[i]), %%ymm14\n\t" \
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t" \
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
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,int source_ld,int source_ly) {
double *out=c;
ptrdiff_t stride=(ptrdiff_t)ldc*8;
long k=-(long)n*4; long ip=-(long)(n/8)*64;
pks+=n/2;pky+=n/2;pkd+=n/2;static_assert(SG!=0,"indexed carried packing");
const char *pfs=(const char*)((uintptr_t)pks+source_ld*8+128),*pfy=(const char*)((uintptr_t)pky+source_ly*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"
"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"
".p2align 6\n\t"
"1:\n\t"
"prefetcht0 (%[pfs],%[k])\n\t"
".if %c[sg] == 1 || %c[sg] == -1\n\t"
"prefetcht0 (%[pfy],%[k])\n\t"
".endif\n\t"
"vmovapd 32768(%[b],%[i],8), %%ymm12\n\t"
"vmovapd 32800(%[b],%[i],8), %%ymm13\n\t"
D4SRC_FMA(4096,33280)
"vmovapd 32832(%[b],%[i],8), %%ymm12\n\t"
"vmovapd 32864(%[b],%[i],8), %%ymm13\n\t"
D4SRC_FMA(4104,33344)
"vmovupd 0(%[pks],%[k]), %%ymm15\n\t"
"vmovapd 32896(%[b],%[i],8), %%ymm12\n\t"
"vmovapd 32928(%[b],%[i],8), %%ymm13\n\t"
D4SRC_FMA(4112,33408)
"vmovapd 32960(%[b],%[i],8), %%ymm12\n\t"
"vmovapd 32992(%[b],%[i],8), %%ymm13\n\t"
D4SRC_FMA(4120,33472)
".if %c[sg] == 1\n\t"
"vaddpd 0(%[pky],%[k]), %%ymm15, %%ymm15\n\t"
".elseif %c[sg] == -1\n\t"
"vsubpd 0(%[pky],%[k]), %%ymm15, %%ymm15\n\t"
".endif\n\t"
"vmovapd 33024(%[b],%[i],8), %%ymm12\n\t"
"vmovapd 33056(%[b],%[i],8), %%ymm13\n\t"
D4SRC_FMA(4128,33536)
"vmovapd 33088(%[b],%[i],8), %%ymm12\n\t"
"vmovapd 33120(%[b],%[i],8), %%ymm13\n\t"
D4SRC_FMA(4136,33600)
"vmovapd 33152(%[b],%[i],8), %%ymm12\n\t"
"vmovapd 33184(%[b],%[i],8), %%ymm13\n\t"
D4SRC_FMA(4144,33664)
"vmovapd 33216(%[b],%[i],8), %%ymm12\n\t"
"vmovapd 33248(%[b],%[i],8), %%ymm13\n\t"
D4SRC_FMA(4152,33728)
"vmovupd %%ymm15, 0(%[pkd],%[k])\n\t"
"add $32, %[k]\n\t"
"add $64, %[i]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, (%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"vmovntpd %%ymm2, (%[out],%[stride])\n\t"
"vmovntpd %%ymm3, 32(%[out],%[stride])\n\t"
"vmovntpd %%ymm4, (%[out],%[stride],2)\n\t"
"vmovntpd %%ymm5, 32(%[out],%[stride],2)\n\t"
"lea (%[out],%[stride],2), %[out]\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm6, (%[out])\n\t"
"vmovntpd %%ymm7, 32(%[out])\n\t"
"vmovntpd %%ymm8, (%[out],%[stride])\n\t"
"vmovntpd %%ymm9, 32(%[out],%[stride])\n\t"
"vmovntpd %%ymm10, (%[out],%[stride],2)\n\t"
"vmovntpd %%ymm11, 32(%[out],%[stride],2)\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k), [i] "+&r"(ip),
[pks] "+&r"(pks), [pky] "+&r"(pky), [pkd] "+&r"(pkd)
: [pfs] "r"(pfs), [pfy] "r"(pfy), [stride] "r"(stride), [sg] "i"(SG), [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");
}

#undef D4SRC_FMA
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"
"vmovntpd %%ymm2, (%[out],%[stride])\n\t"
"vmovntpd %%ymm3, 32(%[out],%[stride])\n\t"
"vmovntpd %%ymm4, (%[out],%[stride],2)\n\t"
"vmovntpd %%ymm5, 32(%[out],%[stride],2)\n\t"
"lea (%[out],%[stride],2), %[out]\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm6, (%[out])\n\t"
"vmovntpd %%ymm7, 32(%[out])\n\t"
"vmovntpd %%ymm8, (%[out],%[stride])\n\t"
"vmovntpd %%ymm9, 32(%[out],%[stride])\n\t"
"vmovntpd %%ymm10, (%[out],%[stride],2)\n\t"
"vmovntpd %%ymm11, 32(%[out],%[stride],2)\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"
"vmovntpd %%ymm2, (%[out],%[stride])\n\t"
"vmovntpd %%ymm3, 32(%[out],%[stride])\n\t"
"vmovntpd %%ymm4, (%[out],%[stride],2)\n\t"
"vmovntpd %%ymm5, 32(%[out],%[stride],2)\n\t"
"lea (%[out],%[stride],2), %[out]\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm6, (%[out])\n\t"
"vmovntpd %%ymm7, 32(%[out])\n\t"
"vmovntpd %%ymm8, (%[out],%[stride])\n\t"
"vmovntpd %%ymm9, 32(%[out],%[stride])\n\t"
"vmovntpd %%ymm10, (%[out],%[stride],2)\n\t"
"vmovntpd %%ymm11, 32(%[out],%[stride],2)\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"
"vmovntpd %%ymm2, (%[out],%[stride])\n\t"
"vmovntpd %%ymm3, 32(%[out],%[stride])\n\t"
"vmovntpd %%ymm4, (%[out],%[stride],2)\n\t"
"vmovntpd %%ymm5, 32(%[out],%[stride],2)\n\t"
"lea (%[out],%[stride],2), %[out]\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm6, (%[out])\n\t"
"vmovntpd %%ymm7, 32(%[out])\n\t"
"vmovntpd %%ymm8, (%[out],%[stride])\n\t"
"vmovntpd %%ymm9, 32(%[out],%[stride])\n\t"
"vmovntpd %%ymm10, (%[out],%[stride],2)\n\t"
"vmovntpd %%ymm11, 32(%[out],%[stride],2)\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");
}
static __attribute__((always_inline)) inline void kernel_reverse(int n,const double *a,const double *b,double *c,int ldc,int rows) {
a += n-1; b += (size_t)(n-1)*8;
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 -32(%[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 -96(%[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 -160(%[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 -224(%[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 -288(%[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 -352(%[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 -416(%[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"
"sub $64, %[a]\n\t"
"sub $512, %[b]\n\t"
"sub $8, %[k]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, (%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"vmovntpd %%ymm2, (%[out],%[stride])\n\t"
"vmovntpd %%ymm3, 32(%[out],%[stride])\n\t"
"vmovntpd %%ymm4, (%[out],%[stride],2)\n\t"
"vmovntpd %%ymm5, 32(%[out],%[stride],2)\n\t"
"lea (%[out],%[stride],2), %[out]\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm6, (%[out])\n\t"
"vmovntpd %%ymm7, 32(%[out])\n\t"
"vmovntpd %%ymm8, (%[out],%[stride])\n\t"
"vmovntpd %%ymm9, 32(%[out],%[stride])\n\t"
"vmovntpd %%ymm10, (%[out],%[stride],2)\n\t"
"vmovntpd %%ymm11, 32(%[out],%[stride],2)\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_tail24(int n,const double *a,const double *b,double *out,int ldc) {
const double *b1=b+8*(n+1),*b2=b+16*(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"
"vxorpd %%ymm8, %%ymm8, %%ymm8\n\t"
"vxorpd %%ymm9, %%ymm9, %%ymm9\n\t"
"vxorpd %%ymm10, %%ymm10, %%ymm10\n\t"
"vxorpd %%ymm11, %%ymm11, %%ymm11\n\t"
".p2align 6\n\t"
"1:\n\t"
"prefetcht0 512(%[b])\n\t"
"prefetcht0 512(%[b1])\n\t"
"prefetcht0 512(%[b2])\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, %%ymm6\n\t"
"vmovapd 32(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"vmovapd 0(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm8\n\t"
"vmovapd 32(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm9\n\t"
"vmovapd 0(%[b2]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm4\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm10\n\t"
"vmovapd 32(%[b2]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm5\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm11\n\t"
"prefetcht0 576(%[b])\n\t"
"prefetcht0 576(%[b1])\n\t"
"prefetcht0 576(%[b2])\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, %%ymm6\n\t"
"vmovapd 96(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"vmovapd 64(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm8\n\t"
"vmovapd 96(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm9\n\t"
"vmovapd 64(%[b2]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm4\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm10\n\t"
"vmovapd 96(%[b2]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm5\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm11\n\t"
"prefetcht0 640(%[b])\n\t"
"prefetcht0 640(%[b1])\n\t"
"prefetcht0 640(%[b2])\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, %%ymm6\n\t"
"vmovapd 160(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"vmovapd 128(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm8\n\t"
"vmovapd 160(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm9\n\t"
"vmovapd 128(%[b2]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm4\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm10\n\t"
"vmovapd 160(%[b2]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm5\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm11\n\t"
"prefetcht0 704(%[b])\n\t"
"prefetcht0 704(%[b1])\n\t"
"prefetcht0 704(%[b2])\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, %%ymm6\n\t"
"vmovapd 224(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"vmovapd 192(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm8\n\t"
"vmovapd 224(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm9\n\t"
"vmovapd 192(%[b2]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm4\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm10\n\t"
"vmovapd 224(%[b2]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm5\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm11\n\t"
"add $32, %[a]\n\t"
"add $256, %[b]\n\t"
"add $256, %[b1]\n\t"
"add $256, %[b2]\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"
"vmovntpd %%ymm4, 128(%[out])\n\t"
"vmovntpd %%ymm5, 160(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm6, 0(%[out])\n\t"
"vmovntpd %%ymm7, 32(%[out])\n\t"
"vmovntpd %%ymm8, 64(%[out])\n\t"
"vmovntpd %%ymm9, 96(%[out])\n\t"
"vmovntpd %%ymm10, 128(%[out])\n\t"
"vmovntpd %%ymm11, 160(%[out])\n\t"
: [a] "+&r"(a),[b] "+&r"(b),[b1] "+&r"(b1),[b2] "+&r"(b2),[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();}
/* ==== d4be_ r12 trail (early decls; the body lives at the end of the engine) ==== */
static bool r12_arm=false;
static int r12_r=0,r12_ph=0,r12_off=0;
static void r12_emit(void);
static void r12_trail_body(void);
static inline void r12_trail(void){ if(r12_arm) r12_trail_body(); }
/* ===== 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();
/* d4be_: the r12 quadrant pass was emitted in-window by the trail; finish the rest here */
for(int i=0;i<0;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);
}
}
while(r12_r<512)r12_emit();
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,int sx,int sy){
if(sg>0) kernel< 1>(n,a,b,c,ldc,mr,pks,pky,pkd,sx,sy);
else if(sg<0) kernel<-1>(n,a,b,c,ldc,mr,pks,pky,pkd,sx,sy);
else kernel< 2>(n,a,b,c,ldc,mr,pks,pky,pkd,sx,sy);
}
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;
static_assert(BASE==512 && MC==30 && MR==6,"exact tail geometry");
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+=24) {
for(int jj=j;jj<j+24 && jj<n;jj+=NR)for(int ri=0;ri<MC;ri+=MR) {
const int r=(jj&NR)?MC-MR-ri:ri;
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,anext->ld,anext->ly);
} else
{ if(jj!=0 && (ri&2)) kernel_reverse(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)jj*(n+1),c+(size_t)(ic+r)*ldc+jj,ldc,MR); 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();r12_trail();}
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();}
}
if(j+24<=n)kernel_tail24(n,pkb[cur]+(size_t)MC*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+MC)*ldc+j,ldc);
else kernel2(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; /* final32 rows use30 full rows plus a dedicated2-row tail */
const int nch=rows2*2;                     /* 2 chunks per row            */
const int qd=n/2;                         /* doubles per chunk           */
static_assert(PERIOD==1,"one reconstruction opportunity per microtile");
const bool is_producer=pending.prog;
const int merge_limit=is_producer?BASE*pending.rowlim:BASE*BASE;
int merge_left=(TASKS&1)&&pending.active&&merge_limit>pending.pos?(merge_limit-pending.pos)/(STEP*(is_producer?2:1)):0;
int ck=0; int skip=(n/NR)*((mc+MR-1)/MR)-nch;
for (int j=0;j<n;j+=NR)for(int ri=0;ri<mc;ri+=MR) {
const int r=(j&NR)?mc-MR-ri:ri;
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) { if(j!=0 && (ri&2)) kernel_reverse(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 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,a.ld,a.ly);
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,a.ld,a.ly);
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,a.ld,a.ly);
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,a.ld,a.ly);
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,a.ld,a.ly);
if((TASKS&1) && merge_left>0){
--merge_left;pending_step_unchecked();if(is_producer)pending_step_unchecked();pending_prefetch();r12_trail();
}
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,const Mat *d4bd_succ=nullptr,bool d4bd_cin=false,bool d4be_r12src=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. */
Mat d4bd_d[5];
d4bd_d[0]=d4bd_first(aa[3],h);                 /* c1 carries into c2 */
d4bd_d[1]=d4bd_first(Mat{pa1,h},h);            /* c3 carries into c4 */
{ Mat d4bd_c5=d4bb_out?Mat{c,h}:Mat{c,ldc}; d4bd_d[2]=d4bd_first(d4bd_c5,h); }   /* c4 -> c5 */
d4bd_d[3]=d4bd_first(aa[0],h);                 /* c5 carries into c6 */
d4bd_d[4]=d4bd_first(Mat{pa2,h},h);            /* c6 carries into c7 */
rec<h>(h,aa[1],bb[2],out[1],ldc,sub,false,&d4bd_d[0],d4bd_cin);
d4ah_pre_fin_to(d4ah_n0+1);
rec<h>(h,aa[3],{pb2,h},out[3],ldc,sub,false,nullptr,true);   /* 2->3 NOT covered (pa0 belt) */
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,false,&d4bd_d[1],false);
rec<h>(h,{pa1,h},{pb1,h},out[5],h,sub,false,&d4bd_d[2],true);
if(d4bb_out) rec<h>(h,{c,h},{c+(size_t)h*h,h},out[6],h,sub,false,&d4bd_d[3],true);   /* d4bb_ s3 面板连续存放 */
else rec<h>(h,{c,ldc},{c+h,ldc},out[6],h,sub,false,&d4bd_d[3],true);
rec<h>(h,aa[0],bb[0],out[0],ldc,sub,false,&d4bd_d[4],true);
if(n==4096) {
root_pending={c,ldc,m,0,false};
rec<h>(h,{pa2,h},bb[3],out[2],ldc,sub,true,d4bd_succ,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,root_last,d4bd_succ,true,root_last);
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];}
}
const bool d4bd_F=((d4bd_np++)&1);   /* per-node parity flip */
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]:d4bd_succ;
const bool oddp=(((rr+d4bd_F)&1)!=0), np0=(rr==0 && !d4bd_cin);
if(rr==6) {
pending_finish();
pending={c,ldc,m,0,true,0,true};
if(d4be_r12src){r12_arm=true;r12_r=0;r12_ph=0;r12_off=0;}
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();
}
/* ==== d4be_ r12 trail body: emit the C12-quadrant endgame one row behind the
pending-combine frontier (see the file header for the derivation). ==== */
static void r12_emit(void){                        /* emit chunk [r12_off,r12_off+128) of row i */
int i=r12_r+(r12_ph?512:0);
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;
int j0=r12_off,j1=r12_off+128; if(j1>1024)j1=1024;
for(int j=j0;j<j1;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+64<1024){
_mm_prefetch((const char*)(mm+j+64),_MM_HINT_T0);
_mm_prefetch((const char*)(c12+j+64),_MM_HINT_T0);
_mm_prefetch((const char*)(pp+j+64),_MM_HINT_T0);
}
}
r12_off=j1;
if(r12_off>=1024){r12_off=0; if(r12_ph){r12_ph=0;r12_r++;} else r12_ph=1;}
}
static void r12_trail_body(void){
for(int rep=0;rep<8;rep++){
if(r12_r>=512)return;
if(pending.pos<(r12_r+1)*BASE)return;
int i=r12_r+(r12_ph?512:0);
int j1=r12_off+128; if(j1>1024)j1=1024;
if(root_pending.pos<=(long long)i*2048+1024+j1-1+64)return;
if(pend2.pos<=i*1024+j1-1+64)return;
r12_emit();
}
}

评测结果

Compilation

N/A

Compile OK

Score: N/A

显示更多

Compile OK

Testcase #1

1.838 s

502 MB + 332 KB

Accepted

Score: 100

显示更多

Time (ms): 1837.54829
Memory (KiB): 514380
Status: Run Finished

ok Accepted