// ===== REFERENCES =====
// [1] 本账号 saffah_cc_v41_agg1, 提交 #122731
// <https://duck.ac/submission/122731>
// 用途：**本提交的基座**。代码区与该件逐字节相同（本席只用定点字符串替换做改动，
// 未改任何算式、操作数、访存地址）。该件的完整 Credits 链
// （saffah_codex_6a_agg3 #122097 / #121415 / #121243 / ... ，pdoom #118506 / ... ，
// 本账号 #105696）原样保留在其正文注释里，未作删改。
//
// [2] 本舰队列题机口径读数（d4ah_ 席自测，探针件未提交）：
// * FMA 地板复核：`vfmadd231pd` 吞吐 = 2.0000 / 核拍（精确），延迟 4；
// 每迭代再加 8 条整数 ALU 仍是 2.0000 ⇒ 12 FMA/k 步 = 6.0 核拍 的地板无掩盖。
// * pre 相 = 纯 DRAM 限：真 pre(1024) 11.35e6 拍 ⇒ 5.9 B/核拍，与同机顺序读墙
// （5.95 B/核拍）持平；同地址形「只读」「只写」各 5.5 / 5.9 B/核拍 ⇒ 读写严格可加。
// * 传送带吸收率：在役 A 面板传送带每调用搬 2048 B，边际价 144.5 拍/调用
// （1x/6x/12x 剂量臂同发测得），同一份工作同步做约 530 拍 ⇒ 吸收比 3.7x。
//
// ===== 思路 =====
// 把上层 pre 相（判题机上 ~69 ms、完全暴露的纯 DRAM 相位）搬进内核相位，
// 用的是引擎现成的传送带，不新增任何汇编：
// * pre() 的 4 个输出本来就能写成 4 个「2 流 ± 融合」pass（每个只做 pre 里已有的
// 那一条 vaddpd/vsubpd，算式与顺序逐条相同 ⇒ 逐位等价）：
// p = c + d ; q = p - a ; r = b - q ; t = a - c （左半边；右半边同构）
// * 这 4 条正是 kernel<SG> 里那条传送带的形状（load [pks+pkd] / add|sub [pky+pkd] /
// store [pkd]），而一次 kernel 调用的行程恰是 2048 B。
// * 于是每个 rec<2h> 帧在入口武装 8 个 pass 任务（每侧 4 个，栈序按各自截止期排：
// 弹出序 = pb0,pb1,pb2,pb3,pa0,pa1,t,pa2），块在叶子的调用循环里发放（与 A 面板
// 传送带共用调用槽，A 面板块优先）；每个消费子节点前用 pre_fin_to() 把该面板
// 自己的兜底做完（child3 前 pb2、child4 前 pa0、child5 前 pa1、child6 前 t、
// child2 前 pa2）⇒ 每个面板都拿到自己的窗口，任何泵送速率下都正确。
// 为什么划算：内核环在 p01 上满（96 FMA / 48 拍），但载入口只用 72/96 槽、p4 空转，
// 而叶相期间 DRAM 近乎全空（~1.5 GB/s vs 21 GB/s 墙）⇒ 传送带每 2048 B 只买 144.5 拍，
// 是同步价的 1/3.7；pre 的字节搬进叶相后不再单独占时间。
// 判题机切片实测（rec<2048>，同二进制只翻一个 static 开关）：-0.84%，逐位 chk 相同；
// 全跑 360448 块 = 738 MB 输出 = 理论值。v1（两侧各 1 个任务）全跑实测 1869.26 ms
// （-7.81 ms）—— 模型显示 root 帧的 131072 块被挤进单个子窗口、每块贵到 ~908 拍
// （帧级只要 ~450）⇒ v2 把任务拆到 pass 粒度、逐面板给窗口。
//
// 与在册封条的关系：不碰 packB（d4ag_ 封）、不碰叶级分解（de Groote 唯一性）、
// 不搬 pre 进钩子（d4u_ 封的是「钩子泵送」形态：钩子步延迟暴露；本件走的是内核内环
// 传送带，吃的是端口余量而非钩子槽）。
// =====================================================================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6a_agg3, 提交 #122097
// <https://duck.ac/submission/122097>
// 用途：**本提交的基座**。本文件代码区与该件**逐字节相同**（本席用
// `python3 -c` 剥去全部注释后逐行比对，diff = 0）；该件自身的
// Credits/References 链（#121415/#121243/#121181/#121167/#121142/#120888/
// #120615/#120547/#120521/#120514/pdoom #118506/#119315/#105696/#104821/...）
// 原样保留在其正文注释里，未作删改。
//
// ===== 思路 =====
// 本件是**对手最新件的取回件**（协调者 §2.18.349「抄对手」通道），目的 = 把本账号的
// 在册最好件从 #105696 换成对手当前榜首件，作为后续加刀的新基座。
//
// 我方在册件的沿革（本席核对过）：
// * #105696 = 对手 #104821 正文 + 我方的一把刀（删掉 `rec(n>2*BASE)` 里位置 1 之后
// 的那处 `pend2_finish()` 屏障，把上层延迟 combine 的窗口跨兄弟节点延长），
// 判题机 1902.881436 ms。
// * 对手随后把**我方这把刀吸收进了他们的件**（#122097 的 Credits 明写
// "saffah_cc_v41_agg1, https://duck.ac/submission/105696: immediate engine base;
// safely extend the upper deferred-combine window across siblings"），并在其上叠了
// 七步新改动（6×8 内核与上层几何的 4096 特化 / 递归尺寸特化 / 叶子任务级特化 /
// 根 M3 子树内的根重建融合 …），判题机 1876.693695 ms。
// ⇒ 我方 #105696 相对它**没有任何落后量**（我方唯一的刀已在其中），只有对手的
// 七步净改进（−26.188 ms = −1.38%）。
//
// 判定（本席判题机口径实测，`tools/tryrun.py`，n=2048 尺，49 叶，双胞胎臂噪声 0.01%）：
// * 相位表：leaf（含 mikro + packA 传送带 + pending/pend2 重建）**93.7%** ·
// packB(d1a2_pk) **3.7%** · pre **2.3%** · 其余 0.3%；
// * 内核本征速率：2.78M k-步 / 叶，实测 **6.41 拍/k 步**，而 FMA 端口地板
// （12 FMA + 0.125 打包 vaddpd → /2 端口）= **6.19 拍/k** ⇒ 内核已在端口墙的
// **96.6%**；且 `k_bsame`（全部 b 载入压到同一行，取数流量降到 1/64）实测
// **+0.02% ⇒ 内核不受访存限制**，`k_nopf`（删 8 条 prefetcht0）+0.15% ⇒ 不受
// 前端限制。（⇒ 本题的 mikro 无刀可下，与全队「mikro 已到端口墙」一致。）
// * ⇒ 本件的价值 = 白捡对手七步（−26.19 ms），把缺口从 44.95 ms 收到 ~18.8 ms。
//
// 合规：引用已列原作者 + 提交 URL + 所用内容（其正文全部）；未改动算术、未加特例、
// 未触碰任何测试点判据；无 /tmp 产物；本件为**真拼写差异**件（不是纯锚）。
// =====================================================================
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/121415: retained frontier-gated FP64 matrix engine and all inherited citations below.
* Idea: Express the two carried-packing source addresses as constant displacements from the packing destination. Advance only the destination pointer inside the FMA loop, replacing two pointer additions per eight-K step with x86 indexed loads. Values, FMA order and packing schedule are unchanged.
* Purpose: Reduce integer updates and live output operands in the carried-packing microkernel.
*/
#include <cstdint>
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/121243,
reused the matrix engine and inherited source credits below. No independent
license declared.
Idea: Gate progressive reconstruction once at each kernel call. A row frontier
is a multiple of 128 elements, and paired 64-element steps also start at that
alignment, so a permitted pair cannot cross the frontier. Use unchecked steps
inside that gate; stalled calls skip both step checks and the prefetch check.
Purpose: Reduce repeated bounds checks while reconstruction waits for more rows.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/121181,
reused the matrix engine and all inherited source references below.
No independent license declared.
Idea: Specialize recursion for the known 4096, 2048, 1024 and 512 sizes. Each
level can discard unreachable size branches and use constant transform strides
and extents, preserving the same product order and reconstruction schedule.
Purpose: Remove dynamic size arithmetic and dispatch above the matrix leaves.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/121167,
reused the matrix engine and all inherited references below. No independent
license declared.
Idea: Every recursive multiplication in this fixed 4096-square task terminates
at a 512-square leaf. Compile that leaf extent as a constant to simplify packing
bounds, panel strides, tail selection and streaming-store conditions.
Purpose: Remove general-size arithmetic from the remaining leaf hot loops.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/121142,
reused the matrix engine and all inherited references below. No independent
license declared.
Idea: Specialize leaf reconstruction independently for each of the three
levels active at entry. Deferred tasks can retire but cannot start within a
leaf; omit checks for the levels that were absent when this leaf began.
Purpose: Extend the measured benefit of inactive-task specialization to leaves
with only one or two active reconstruction levels.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/120888,
reused the current pipelined matrix engine and all inherited source credits below.
No independent license declared.
Idea: Select a leaf implementation once based on whether deferred reconstruction
is active. No reconstruction task can begin inside a leaf, so a leaf entering
with none can omit every per-microkernel task check and row-frontier update.
The other implementation retains all existing task progression and ordering.
Purpose: Remove repeated inactive-task dispatch from most microkernel calls.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/120615,
reused the current matrix engine and plain-kernel software B-load pipeline.
All earlier source attributions remain below; no independent license declared.
Idea: Apply the one-step B-load pipeline to the A-packing microkernel too.
Its spare register is only needed for packing after the last inner FMA, so
alternating that register with the current B vector preserves all values.
Purpose: Hide B load latency in the remaining carried-packing kernel calls.
*/
/* References: pdoom https://duck.ac/submission/118506: reused FP64
recursive matrix engine. saffah_cc_v41_agg1 https://duck.ac/submission/119315:
reused plain-kernel B-load software pipeline and deferred prefetch scheduling.
saffah_codex_6a_agg3 https://duck.ac/submission/120547: B-panel reuse and merge64.
All inherited source credits retained below; no independent licenses declared.
Idea: Apply the one-step-ahead B-load scheduling to the larger matrix kernel,
using the fixed-size output path and preserving each FMA accumulation order.
Purpose: Hide B load-to-use latency across more than100000 microkernel calls.
*/
/* References:
[1] saffah_codex_6a_agg3, https://duck.ac/submission/120521.
Reused the complete FP64 matrix engine; all inherited source citations remain
below. No separate license is declared on the public source.
Idea: Reconstruct64 positions per bottom merge step instead of32. Both divide
the512-wide rows, preserving frontiers and operation grouping while amortizing
step checks and cursor updates across twice as much contiguous memory traffic.
Purpose: Measure lower reconstruction bookkeeping in the fused compute schedule.
*/
/* References:
[1] pdoom, https://duck.ac/submission/118506.
Reused the4k FP64 engine and its deferred root reconstruction, preserving its
complete inherited attribution list below. No separate license is declared.
[2] saffah_cc_v41_agg1, https://duck.ac/submission/119315.
Copied the1k fused packed-B helpers and two-reuse leaf ordering.
[3] saffah_codex_6a_agg3, https://duck.ac/submission/120514.
Reused our validated2k integration of those helpers; no separate license.
Idea: Apply two packed-B reuse positions to every bottom group of the4k
recursion, removing repeated source traversal when a pure panel is still packed.
Purpose: Measure the accumulated traffic reduction at the largest FP64 size.
*/
/*
Credits for inherited implementation and ideas:
- saffah_cc_v41_agg1, https://duck.ac/submission/105696: immediate engine
base; safely extend the upper deferred-combine window across siblings.
- saffah_codex_6s_agg2, https://duck.ac/submission/104821: preceding engine,
compiler model and upper workspace layout, including its inherited work.
- saffah_cc_v41_agg1, https://duck.ac/submission/101193: two-level deferred
reconstruction, product order and panel-layout improvements.
- saffah_cc_v41_agg1, https://duck.ac/submission/103010: produced-row
frontier and early reconstruction; saffah_codex_6s_agg2,
https://duck.ac/submission/103224: its4k integration and leaf ordering.
- saffah_codex_6s_agg2, https://duck.ac/submission/103194: eight-K plain
microkernel; https://duck.ac/submission/103275: power-of-two addressing;
https://duck.ac/submission/104809: leaf scratch-plane layout.
- saffah_cc_v41_agg1, https://duck.ac/submission/100776: dedicated two-row
tail; https://duck.ac/submission/100511 and https://duck.ac/submission/100475:
carry next-A packing in the FMA loop, double-buffer A and grouped B packing.
- saffah_cc_v41_260924, https://duck.ac/submission/99846: fused operand sums
and packing, one-pass sums and reconstruction.
- saffah_cc_v41_agg1, https://duck.ac/submission/99875 and
https://duck.ac/submission/99950: non-temporal intermediate stores.
- saffah_cc_v41_260924, https://duck.ac/submission/99141, and
saffah_cc_v41_agg1, https://duck.ac/submission/99167:6x8 AVX2 blocking,
explicit vector memory instructions, prefetch and packing-layout baseline.
- Earlier pdoom work: hybrid Winograd/Strassen recurrence, output-backed
scratch placement, delayed reconstruction, merged32-row tail and2x16
kernel, and page-spacing experiments.
New changes in this version:
1. Specialize the full6x8 kernels and hot upper reconstruction geometry for
the fixed4096 problem. Two-row output has its separate2x16 kernel.
2. The root's last product M3 only contributes to C12. Once the preceding
M1 is ready, reconstruct C11,C21,C22 and the C12 base (M1+M6)+M5 during
M3's remaining five1024 subproblems. Retired M5 storage holds that base.
3. Add the base while reconstructing M3's four quadrants. The root task
processes32 elements per step only when bottom reconstruction is inactive.
This moves320MiB of logical root traffic into the last multiplication
subtree; total nominal traffic and arithmetic count do not decrease.
All arithmetic remains binary64, with the original expression grouping.
No input quantization. Specialized to the problem's n=4096 contract.
*/
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing")
#include <sys/mman.h>
#pragma GCC target("avx2,fma,arch=skylake")
#include <immintrin.h>
#ifndef B2M_NO_BARA
#define B2M_NO_BARA 1
#endif
#define MD9_T2BASE 1
#include <stddef.h>
#ifndef BASE
#define BASE 512
#endif
#ifndef MC
#define MC 30
#endif
#ifndef NTMIN
#define NTMIN 512
#endif
#ifndef UPPER_PAGE_PAD
#define UPPER_PAGE_PAD 15872
#endif
#ifndef LEAF_PAD
#define LEAF_PAD 1152
#endif
#ifndef SUMNT
#define SUMNT 1024
#endif
static constexpr int MR=6,NR=8;
alignas(4096) static double workspace[48*1024*1024];
alignas(64) static double ap[(MC+MR)*(BASE+24)],ap2[(MC+MR)*(BASE+24)],bp[BASE*(BASE+24)];
/* m2s/md5 fused-pack support: a zero row (for pad rows and for pure copies) and the
two sign masks (+0.0 / -0.0) used to apply Strassen's operand sign in-register. */
alignas(64) static double pk_zero[BASE+8];
alignas(64) static double pk_mp[4]={0.0,0.0,0.0,0.0}, pk_mn[4]={-0.0,-0.0,-0.0,-0.0};
struct Mat { const double *x; int ld; const double *y=nullptr; int ly=0; int sign=0; };
static inline __m256d ld(Mat a,int i,int j) {
__m256d x=_mm256_load_pd(a.x+(size_t)i*a.ld+j);
if(a.sign) {
__m256d y=_mm256_load_pd(a.y+(size_t)i*a.ly+j);
x=a.sign>0?_mm256_add_pd(x,y):_mm256_sub_pd(x,y);
}
return x;
}
static inline double scalar(Mat a,int i,int k) {
double x=a.x[(size_t)i*a.ld+k];
if(a.sign)x+=a.sign>0?a.y[(size_t)i*a.ly+k]:-a.y[(size_t)i*a.ly+k];
return x;
}
template<int SG> static __attribute__((always_inline)) inline void kernel(int n,const double *a,const double *b,double *c,int ldc,int rows, const double *pks,const double *pky,double *pkd) {
double *out=c;
ptrdiff_t stride=(ptrdiff_t)ldc*8;
int k=n;
pks=(const double*)((uintptr_t)pks-(uintptr_t)pkd);
pky=(const double*)((uintptr_t)pky-(uintptr_t)pkd);
asm volatile (
"vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
"vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
"vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
"vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
"vxorpd %%ymm4, %%ymm4, %%ymm4\n\t"
"vxorpd %%ymm5, %%ymm5, %%ymm5\n\t"
"vxorpd %%ymm6, %%ymm6, %%ymm6\n\t"
"vxorpd %%ymm7, %%ymm7, %%ymm7\n\t"
"vxorpd %%ymm8, %%ymm8, %%ymm8\n\t"
"vxorpd %%ymm9, %%ymm9, %%ymm9\n\t"
"vxorpd %%ymm10, %%ymm10, %%ymm10\n\t"
"vxorpd %%ymm11, %%ymm11, %%ymm11\n\t"
"vmovapd 0(%[b]), %%ymm12\n\t"
".p2align 6\n\t"
"1:\n\t"
"vmovapd 64(%[b]), %%ymm15\n\t"
"vmovapd 32(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 128(%[b]), %%ymm12\n\t"
"vmovapd 96(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 512(%[b])\n\t"
"prefetcht0 576(%[b])\n\t"
"prefetcht0 640(%[b])\n\t"
"prefetcht0 704(%[b])\n\t"
"vmovapd 192(%[b]), %%ymm15\n\t"
"vmovapd 160(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 256(%[b]), %%ymm12\n\t"
"vmovapd 224(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 768(%[b])\n\t"
"prefetcht0 832(%[b])\n\t"
"prefetcht0 896(%[b])\n\t"
"prefetcht0 960(%[b])\n\t"
"vmovapd 320(%[b]), %%ymm15\n\t"
"vmovapd 288(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 384(%[b]), %%ymm12\n\t"
"vmovapd 352(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 448(%[b]), %%ymm15\n\t"
"vmovapd 416(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 512(%[b]), %%ymm12\n\t"
"vmovapd 480(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovupd 0(%[pkd],%[pks]), %%ymm15\n\t"
".if %c[sg] == 1\n\t"
"vaddpd 0(%[pkd],%[pky]), %%ymm15, %%ymm15\n\t"
".elseif %c[sg] == -1\n\t"
"vsubpd 0(%[pkd],%[pky]), %%ymm15, %%ymm15\n\t"
".endif\n\t"
"vmovupd %%ymm15, 0(%[pkd])\n\t"
"add $32, %[pkd]\n\t"
"add $64, %[a]\n\t"
"add $512, %[b]\n\t"
"sub $8, %[k]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, (%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm2, (%[out])\n\t"
"vmovntpd %%ymm3, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm4, (%[out])\n\t"
"vmovntpd %%ymm5, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm6, (%[out])\n\t"
"vmovntpd %%ymm7, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm8, (%[out])\n\t"
"vmovntpd %%ymm9, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm10, (%[out])\n\t"
"vmovntpd %%ymm11, 32(%[out])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k), [pkd] "+&r"(pkd)
: [pks] "r"(pks), [pky] "r"(pky), [sg] "i"(SG), [stride] "r"(stride), [o0] "i"(0*(BASE+24)*8), [o1] "i"(1*(BASE+24)*8), [o2] "i"(2*(BASE+24)*8), [o3] "i"(3*(BASE+24)*8), [o4] "i"(4*(BASE+24)*8), [o5] "i"(5*(BASE+24)*8)
: "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
static __attribute__((always_inline)) inline void kernel_plain(int n,const double *a,const double *b,double *c,int ldc,int rows) {
double *out=c;
ptrdiff_t stride=(ptrdiff_t)ldc*8;
int k=n;
asm volatile (
"vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
"vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
"vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
"vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
"vxorpd %%ymm4, %%ymm4, %%ymm4\n\t"
"vxorpd %%ymm5, %%ymm5, %%ymm5\n\t"
"vxorpd %%ymm6, %%ymm6, %%ymm6\n\t"
"vxorpd %%ymm7, %%ymm7, %%ymm7\n\t"
"vxorpd %%ymm8, %%ymm8, %%ymm8\n\t"
"vxorpd %%ymm9, %%ymm9, %%ymm9\n\t"
"vxorpd %%ymm10, %%ymm10, %%ymm10\n\t"
"vxorpd %%ymm11, %%ymm11, %%ymm11\n\t"
"vmovapd 0(%[b]), %%ymm12\n\t"
".p2align 6\n\t"
"1:\n\t"
"vmovapd 64(%[b]), %%ymm15\n\t"
"vmovapd 32(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 128(%[b]), %%ymm12\n\t"
"vmovapd 96(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 512(%[b])\n\t"
"prefetcht0 576(%[b])\n\t"
"prefetcht0 640(%[b])\n\t"
"prefetcht0 704(%[b])\n\t"
"vmovapd 192(%[b]), %%ymm15\n\t"
"vmovapd 160(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 256(%[b]), %%ymm12\n\t"
"vmovapd 224(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 768(%[b])\n\t"
"prefetcht0 832(%[b])\n\t"
"prefetcht0 896(%[b])\n\t"
"prefetcht0 960(%[b])\n\t"
"vmovapd 320(%[b]), %%ymm15\n\t"
"vmovapd 288(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 384(%[b]), %%ymm12\n\t"
"vmovapd 352(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 448(%[b]), %%ymm15\n\t"
"vmovapd 416(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 512(%[b]), %%ymm12\n\t"
"vmovapd 480(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"add $64, %[a]\n\t"
"add $512, %[b]\n\t"
"sub $8, %[k]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, (%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm2, (%[out])\n\t"
"vmovntpd %%ymm3, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm4, (%[out])\n\t"
"vmovntpd %%ymm5, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm6, (%[out])\n\t"
"vmovntpd %%ymm7, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm8, (%[out])\n\t"
"vmovntpd %%ymm9, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm10, (%[out])\n\t"
"vmovntpd %%ymm11, 32(%[out])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k)
: [stride] "r"(stride), [o0] "i"(0*(BASE+24)*8), [o1] "i"(1*(BASE+24)*8), [o2] "i"(2*(BASE+24)*8), [o3] "i"(3*(BASE+24)*8), [o4] "i"(4*(BASE+24)*8), [o5] "i"(5*(BASE+24)*8)
: "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
/* ===== mdA_dp2: deferred h=1024 combine (Winograd form, exact same arithmetic) ===== */
struct Pend2 { double *c; int ldc; const double *m; int h; int pz; int pos; bool active; const double *add=nullptr; };
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;
int i=pend2.pos >> __builtin_ctz((unsigned)H), j=pend2.pos & (H-1);
const double *c=pend2.c+(size_t)i*pend2.ldc+j;
const double *m=pend2.m+pend2.pos;
for(int t=0;t<STEP2;t+=8) {
_mm_prefetch((const char*)(c+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+H+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)H*pend2.ldc+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)H*pend2.ldc+H+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+(H*H+UPPER_PAGE_PAD)+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+2*(H*H+UPPER_PAGE_PAD)+t),_MM_HINT_T1);
}
}
template<int H=1024,bool ADD=false> static inline void pend2_step(void){
int i=pend2.pos >> __builtin_ctz((unsigned)H), j=pend2.pos & (H-1);
for(int t=0;t<STEP2;t+=4) {
double *c11=pend2.c+(size_t)i*pend2.ldc+j+t,*c12=c11+H,
*c21=c11+(size_t)H*pend2.ldc,*c22=c21+H;
const double *m=pend2.m+pend2.pos+t;
__m256d m1=_mm256_load_pd(c11),m2=_mm256_load_pd(c21),
m3=_mm256_load_pd(c12),m4=_mm256_load_pd(c22),
m5=_mm256_load_pd(m),m6=_mm256_load_pd(m+(H*H+UPPER_PAGE_PAD)),m7=_mm256_load_pd(m+2*(H*H+UPPER_PAGE_PAD));
__m256d u2=_mm256_add_pd(m1,m6),u3=_mm256_add_pd(u2,m7),u4=_mm256_add_pd(u2,m5);
__m256d r11=_mm256_add_pd(m1,m2),r12=_mm256_add_pd(u4,m3),
r21=_mm256_sub_pd(u3,m4),r22=_mm256_add_pd(u3,m5);
if(ADD) {
const double *p=pend2.add+(size_t)i*2048+j+t;
r11=_mm256_add_pd(_mm256_load_pd(p),r11);
r12=_mm256_add_pd(_mm256_load_pd(p+H),r12);
r21=_mm256_add_pd(_mm256_load_pd(p+(size_t)H*2048),r21);
r22=_mm256_add_pd(_mm256_load_pd(p+(size_t)H*2048+H),r22);
}
_mm256_store_pd(c11,r11);_mm256_store_pd(c12,r12);
_mm256_store_pd(c21,r21);_mm256_store_pd(c22,r22);
}
pend2.pos+=STEP2;
if(pend2.pos>=H*H)pend2.active=false;
}
static void pend2_finish(void){
if(pend2.add) {while(pend2.active)pend2_step<1024,true>();}
else if(pend2.h==1024) {while(pend2.active)pend2_step<1024>();}
else {while(pend2.active)pend2_step<2048>();}
}
struct Pending { double *c; int ldc; const double *m; int pos=0; bool active=false; int rowlim=0; bool prog=false; };
static Pending pending;
alignas(4096) static double pending_work[3*(BASE*BASE+LEAF_PAD)+LEAF_PAD];
static unsigned pending_tick=0;
static constexpr int STEP=64, PERIOD=1;
static inline void pending_prefetch() {
if(!pending.active)return;
if(pending.prog && pending.pos/BASE>=pending.rowlim)return;
int i=pending.pos/BASE,j=pending.pos%BASE;
const double *c=pending.c+(size_t)i*pending.ldc+j;
const double *m=pending.m+pending.pos;
for(int t=0;t<STEP;t+=8) {
_mm_prefetch((const char*)(c+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+BASE+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)BASE*pending.ldc+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)BASE*pending.ldc+BASE+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+(BASE*BASE+LEAF_PAD)+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+2*(BASE*BASE+LEAF_PAD)+t),_MM_HINT_T1);
}
}
static inline void pending_step_unchecked() {
int i=pending.pos/BASE,j=pending.pos%BASE;
for(int t=0;t<STEP;t+=4) {
double *c11=pending.c+(size_t)i*pending.ldc+j+t,
*c12=c11+BASE,*c21=c11+(size_t)BASE*pending.ldc,*c22=c21+BASE;
const double *m=pending.m+pending.pos+t;
__m256d m1=_mm256_load_pd(c11),m2=_mm256_load_pd(c21),
m3=_mm256_load_pd(c12),m4=_mm256_load_pd(c22),
m5=_mm256_load_pd(m),m6=_mm256_load_pd(m+(BASE*BASE+LEAF_PAD)),m7=_mm256_load_pd(m+2*(BASE*BASE+LEAF_PAD));
__m256d r11=_mm256_add_pd(_mm256_sub_pd(_mm256_add_pd(m1,m4),m5),m7),
r12=_mm256_add_pd(m3,m5),r21=_mm256_add_pd(m2,m4),
r22=_mm256_add_pd(_mm256_add_pd(_mm256_sub_pd(m1,m2),m3),m6);
_mm256_store_pd(c11,r11);_mm256_store_pd(c12,r12);
_mm256_store_pd(c21,r21);_mm256_store_pd(c22,r22);
}
pending.pos+=STEP;
if(pending.pos==BASE*BASE)pending.active=false;
}
static inline void pending_step(){if(pending.prog && pending.pos/BASE>=pending.rowlim)return;pending_step_unchecked();}
static void pending_finish() {
pending.prog=false;
while(pending.active)pending_step();
}
static __attribute__((always_inline)) inline void kernel2(int n,const double *a,const double *b,double *c,int ldc) {
int k=n;
double *out=c;
ptrdiff_t stride=(ptrdiff_t)ldc*8;
asm volatile (
"vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
"vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
"vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
"vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
".p2align 6\n\t" "1:\n\t"
"prefetcht0 512(%[b])\n\t"
"prefetcht0 576(%[b])\n\t"
"prefetcht0 640(%[b])\n\t"
"prefetcht0 704(%[b])\n\t"
"vmovapd 0(%[b]), %%ymm12\n\t"
"vmovapd 32(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vmovapd 64(%[b]), %%ymm12\n\t"
"vmovapd 96(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vmovapd 128(%[b]), %%ymm12\n\t"
"vmovapd 160(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vmovapd 192(%[b]), %%ymm12\n\t"
"vmovapd 224(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"add $32, %[a]\n\t"
"add $256, %[b]\n\t"
"sub $4, %[k]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, (%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm2, (%[out])\n\t"
"vmovntpd %%ymm3, 32(%[out])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k)
: [stride] "r"(stride), [o0] "i"(0*(BASE+24)*8), [o1] "i"(1*(BASE+24)*8)
: "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm12","ymm13","ymm14");
}
static __attribute__((noinline)) void kernel_tail16(int n,const double *a,const double *b,double *out,int ldc) {
const double *b1=b+8*(n+1);
ptrdiff_t stride=(ptrdiff_t)ldc*8;
int kk=n;
asm volatile (
"vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
"vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
"vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
"vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
"vxorpd %%ymm4, %%ymm4, %%ymm4\n\t"
"vxorpd %%ymm5, %%ymm5, %%ymm5\n\t"
"vxorpd %%ymm6, %%ymm6, %%ymm6\n\t"
"vxorpd %%ymm7, %%ymm7, %%ymm7\n\t"
".p2align 6\n\t"
"1:\n\t"
"prefetcht0 512(%[b])\n\t"
"prefetcht0 512(%[b1])\n\t"
"vbroadcastsd %c[o0]+0(%[a]), %%ymm12\n\t"
"vbroadcastsd %c[o1]+0(%[a]), %%ymm13\n\t"
"vmovapd 0(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm0\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm4\n\t"
"vmovapd 32(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm5\n\t"
"vmovapd 0(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm6\n\t"
"vmovapd 32(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"prefetcht0 576(%[b])\n\t"
"prefetcht0 576(%[b1])\n\t"
"vbroadcastsd %c[o0]+8(%[a]), %%ymm12\n\t"
"vbroadcastsd %c[o1]+8(%[a]), %%ymm13\n\t"
"vmovapd 64(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm0\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm4\n\t"
"vmovapd 96(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm5\n\t"
"vmovapd 64(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm6\n\t"
"vmovapd 96(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"prefetcht0 640(%[b])\n\t"
"prefetcht0 640(%[b1])\n\t"
"vbroadcastsd %c[o0]+16(%[a]), %%ymm12\n\t"
"vbroadcastsd %c[o1]+16(%[a]), %%ymm13\n\t"
"vmovapd 128(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm0\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm4\n\t"
"vmovapd 160(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm5\n\t"
"vmovapd 128(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm6\n\t"
"vmovapd 160(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"prefetcht0 704(%[b])\n\t"
"prefetcht0 704(%[b1])\n\t"
"vbroadcastsd %c[o0]+24(%[a]), %%ymm12\n\t"
"vbroadcastsd %c[o1]+24(%[a]), %%ymm13\n\t"
"vmovapd 192(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm0\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm4\n\t"
"vmovapd 224(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm5\n\t"
"vmovapd 192(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm6\n\t"
"vmovapd 224(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"add $32, %[a]\n\t"
"add $256, %[b]\n\t"
"add $256, %[b1]\n\t"
"sub $4, %[kk]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, 0(%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"vmovntpd %%ymm2, 64(%[out])\n\t"
"vmovntpd %%ymm3, 96(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm4, 0(%[out])\n\t"
"vmovntpd %%ymm5, 32(%[out])\n\t"
"vmovntpd %%ymm6, 64(%[out])\n\t"
"vmovntpd %%ymm7, 96(%[out])\n\t"
: [a] "+&r"(a),[b] "+&r"(b),[b1] "+&r"(b1),[out] "+&r"(out),[kk] "+&r"(kk)
: [stride] "r"(stride),[o0] "i"(0),[o1] "i"((BASE+24)*8)
: "memory","cc","ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
// Root M3 is last and contributes only to C12. The other three outputs and
// the C12 base (M1+M6)+M5 are independent of M3.
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();}
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 相传送带吸收（机制；每 pass 一个任务 + 按截止期兜底） ===== */
static void pre(int h,Mat a,double *s1,double *s2,double *s4,double *s3,int ld3,bool right);
struct D4ahPT { int h; int i, j; const double *s1; int l1; const double *s2; int l2; double *d; int l3; int sg; };
static D4ahPT d4ah_pt[24];
static int d4ah_ptn=0;
static int d4ah_on=1;
static long long d4ah_pchunks=0;
static void d4ah_push(int h,const double *s1,int l1,const double *s2,int l2,double *d,int l3,int sg){
if(d4ah_ptn>=24) return;
D4ahPT *t=&d4ah_pt[d4ah_ptn++];
t->h=h;t->i=0;t->j=0;t->s1=s1;t->l1=l1;t->s2=s2;t->l2=l2;t->d=d;t->l3=l3;t->sg=sg;
}
static void d4ah_run(D4ahPT *t){
int i0=t->i, j0=t->j;
for(int i=i0;i<t->h;i++){
const double *r1=t->s1+(size_t)i*t->l1, *r2=t->s2+(size_t)i*t->l2;
double *rr=t->d+(size_t)i*t->l3;
int jb=(i==i0)?j0:0;
for(int j=jb;j<t->h;j+=4){
__m256d u=_mm256_loadu_pd(r1+j), v=_mm256_loadu_pd(r2+j);
_mm256_storeu_pd(rr+j, t->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;
if(right){
d4ah_push(h, s2,h, x+(size_t)h*ld,ld, s4,h, -1); /* pb2 = q - c -> child3 */
d4ah_push(h, x+(size_t)h*ld+h,ld, x+h,ld, s3,ld3, -1); /* t = d - b -> child6 */
d4ah_push(h, x+(size_t)h*ld+h,ld, s1,h, s2,h, -1); /* pb1 = d - p -> child5 */
d4ah_push(h, x+h,ld, x,ld, s1,h, -1); /* pb0 = b - a -> child4 */
} else {
d4ah_push(h, x+h,ld, s2,h, s4,h, -1); /* pa2 = b - q -> child2(末) */
d4ah_push(h, x,ld, x+(size_t)h*ld,ld, s3,ld3, -1); /* t = a - c -> child6 */
d4ah_push(h, s1,h, x,ld, s2,h, -1); /* pa1 = p - a -> child5 */
d4ah_push(h, x+(size_t)h*ld,ld, x+(size_t)h*ld+h,ld, s1,h, +1); /* pa0 = c + d -> child4 */
}
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];
*q1=t->s1+(size_t)t->i*t->l1+t->j;
*q2=t->s2+(size_t)t->i*t->l2+t->j;
*qd=t->d +(size_t)t->i*t->l3+t->j;
*qsg=t->sg;
d4ah_pchunks++;
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--; }
}
template<int TASKS> static void leaf_impl(int leaf_n,Mat a,Mat b,double *c,int ldc) {
constexpr int n=BASE;
d1a2_pk(n,b,bp);
if constexpr(TASKS&1)pending_prefetch();
double *pkb[2]={ap,ap2};
int cur=0;
{ /* 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[0]+(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[0]+(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) {
for(int j=0;j<n;j+=16) {
for(int jj=j;jj<j+16;jj+=NR)for(int r=0;r<MC;r+=MR) {
kernel_plain(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)jj*(n+1),c+(size_t)(ic+r)*ldc+jj,ldc,MR);
if((TASKS&1) && pending.active && (++pending_tick%PERIOD)==0 && (!pending.prog || pending.pos<BASE*pending.rowlim)) {pending_step_unchecked();if(pending.prog)pending_step_unchecked();pending_prefetch();}
if((TASKS&2) && pend2.active) {pend2_step();pend2_prefetch();}
if((TASKS&4) && root_pending.active && !pending.active) {root_step();root_prefetch();}
}
kernel_tail16(n,pkb[cur]+(size_t)MC*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+MC)*ldc+j,ldc);
}
ic+=mc;
if((TASKS&1) && pending.prog)pending.rowlim=ic;
continue;
}
/* ---- the NEXT ic-block's A panel is carried by this block's kernel calls ---- */
const int ic2=ic+mc;
const int nxt=(ic2<n)?(n-ic2<=MC+MR?n-ic2:MC):0;
const int rows2=((nxt+MR-1)/MR)*MR; /* rows incl. the zero pad rows */
const int nch=rows2*2; /* 2 chunks per row */
const int qd=n/2; /* doubles per chunk */
int ck=0; int skip=(n/NR)*((mc+MR-1)/MR)-nch;
for(int j=0;j<n;j+=NR)for(int r=0;r<mc;r+=MR) {
if(mc-r==2) { /* md9_t2: 2-row tail body (was 6-row body + tmp copy) */
kernel2(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc);
continue;
}
const double *pks=pk_zero,*pky=pk_zero,*pkm=pk_mp;
double *pkd=pk_zero;
int pkc=0; int psg=0;
if(ck<nch && --skip<0) {
const int row=ck/2, off=(ck%2)*qd;
if(row<nxt) {
pks=a.x+(size_t)(ic2+row)*a.ld+(size_t)off;
pky=a.sign?(a.y+(size_t)(ic2+row)*a.ly+(size_t)off):pk_zero;
pkm=a.sign>0?pk_mp:pk_mn;
}
pkd=pkb[cur^1]+(size_t)row*(BASE+24)+(size_t)off;
pkc=8;
ck++;
}
if(!pkc && d4ah_ptn){
if(d4ah_pre_chunk(&pks,&pky,&pkd,&psg)) pkc=1;
}
if(!pkc) kernel_plain(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR);
else if(psg>0) kernel<1>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,pks,pky,pkd);
else if(psg<0) kernel<-1>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,pks,pky,pkd);
else if(a.sign==0) kernel<2>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,pks,pky,pkd);
else if(a.sign>0) kernel<1>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,pks,pky,pkd);
else kernel<-1>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,pks,pky,pkd);
if((TASKS&1) && pending.active && (++pending_tick%PERIOD)==0 && (!pending.prog || pending.pos<BASE*pending.rowlim)) {
pending_step_unchecked(); if(pending.prog)pending_step_unchecked(); pending_prefetch();
}
if((TASKS&2) && pend2.active) { 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){int mode=pending.active+2*pend2.active+4*root_pending.active;switch(mode){case 0:leaf_impl<0>(n,a,b,c,ldc);break;case 1:leaf_impl<1>(n,a,b,c,ldc);break;case 2:leaf_impl<2>(n,a,b,c,ldc);break;case 3:leaf_impl<3>(n,a,b,c,ldc);break;case 4:leaf_impl<4>(n,a,b,c,ldc);break;case 5:leaf_impl<5>(n,a,b,c,ldc);break;case 6:leaf_impl<6>(n,a,b,c,ldc);break;case 7:leaf_impl<7>(n,a,b,c,ldc);break;}}
static void sum(int h,double *dst,Mat a,Mat b,int sign) {
for(int i=0;i<h;i++)for(int j=0;j<h;j+=8) {
__m256d v=sign>0?_mm256_add_pd(ld(a,i,j),ld(b,i,j)):_mm256_sub_pd(ld(a,i,j),ld(b,i,j));
__m256d w=sign>0?_mm256_add_pd(ld(a,i,j+4),ld(b,i,j+4)):_mm256_sub_pd(ld(a,i,j+4),ld(b,i,j+4));
asm volatile("" : "+x"(v), "+x"(w) : : "memory");
if(h>=SUMNT) {
_mm256_stream_pd(dst+(size_t)i*h+j,v);_mm256_stream_pd(dst+(size_t)i*h+j+4,w);
} else {
_mm256_store_pd(dst+(size_t)i*h+j,v);_mm256_store_pd(dst+(size_t)i*h+j+4,w);
}
}
if(h>=SUMNT)_mm_sfence();
}
static void combine(int h,double *c,int ldc,const double *m,bool wino=false) {
size_t z=(size_t)h*h+(wino?UPPER_PAGE_PAD:0);
for(int i=0;i<h;i++)for(int j=0;j<h;j+=4) {
size_t t=(size_t)i*h+j;
double *c11=c+(size_t)i*ldc+j,*c12=c11+h,*c21=c+(size_t)(i+h)*ldc+j,*c22=c21+h;
__m256d m1=_mm256_load_pd(c11),m2=_mm256_load_pd(c21),m3=_mm256_load_pd(c12),m4=_mm256_load_pd(c22),
m5=_mm256_load_pd(m+t),m6=_mm256_load_pd(m+z+t),m7=_mm256_load_pd(m+2*z+t);
__m256d r11,r12,r21,r22;
if(wino) {
__m256d u2=_mm256_add_pd(m1,m6),u3=_mm256_add_pd(u2,m7),u4=_mm256_add_pd(u2,m5);
r11=_mm256_add_pd(m1,m2);r12=_mm256_add_pd(u4,m3);r21=_mm256_sub_pd(u3,m4);r22=_mm256_add_pd(u3,m5);
} else {
r11=_mm256_add_pd(_mm256_sub_pd(_mm256_add_pd(m1,m4),m5),m7);r12=_mm256_add_pd(m3,m5);r21=_mm256_add_pd(m2,m4);
r22=_mm256_add_pd(_mm256_add_pd(_mm256_sub_pd(m1,m2),m3),m6);
}
_mm256_store_pd(c11,r11);_mm256_store_pd(c12,r12);_mm256_store_pd(c21,r21);_mm256_store_pd(c22,r22);
}
}
static inline __m256d raw_pre_ld(Mat a,int i,int j){return _mm256_load_pd(a.x+(size_t)i*a.ld+j);}
static void pre(int h,Mat a,double *s1,double *s2,double *s4,double *s3,int ld3,bool right) {
for(int i=0;i<h;i++)for(int j=0;j<h;j+=8) {
__m256d a0=raw_pre_ld(a,i,j),a1=raw_pre_ld(a,i,j+4),b0=raw_pre_ld(a,i,j+h),b1=raw_pre_ld(a,i,j+h+4),
c0=raw_pre_ld(a,i+h,j),c1=raw_pre_ld(a,i+h,j+4),d0=raw_pre_ld(a,i+h,j+h),d1=raw_pre_ld(a,i+h,j+h+4);
__m256d p0,p1,q0,q1,r0,r1,t0,t1;
if(!right) {
p0=_mm256_add_pd(c0,d0);p1=_mm256_add_pd(c1,d1);
q0=_mm256_sub_pd(p0,a0);q1=_mm256_sub_pd(p1,a1);
r0=_mm256_sub_pd(b0,q0);r1=_mm256_sub_pd(b1,q1);
t0=_mm256_sub_pd(a0,c0);t1=_mm256_sub_pd(a1,c1);
}else {
p0=_mm256_sub_pd(b0,a0);p1=_mm256_sub_pd(b1,a1);
q0=_mm256_sub_pd(d0,p0);q1=_mm256_sub_pd(d1,p1);
r0=_mm256_sub_pd(q0,c0);r1=_mm256_sub_pd(q1,c1);
t0=_mm256_sub_pd(d0,b0);t1=_mm256_sub_pd(d1,b1);
}
size_t off=(size_t)i*h+j,off3=(size_t)i*ld3+j;
_mm256_stream_pd(s1+off,p0);_mm256_stream_pd(s1+off+4,p1);
_mm256_stream_pd(s2+off,q0);_mm256_stream_pd(s2+off+4,q1);
_mm256_stream_pd(s4+off,r0);_mm256_stream_pd(s4+off+4,r1);
_mm256_stream_pd(s3+off3,t0);_mm256_stream_pd(s3+off3+4,t1);
}
_mm_sfence();
}
template<int N> static void rec(int rec_n,Mat a,Mat b,double *c,int ldc,double *work,bool root_last=false) {
constexpr int n=N;
if constexpr(n<=BASE){leaf(n,a,b,c,ldc);return;} else {
constexpr int h=n/2;size_t z=(size_t)h*h;
Mat aa[4]={{a.x,a.ld},{a.x+h,a.ld},{a.x+(size_t)h*a.ld,a.ld},{a.x+(size_t)h*a.ld+h,a.ld}};
Mat bb[4]={{b.x,b.ld},{b.x+h,b.ld},{b.x+(size_t)h*b.ld,b.ld},{b.x+(size_t)h*b.ld+h,b.ld}};
double *m=n==2*BASE?work+LEAF_PAD:work;
if(n==2*BASE && pending.active && pending.m==m)m=pending_work+LEAF_PAD;
size_t pz=z+(n==2*BASE?LEAF_PAD:UPPER_PAGE_PAD);
double *p=work+3*z,*q=p+z,*sub=q+z;
double *out[7]={c,c+(size_t)h*ldc,c+h,c+(size_t)h*ldc+h,m,m+pz,m+2*pz};
if(n>BASE*2) {
pending_finish();
double *pa0=m+3*pz,*pa1=pa0+pz,*pa2=pa1+pz,*pb0=pa2+pz,*pb1=pb0+pz,*pb2=pb1+pz;
sub=pb2+pz;
int d4ah_n0=d4ah_ptn;
d4ah_pre_arm(h,a,pa0,pa1,pa2,c,ldc,false);
d4ah_pre_arm(h,b,pb0,pb1,pb2,c+h,ldc,true);
/* mdA_dp2: child order [1,3,4,5,6,0,2] -- constraint "6 before 0 and 2" is the only
hard one; this order puts the three m-panel writers (4,5,6) as late as possible,
which IS the deferral window for the previous sibling's deferred h=1024 combine. */
rec<h>(h,aa[1],bb[2],out[1],ldc,sub);
d4ah_pre_fin_to(d4ah_n0+4);
rec<h>(h,aa[3],{pb2,h},out[3],ldc,sub);
if(root_last) {
pend2_finish(); // M1 of the root is now fully reconstructed.
root_pending.active=true;
root_prefetch();
}
#if !B2M_NO_BARA
pend2_finish(); /* barrier: next call writes m, the deferred products */
#endif
d4ah_pre_fin_to(d4ah_n0+3);
rec<h>(h,{pa0,h},{pb0,h},out[4],h,sub);
d4ah_pre_fin_to(d4ah_n0+2);
rec<h>(h,{pa1,h},{pb1,h},out[5],h,sub);
d4ah_pre_fin_to(d4ah_n0+1);
rec<h>(h,{c,ldc},{c+h,ldc},out[6],h,sub);
rec<h>(h,aa[0],bb[0],out[0],ldc,sub);
if(n==4096) {
root_pending={c,ldc,m,0,false};
d4ah_pre_fin_to(d4ah_n0);
rec<h>(h,{pa2,h},bb[3],out[2],ldc,sub,true);
pending_finish();root_finish();pend2_finish();
return;
}
d4ah_pre_fin_to(d4ah_n0);
rec<h>(h,{pa2,h},bb[3],out[2],ldc,sub);
pending_finish();
pend2_finish(); /* finish an older deferred one */
if(root_last)root_finish();
pend2={c,ldc,m,h,(int)pz,0,true}; /* defer THIS combine (wino form);
m-plane stride = the engine's pz */
if(root_last)pend2.add=root_pending.m;
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};
for(int r=0;r<7;r++) {
Mat x=aa[ai[r]],y=bb[bi[r]];
if(h<=BASE) {
/* Independent product order; M7 is last so reconstruction may
advance behind the produced-row frontier of that leaf. */
const int q=lperm[r];
g_fuse_sd_cur=g_fuse_sd[r];g_fuse_sgn_cur=g_fuse_sgn[r];
Mat x=aa[ai[q]],y=bb[bi[q]];
if(as[q]){x.y=aa[aj[q]].x;x.ly=a.ld;x.sign=as[q];}
if(bs[q]){y.y=bb[bj[q]].x;y.ly=b.ld;y.sign=bs[q];}
if(r==6) {
pending_finish();
pending={c,ldc,m,0,true,0,true};
leaf(h,x,y,out[q],q<4?ldc:h);
pending.prog=false;
} else leaf(h,x,y,out[q],q<4?ldc:h);
} else {
if(as[r]){sum(h,p,x,aa[aj[r]],as[r]);x={p,h};}
if(bs[r]){sum(h,q,y,bb[bj[r]],bs[r]);y={q,h};}
rec<h>(h,x,y,out[r],r<4?ldc:h,sub);
}
}
if(h==BASE) {
/* Keep any final portion deferred into the next sibling. */
} else {
pending_finish();
combine(h,c,ldc,m);
}
}
}
void matrix_multiply(int n,const double *a,const double *b,double *c) {
d1a2_fuse_tab();
madvise((void*)workspace,sizeof(workspace),MADV_HUGEPAGE);
madvise((void*)pending_work,sizeof(pending_work),MADV_HUGEPAGE);
rec<4096>(n,{a,n},{b,n},c,n,workspace);
d4ah_pre_fin_to(0);
pending_finish();
pend2_finish();
_mm_sfence();
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 1.869 s | 502 MB + 336 KB | Accepted | Score: 100 | 显示更多 |