// References:
// [1] saffah_codex_6s_agg2, https://duck.ac/submission/110281:
// Directly reused our accepted integer Winograd/AVX2 source and all
// inherited citations. No separate license notice appeared.
// Approach:
// Form packed B low and rounded-high limb pairs with 16-bit lane blends.
// Each blend selects an original first-row low limb and the second-row
// limb in its high half, replacing masks, shifts, and ORs with fewer AVX2
// instructions while keeping the exact packed bytes.
// Purpose:
// Experimental official timing of the B leaf-packing instruction path.
// References:
// [1] saffah_codex_6s_agg2, https://duck.ac/submission/110211: directly
// reused our accepted 4k Winograd/AVX2 engine and citations below.
// [2] saffah_codex_6s_agg2, https://duck.ac/submission/110277: adapted our
// officially faster 2k top-level row-interleaved result layout. No
// separate license notice was attached to either submission.
// Approach:
// At the top 4k Winograd node, place seven 2k-square result rows adjacent,
// with each child writing at a 7h row stride. This lets the final combine
// read the seven products from nearby addresses. Other levels keep their
// original layout. Modular arithmetic and multiplication order stay exact.
// Purpose:
// Experimental official timing of the 2k layout improvement at 4k scale.
// References:
// [1] saffah_codex_6s_agg2, https://duck.ac/submission/110201: copied our accepted combined integer Winograd and AVX2 engine from mmmi2k, retaining all inherited citations below; no separate license notice.
// [2] saffah_cc_v41_agg1, https://duck.ac/submission/105540: compared its public 4k matrix implementation; this candidate does not copy that source, and no separate license notice appeared.
// Approach: Reuse the 768-int packed-A stride, paired Strassen operand pass, and P2-first Winograd order at n=4096. The matrix recursion and ring arithmetic are generic, with no scale-specific behavior change.
// Purpose: Officially measure whether the accepted 2k optimization reduces 4k runtime below the strict target.
// References:
// [1] saffah_codex_6s_agg2, https://duck.ac/submission/106302:
// Reused our accepted Winograd/AVX2 engine and inherited citations.
// Approach:
// Set the packed-A row stride to 768 int32 values, updating all inline
// assembly offsets consistently. Packed values and ring products are identical.
// Purpose:
// Experimental cache-set and alignment sweep of the packed-A microkernel.
// References:
// [1] Our account saffah_codex_6s_agg2, https://duck.ac/submission/106270: copied the accepted integer Winograd/Strassen engine and all inherited source citations below.
// [2] duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/105523: compared its published mmmi2k B-panel layout; no separate license notice accompanied that public source.
// Approach:
// Port our 1k engine unchanged to the 2k case. It retains the packed-A branch hoist, 24-int packed-B gap, 6-2-4-1-3-5-7 bottom product order, 768-byte packed-B prefetch, and 32-byte kernel loop alignment. The generic recursion accepts n=2048 with exact modulo-2^32 arithmetic.
// Purpose:
// Experimental official submission to price the combined 1k improvements on mmmi2k.
// References:
// [1] Our account saffah_codex_6s_agg2, https://duck.ac/submission/106266: copied our accepted integer microkernel with packed-B prefetch and retained all inherited citations below.
// Approach:
// Adjust the packed-B prefetch loop with al5; the same input data, modular products and output writes are preserved.
// Purpose:
// Experimental official timing of integer microkernel prefetch scheduling on mmmi1k.
// References:
// [1] Our account saffah_codex_6s_agg2, https://duck.ac/submission/106255: copied our accepted bottom-product schedule and its inherited source citations.
// Approach:
// Prefetch the packed B panel 768 bytes ahead once per four-K-pair AVX2 kernel iteration. The packed values, modular arithmetic and stores remain unchanged.
// Purpose:
// Experimental official timing of B-panel cache lead distance on mmmi1k.
// References:
// [1] Our account saffah_codex_6s_agg2, https://duck.ac/submission/106244: copied our accepted integer Winograd engine with all inherited citations retained.
// [2] Our account saffah_codex_6s_agg2, https://duck.ac/submission/103088: prior matrix-product ordering experiment motivated testing adjacent input-block reuse.
// Approach:
// Visit the seven independent bottom Strassen products in the order 6241357. Every product writes its own scratch plane; the output combination still uses the original product operands and integer arithmetic.
// Purpose:
// Experimental official timing of input-block reuse at the 256-square recursion level on mmmi1k.
// References:
// [1] Our account saffah_codex_6s_agg2, https://duck.ac/submission/105115: copied our accepted integer matrix engine, retaining its accumulated source citations below.
// [2] duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/105523: adapted the B-panel gap sweep and cache-set rationale; no separate license notice accompanied that public source.
// Approach:
// Set the B-panel gap to 24 int32 values per eight-column panel, while preserving every packed value and matrix arithmetic.
// Purpose:
// Experimental official timing of B-panel cache placement on mmmi1k.
// References:
// [1] saffah_cc_v41_agg1, https://duck.ac/submission/104850.
// Full public packA-hoisted integer Strassen engine is copied, with its
// upstream attributions retained. No separate license notice appeared.
// [2] saffah_codex_6s_agg2, https://duck.ac/submission/105088.
// Its packed-A stride adjustment from 544 to 528 motivated this sweep.
// Approach:
// Set the packed-A row stride to 512 int32 values and adjust all three
// inline assembly row-offset families together. The 512-value row payload,
// aligned panel access, operation order, and outputs remain identical.
// Purpose:
// Experimental official mmmi1k cache layout and stride measurement.
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #103173 <https://duck.ac/submission/103173>
// 用途：本题三档（mmmi1k/2k/4k）最好件的正文来源（与 mmmi2k 的 #103174、mmmi4k 的 #103175 逐字节相同）。
// [2] 本账号 saffah_cc_v41_agg1，提交 #104472 <https://duck.ac/submission/104472>
// 用途：本件的**直接基座**（`work/mx_pam.cpp`：packA2 行主序 + 死标量尾 + 死函数清理，三档在榜最好件）。
// ======================
// ===== 思路 =====
// 单变量：**把 `packA2` 里的两个"逐行"判定（`if (src2)`、`sg>0 ? add : sub`）提到函数顶层**，
// 让三条内层循环各自**无分支、无数组**（原来 `srcs[4]/src2s[4]` 数组逼得 gcc 把 8 个源指针
// 打包进 ymm 再**溢出到栈上、每个 u 迭代重新加载 8 条**：判题端 gcc 9.3 的 `.isra.0` 克隆里
// 每个 u 迭代 = 8 条指针重载 + 3 条常量向量重载（.LC0/1/2）+ 4 组 `test/je` + 8 条
// `vmovdqa` 寄存器互拷）。
// * 现在（判题端 gcc 9.3 复验，`tools/gcc93check.py`）：packA2 的 2 源 ADD 循环体 =
// `vmovdqa [rcx+rax] / vpaddd ymm1,ymm5,[r9+rax] / vpand / vpaddd / vpackusdw / vpsrld /
// vpermd / vpackusdw / vpermd` ×4 行 + 16 条 Apk 存 + `add rax,32; cmp; jne`——
// **全函数 stackrefs = 0、循环内 test/cmp = 0、指针全在寄存器（rcx/r9/rsi/rdi/r10/r11/rbx/rdx）**。
// * 判据（判题机 in-situ 双腕 A/B，同一二进制、轮流跑、取 min、两臂逐位对拍 bad=0）：
// **n=512 −0.275%、n=1024 −0.363%**（两个规模同号）。
// * 正确性：本地双引擎逐位闸门 n=128/256/512/1024/2048 **C 全元素 FNV 相同**（同一份
// `mx_ver_drv.inc` 驱动两代引擎），且 A/B 二进制里两臂 `bad=0`。
// * 语义：三条循环写的字节、写的次序、算的表达式都与原版**逐位相同**（只是把判定搬家）。
// * 同一轮试过的同类改动（**判题机中性或负，未采纳**）：`packB2` 同样解分支（+0.11%/+0.02%）、
// `tiled2_128` 里把 `first` 固化成常量 1（+0.07%/−0.03%）、`packB2` 源行前瞻预取（+0.40%/+0.55%）。
// ================
#include <immintrin.h>
#include <string.h>
#include <stdint.h>
#include <stdlib.h>
#define MR 4
#define NCG 1
#define KC 128
#define MC 4
#define NC 128
#define NPAIR (KC / 2)
#define APSTR 768 /* 2048 B per row = 16 units x 128 B */
#define NR (8 * NCG)
#define BASE 128
typedef int32_t T;
#define AVX2 __attribute__((target("avx2")))
#define LD(p) _mm256_load_si256((const __m256i *)(p))
#define LDA(p) _mm256_load_si256((const __m256i *)(p))
#define STA(p, v) _mm256_store_si256((__m256i *)(p), (v))
#define ST(p, v) _mm256_store_si256((__m256i *)(p), (v))
static int32_t Apk[(MC + 4) * APSTR + 64] __attribute__((aligned(64)));
static int32_t Bpk[16896] __attribute__((aligned(64))); /* [jg][m][B0|B1] */
/* MIKRO-BEGIN */
#define AVX2 __attribute__((target("avx2")))
AVX2 static inline __attribute__((always_inline)) void mikro(int npair, const int32_t *ap, const int32_t *bp,
int32_t *c, int ldc, int first) {
(void)npair;
__m256i q0, q1, q2, q3;
__asm__ volatile(
"vmovdqu 0(%[bp]), %%ymm8\n\t" "vmovdqu 128(%[bp]), %%ymm9\n\t" "vpmaddwd 0(%[ap]), %%ymm8, %%ymm0\n\t" "vpmaddwd 3072(%[ap]), %%ymm8, %%ymm1\n\t" "vpmaddwd 6144(%[ap]), %%ymm8, %%ymm2\n\t" "vpmaddwd 9216(%[ap]), %%ymm8, %%ymm3\n\t" "vpmaddwd 64(%[ap]), %%ymm8, %%ymm4\n\t" "vpmaddwd 3136(%[ap]), %%ymm8, %%ymm5\n\t" "vpmaddwd 6208(%[ap]), %%ymm8, %%ymm6\n\t" "vpmaddwd 9280(%[ap]), %%ymm8, %%ymm7\n\t" "vpmaddwd 0(%[ap]), %%ymm9, %%ymm10\n\t" "vpmaddwd 3072(%[ap]), %%ymm9, %%ymm11\n\t" "vpmaddwd 6144(%[ap]), %%ymm9, %%ymm12\n\t" "vpmaddwd 9216(%[ap]), %%ymm9, %%ymm13\n\t" "vpaddd %%ymm10, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm11, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm12, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm13, %%ymm7, %%ymm7\n\t" "vmovdqu 32(%[bp]), %%ymm8\n\t" "vmovdqu 160(%[bp]), %%ymm9\n\t" "vpmaddwd 4(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 3076(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 6148(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 9220(%[ap]), %%ymm8, %%ymm13\n\t" "vpaddd %%ymm10, %%ymm0, %%ymm0\n\t" "vpaddd %%ymm11, %%ymm1, %%ymm1\n\t" "vpaddd %%ymm12, %%ymm2, %%ymm2\n\t" "vpaddd %%ymm13, %%ymm3, %%ymm3\n\t" "vpmaddwd 68(%[ap]), %%ymm8, %%ymm14\n\t" "vpmaddwd 3140(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 6212(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 9284(%[ap]), %%ymm8, %%ymm11\n\t" "vpaddd %%ymm14, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm15, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm10, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm11, %%ymm7, %%ymm7\n\t" "vpmaddwd 4(%[ap]), %%ymm9, %%ymm12\n\t" "vpmaddwd 3076(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 6148(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 9220(%[ap]), %%ymm9, %%ymm15\n\t" "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm14, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t" "vmovdqu 64(%[bp]), %%ymm8\n\t" "vmovdqu 192(%[bp]), %%ymm9\n\t" "vpmaddwd 8(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 3080(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 6152(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 9224(%[ap]), %%ymm8, %%ymm13\n\t" "vpaddd %%ymm10, %%ymm0, %%ymm0\n\t" "vpaddd %%ymm11, %%ymm1, %%ymm1\n\t" "vpaddd %%ymm12, %%ymm2, %%ymm2\n\t" "vpaddd %%ymm13, %%ymm3, %%ymm3\n\t" "vpmaddwd 72(%[ap]), %%ymm8, %%ymm14\n\t" "vpmaddwd 3144(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 6216(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 9288(%[ap]), %%ymm8, %%ymm11\n\t" "vpaddd %%ymm14, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm15, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm10, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm11, %%ymm7, %%ymm7\n\t" "vpmaddwd 8(%[ap]), %%ymm9, %%ymm12\n\t" "vpmaddwd 3080(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 6152(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 9224(%[ap]), %%ymm9, %%ymm15\n\t" "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm14, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t" "vmovdqu 96(%[bp]), %%ymm8\n\t" "vmovdqu 224(%[bp]), %%ymm9\n\t" "vpmaddwd 12(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 3084(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 6156(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 9228(%[ap]), %%ymm8, %%ymm13\n\t" "vpaddd %%ymm10, %%ymm0, %%ymm0\n\t" "vpaddd %%ymm11, %%ymm1, %%ymm1\n\t" "vpaddd %%ymm12, %%ymm2, %%ymm2\n\t" "vpaddd %%ymm13, %%ymm3, %%ymm3\n\t" "vpmaddwd 76(%[ap]), %%ymm8, %%ymm14\n\t" "vpmaddwd 3148(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 6220(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 9292(%[ap]), %%ymm8, %%ymm11\n\t" "vpaddd %%ymm14, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm15, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm10, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm11, %%ymm7, %%ymm7\n\t" "vpmaddwd 12(%[ap]), %%ymm9, %%ymm12\n\t" "vpmaddwd 3084(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 6156(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 9228(%[ap]), %%ymm9, %%ymm15\n\t" "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm14, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t" "add $128, %[ap]\n\t"
"add $256, %[bp]\n\t"
"sub $4, %[cnt]\n\t"
"jle 2f\n\t"
".p2align 5\n\t"
"1:\n\t" "prefetcht0 768(%[bp])\n\t"
"vmovdqu 0(%[bp]), %%ymm8\n\t" "vmovdqu 128(%[bp]), %%ymm9\n\t" "vpmaddwd 0(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 3072(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 6144(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 9216(%[ap]), %%ymm8, %%ymm13\n\t" "vpaddd %%ymm10, %%ymm0, %%ymm0\n\t" "vpaddd %%ymm11, %%ymm1, %%ymm1\n\t" "vpaddd %%ymm12, %%ymm2, %%ymm2\n\t" "vpaddd %%ymm13, %%ymm3, %%ymm3\n\t" "vpmaddwd 64(%[ap]), %%ymm8, %%ymm14\n\t" "vpmaddwd 3136(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 6208(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 9280(%[ap]), %%ymm8, %%ymm11\n\t" "vpaddd %%ymm14, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm15, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm10, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm11, %%ymm7, %%ymm7\n\t" "vpmaddwd 0(%[ap]), %%ymm9, %%ymm12\n\t" "vpmaddwd 3072(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 6144(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 9216(%[ap]), %%ymm9, %%ymm15\n\t" "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm14, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t" "vmovdqu 32(%[bp]), %%ymm8\n\t" "vmovdqu 160(%[bp]), %%ymm9\n\t" "vpmaddwd 4(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 3076(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 6148(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 9220(%[ap]), %%ymm8, %%ymm13\n\t" "vpaddd %%ymm10, %%ymm0, %%ymm0\n\t" "vpaddd %%ymm11, %%ymm1, %%ymm1\n\t" "vpaddd %%ymm12, %%ymm2, %%ymm2\n\t" "vpaddd %%ymm13, %%ymm3, %%ymm3\n\t" "vpmaddwd 68(%[ap]), %%ymm8, %%ymm14\n\t" "vpmaddwd 3140(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 6212(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 9284(%[ap]), %%ymm8, %%ymm11\n\t" "vpaddd %%ymm14, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm15, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm10, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm11, %%ymm7, %%ymm7\n\t" "vpmaddwd 4(%[ap]), %%ymm9, %%ymm12\n\t" "vpmaddwd 3076(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 6148(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 9220(%[ap]), %%ymm9, %%ymm15\n\t" "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm14, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t" "vmovdqu 64(%[bp]), %%ymm8\n\t" "vmovdqu 192(%[bp]), %%ymm9\n\t" "vpmaddwd 8(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 3080(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 6152(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 9224(%[ap]), %%ymm8, %%ymm13\n\t" "vpaddd %%ymm10, %%ymm0, %%ymm0\n\t" "vpaddd %%ymm11, %%ymm1, %%ymm1\n\t" "vpaddd %%ymm12, %%ymm2, %%ymm2\n\t" "vpaddd %%ymm13, %%ymm3, %%ymm3\n\t" "vpmaddwd 72(%[ap]), %%ymm8, %%ymm14\n\t" "vpmaddwd 3144(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 6216(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 9288(%[ap]), %%ymm8, %%ymm11\n\t" "vpaddd %%ymm14, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm15, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm10, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm11, %%ymm7, %%ymm7\n\t" "vpmaddwd 8(%[ap]), %%ymm9, %%ymm12\n\t" "vpmaddwd 3080(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 6152(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 9224(%[ap]), %%ymm9, %%ymm15\n\t" "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm14, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t" "vmovdqu 96(%[bp]), %%ymm8\n\t" "vmovdqu 224(%[bp]), %%ymm9\n\t" "vpmaddwd 12(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 3084(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 6156(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 9228(%[ap]), %%ymm8, %%ymm13\n\t" "vpaddd %%ymm10, %%ymm0, %%ymm0\n\t" "vpaddd %%ymm11, %%ymm1, %%ymm1\n\t" "vpaddd %%ymm12, %%ymm2, %%ymm2\n\t" "vpaddd %%ymm13, %%ymm3, %%ymm3\n\t" "vpmaddwd 76(%[ap]), %%ymm8, %%ymm14\n\t" "vpmaddwd 3148(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 6220(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 9292(%[ap]), %%ymm8, %%ymm11\n\t" "vpaddd %%ymm14, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm15, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm10, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm11, %%ymm7, %%ymm7\n\t" "vpmaddwd 12(%[ap]), %%ymm9, %%ymm12\n\t" "vpmaddwd 3084(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 6156(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 9228(%[ap]), %%ymm9, %%ymm15\n\t" "vpaddd %%ymm12, %%ymm4, %%ymm4\n\t" "vpaddd %%ymm13, %%ymm5, %%ymm5\n\t" "vpaddd %%ymm14, %%ymm6, %%ymm6\n\t" "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t" "add $128, %[ap]\n\t"
"add $256, %[bp]\n\t"
"sub $4, %[cnt]\n\t"
"jnz 1b\n\t"
"2:\n\t"
"vpslld $16, %%ymm4, %%ymm4\n\t"
"vpaddd %%ymm4, %%ymm0, %[q0]\n\t"
"vpslld $16, %%ymm5, %%ymm5\n\t"
"vpaddd %%ymm5, %%ymm1, %[q1]\n\t"
"vpslld $16, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm6, %%ymm2, %[q2]\n\t"
"vpslld $16, %%ymm7, %%ymm7\n\t"
"vpaddd %%ymm7, %%ymm3, %[q3]\n\t"
: [ap] "+r" (ap), [bp] "+r" (bp), [cnt] "+r" (npair),
[q0] "=x" (q0), [q1] "=x" (q1), [q2] "=x" (q2), [q3] "=x" (q3)
:
: "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7", "ymm8", "ymm9", "ymm10", "ymm15", "memory"); /* ymm15 kept DECLARED on purpose -- see notes */
{ int32_t *cp = c;
if (first) {
_mm256_store_si256((__m256i *)(cp + (size_t)0 * ldc), q0);
_mm256_store_si256((__m256i *)(cp + (size_t)1 * ldc), q1);
_mm256_store_si256((__m256i *)(cp + (size_t)2 * ldc), q2);
_mm256_store_si256((__m256i *)(cp + (size_t)3 * ldc), q3);
} else {
int32_t *cpr0 = cp + (size_t)0 * ldc;
_mm256_store_si256((__m256i *)cpr0, _mm256_add_epi32(_mm256_load_si256((const __m256i *)cpr0), q0));
int32_t *cpr1 = cp + (size_t)1 * ldc;
_mm256_store_si256((__m256i *)cpr1, _mm256_add_epi32(_mm256_load_si256((const __m256i *)cpr1), q1));
int32_t *cpr2 = cp + (size_t)2 * ldc;
_mm256_store_si256((__m256i *)cpr2, _mm256_add_epi32(_mm256_load_si256((const __m256i *)cpr2), q2));
int32_t *cpr3 = cp + (size_t)3 * ldc;
_mm256_store_si256((__m256i *)cpr3, _mm256_add_epi32(_mm256_load_si256((const __m256i *)cpr3), q3));
} }
}
/* MIKRO-END */
/* ---- ld-aware packing + tiled macro-kernel (base case) ---- */
#pragma GCC push_options
#pragma GCC optimize("O3")
AVX2 static inline void packA2(int ic, int kk, int n, int lda, const T * __restrict A, int sg, const T * __restrict A2) {
const __m256i m16 = _mm256_set1_epi32(0xFFFF);
const __m256i half = _mm256_set1_epi32(0x8000);
const __m256i ridx = _mm256_setr_epi32(0, 1, 4, 5, 0, 1, 4, 5);
(void)n;
/* UNSWITCHED: the (A2!=0) and (sg>0) decisions are hoisted to the top so each
inner loop is branch-free and no source pointer has to live in an array.
Same stores, same order of arithmetic -> bit-identical panels. */
const size_t st = (size_t)lda;
const T *s0 = A + (size_t)ic * st + kk, *s1 = s0 + st, *s2 = s1 + st, *s3 = s2 + st;
#define PA_EMIT(RO, XV) \
do { __m256i al = _mm256_and_si256((XV), m16); \
__m256i ah = _mm256_srli_epi32(_mm256_add_epi32((XV), half), 16); \
__m256i rl = _mm256_permutevar8x32_epi32(_mm256_packus_epi32(al, al), ridx); \
__m256i rh = _mm256_permutevar8x32_epi32(_mm256_packus_epi32(ah, ah), ridx); \
T *d = Apk + (size_t)(RO) * APSTR + 32 * u; \
ST(d + 0, rl); _mm_store_si128((__m128i *)(d + 8), _mm256_castsi256_si128(rl)); \
ST(d + 16, rh); _mm_store_si128((__m128i *)(d + 24), _mm256_castsi256_si128(rh)); } while (0)
if (A2) {
const T *t0 = A2 + (size_t)ic * st + kk, *t1 = t0 + st, *t2 = t1 + st, *t3 = t2 + st;
if (sg > 0) {
for (int u = 0; u < KC / 8; u++) {
PA_EMIT(0, _mm256_add_epi32(LD(s0 + 8 * u), LD(t0 + 8 * u)));
PA_EMIT(1, _mm256_add_epi32(LD(s1 + 8 * u), LD(t1 + 8 * u)));
PA_EMIT(2, _mm256_add_epi32(LD(s2 + 8 * u), LD(t2 + 8 * u)));
PA_EMIT(3, _mm256_add_epi32(LD(s3 + 8 * u), LD(t3 + 8 * u)));
}
} else {
for (int u = 0; u < KC / 8; u++) {
PA_EMIT(0, _mm256_sub_epi32(LD(s0 + 8 * u), LD(t0 + 8 * u)));
PA_EMIT(1, _mm256_sub_epi32(LD(s1 + 8 * u), LD(t1 + 8 * u)));
PA_EMIT(2, _mm256_sub_epi32(LD(s2 + 8 * u), LD(t2 + 8 * u)));
PA_EMIT(3, _mm256_sub_epi32(LD(s3 + 8 * u), LD(t3 + 8 * u)));
}
}
} else {
for (int u = 0; u < KC / 8; u++) {
PA_EMIT(0, LD(s0 + 8 * u)); PA_EMIT(1, LD(s1 + 8 * u));
PA_EMIT(2, LD(s2 + 8 * u)); PA_EMIT(3, LD(s3 + 8 * u));
}
}
#undef PA_EMIT
}
#pragma GCC pop_options
#pragma GCC push_options
#pragma GCC optimize("O3")
AVX2 static inline void packB2(int jc, int kk, int n, int ldb, const T * __restrict B, int sg, const T * __restrict B2) {
const __m256i half = _mm256_set1_epi32(0x8000);
(void)n;
for (int u = 0; u < NPAIR / 4; u++) {
const T *s = B + (size_t)(kk + 8 * u) * ldb + jc;
const T *t = B2 ? (B2 + (size_t)(kk + 8 * u) * ldb + jc) : 0;
for (int j = 0; j < NC; j += 8) {
T *d = Bpk + (size_t)(j / 8) * (NPAIR * 16 + 24) + (size_t)u * 64;
__m256i L0, L1, L2, L3, H0, H1, H2, H3;
#define F11_PAIR(P, LO, HI) \
{ const T *sp = s + (size_t)(2 * P) * ldb + j; \
__m256i x0 = LD(sp), x1 = LD(sp + ldb); \
if (t) { const T *tp = t + (size_t)(2 * P) * ldb + j; \
__m256i y0 = LD(tp), y1 = LD(tp + ldb); \
if (sg > 0) { x0 = _mm256_add_epi32(x0, y0); x1 = _mm256_add_epi32(x1, y1); } \
else { x0 = _mm256_sub_epi32(x0, y0); x1 = _mm256_sub_epi32(x1, y1); } } \
LO = _mm256_blend_epi16(x0, _mm256_slli_epi32(x1, 16), 0xaa);\
HI = _mm256_blend_epi16(_mm256_srli_epi32(_mm256_add_epi32(x0, half), 16),\
_mm256_add_epi32(x1, half), 0xaa); }
F11_PAIR(0, L0, H0)
F11_PAIR(1, L1, H1)
F11_PAIR(2, L2, H2)
F11_PAIR(3, L3, H3)
#undef F11_PAIR
ST(d + 0, _mm256_blend_epi32(_mm256_blend_epi32(L0, L2, 0x44), _mm256_blend_epi32(L1, L3, 0x88), 0xaa));
ST(d + 32, _mm256_blend_epi32(_mm256_blend_epi32(H0, H2, 0x44), _mm256_blend_epi32(H1, H3, 0x88), 0xaa));
ST(d + 8, _mm256_blend_epi32(_mm256_blend_epi32(L0, L2, 0x22), _mm256_blend_epi32(L1, L3, 0x44), 0x55));
ST(d + 40, _mm256_blend_epi32(_mm256_blend_epi32(H0, H2, 0x22), _mm256_blend_epi32(H1, H3, 0x44), 0x55));
ST(d + 16, _mm256_blend_epi32(_mm256_blend_epi32(L0, L2, 0x11), _mm256_blend_epi32(L1, L3, 0x22), 0xaa));
ST(d + 48, _mm256_blend_epi32(_mm256_blend_epi32(H0, H2, 0x11), _mm256_blend_epi32(H1, H3, 0x22), 0xaa));
ST(d + 24, _mm256_blend_epi32(_mm256_blend_epi32(L0, L2, 0x88), _mm256_blend_epi32(L1, L3, 0x11), 0x55));
ST(d + 56, _mm256_blend_epi32(_mm256_blend_epi32(H0, H2, 0x88), _mm256_blend_epi32(H1, H3, 0x11), 0x55));
}
}
}
#pragma GCC pop_options
AVX2 static inline __attribute__((always_inline)) void block_mikro(int ic, int jc, int n, int ldc, T *C, int first) {
T *cp0 = C + (size_t)ic * ldc + jc;
for (int r = 0; r < MC; r += MR) {
const int32_t *ap = Apk + (size_t)r * APSTR;
T *cp = cp0 + (size_t)r * ldc;
const int32_t *bp = Bpk;
for (int jg = 0; jg < NC; jg += NR) {
mikro(NPAIR, ap, bp, cp + jg, ldc, first);
bp += (NPAIR * 16 + 24);
}
}
}
/* LEAF SPECIALISED: with BASE==KC==NC==128 every tiled2 leaf call has h==128 and
ldc==h==128, so both become compile-time constants -> the jc/kk/ic loops fully fold. */
AVX2 static void tiled2_128(const T *A1, const T *A2, int sgA, int lda,
const T *B1, const T *B2, int sgB, int ldb, T *C) {
const int n = 128;
const int ldc = 128;
for (int jc = 0; jc < n; jc += NC) {
for (int kk = 0; kk < n; kk += KC) {
packB2(jc, kk, n, ldb, B1, sgB, B2);
for (int ic = 0; ic < n; ic += MC) {
packA2(ic, kk, n, lda, A1, sgA, A2);
if (ic + MC < n) {
const T *n1 = A1 + (size_t)(ic + MC) * lda + kk;
const T *n2 = A2 ? (A2 + (size_t)(ic + MC) * lda + kk) : 0;
for (int r = 0; r < MC; r++) {
const T *p1 = n1 + (size_t)r * lda;
const T *p2 = n2 ? (n2 + (size_t)r * lda) : 0;
__builtin_prefetch(p1, 0, 2);
if (p2) __builtin_prefetch(p2, 0, 2);
}
}
block_mikro(ic, jc, n, ldc, C, kk == 0);
}
}
}
}
/* correct for any n, just slow */
AVX2 static void generic(int n, const T *A, int lda, const T *B, int ldb, T *C, int ldc) {
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++) {
uint32_t s = 0;
for (int k = 0; k < n; k++)
s += (uint32_t)A[(size_t)i * lda + k] * (uint32_t)B[(size_t)k * ldb + j];
C[(size_t)i * ldc + j] = (T)s;
}
}
/* tiled base case: requires n % NC == 0, n % KC == 0, n >= MC */
AVX2 static void tiled(int n, const T *A, int lda, const T *B, int ldb, T *C, int ldc) {
/* C is fully written by the kk==0 pass (micro stores instead of RMWs) */
for (int jc = 0; jc < n; jc += NC) {
for (int kk = 0; kk < n; kk += KC) {
packB2(jc, kk, n, ldb, B, 1, (const T *)0);
for (int ic = 0; ic < n; ic += MC) {
packA2(ic, kk, n, lda, A, 1, (const T *)0);
block_mikro(ic, jc, n, ldc, C, kk == 0);
}
}
}
}
AVX2 static int ok_tiled(int n) {
return (n % NC) == 0 && (n % KC) == 0 && (n % MC) == 0 && (MC % MR) == 0 && (NC % NR) == 0;
}
/* ---- Strassen over Z/2^32. Exact: all ops are ring ops of Z/2^32. ---- */
#define STR_DEPTH 6
/* P = X + s*Y in ONE pass over h x h (replaces copyblk + axpy: one write instead of two) */
#pragma GCC push_options
#pragma GCC optimize("O3")
AVX2 static void addblk(int h, T * __restrict P, const T * __restrict X, int ldx, const T * __restrict Y, int ldy, int s) {
for (int i = 0; i < h; i++) {
const T *xp = X + (size_t)i * ldx, *yp = Y + (size_t)i * ldy;
T *pp = P + (size_t)i * h;
__builtin_prefetch(X + (size_t)(i+1)*ldx, 0, 3); __builtin_prefetch(Y + (size_t)(i+1)*ldy, 0, 3);
int j = 0;
if (s > 0) for (; j + 8 <= h; j += 8) STA(pp + j, _mm256_add_epi32(LDA(xp + j), LDA(yp + j)));
else for (; j + 8 <= h; j += 8) STA(pp + j, _mm256_sub_epi32(LDA(xp + j), LDA(yp + j)));
/* scalar tail removed: h is always a multiple of 8 (h = n/2, n a power of two >= 256);
mutate_blk in this engine already relies on the same invariant. */
}
}
#pragma GCC pop_options
#pragma GCC push_options
#pragma GCC optimize("O3")
template<int SA, int SB>
AVX2 static void addblk_pair(int h, T * __restrict P,
const T * __restrict A, int lda, const T * __restrict A2,
T * __restrict Q,
const T * __restrict B, int ldb, const T * __restrict B2) {
for (int i=0;i<h;i++) {
const T *a=A+(size_t)i*lda,*a2=A2+(size_t)i*lda;
const T *b=B+(size_t)i*ldb,*b2=B2+(size_t)i*ldb;
T *p=P+(size_t)i*h,*q=Q+(size_t)i*h;
__builtin_prefetch(a+lda,0,3); __builtin_prefetch(a2+lda,0,3);
__builtin_prefetch(b+ldb,0,3); __builtin_prefetch(b2+ldb,0,3);
for(int j=0;j<h;j+=8) {
__m256i x=LD(a+j), y=LD(a2+j), z=LD(b+j), w=LD(b2+j);
STA(p+j,SA>0?_mm256_add_epi32(x,y):_mm256_sub_epi32(x,y));
STA(q+j,SB>0?_mm256_add_epi32(z,w):_mm256_sub_epi32(z,w));
}
}
}
#pragma GCC pop_options
/* dst = s1*m1 + s2*m2 , ONE pass (M blocks are contiguous, row stride h) */
/* dst = s1*m1 + s2*m2 + s3*m3 + s4*m4 , ONE pass */
#define B11 (A)
#define B12 (A + h)
#define B21 (A + (size_t)h * lda)
#define B22 (A + (size_t)h * lda + h)
#define D11 (B)
#define D12 (B + h)
#define D21 (B + (size_t)h * ldb)
#define D22 (B + (size_t)h * ldb + h)
/* ---- FUSED combine: one pass instead of four (l_sub1pct) --------------------
The shipped rec() makes FOUR passes over the seven M blocks (12 block-reads
for 4 outputs: M1..M5 are each read twice). Holding the seven row vectors in
registers lets all four quadrants be formed from 7 loads, so the combine's
issue cost per 8 outputs goes 24 -> 19 uops and its read volume 12h^2 -> 7h^2.
Same stores, same NT form, byte-exact. */
#pragma GCC push_options
#pragma GCC optimize("O3")
AVX2 static void comb_all(int h, T * __restrict C, int ldc,
const T * __restrict M1, const T * __restrict M2, const T * __restrict M3, const T * __restrict M4,
const T * __restrict M5, const T * __restrict M6, const T * __restrict M7) {
for (int i = 0; i < h; i++) {
const T *a = M1 + (size_t)i * h, *b = M2 + (size_t)i * h, *c = M3 + (size_t)i * h,
*d = M4 + (size_t)i * h, *e = M5 + (size_t)i * h, *f = M6 + (size_t)i * h,
*g = M7 + (size_t)i * h;
T *d1 = C + (size_t)i * ldc, *d2 = d1 + h,
*d3 = d1 + (size_t)h * ldc, *d4 = d3 + h;
int j = 0;
for (; j + 8 <= h; j += 8) {
__m256i m1 = LDA(a + j), m2 = LDA(b + j), m3 = LDA(c + j), m4 = LDA(d + j),
m5 = LDA(e + j), m6 = LDA(f + j), m7 = LDA(g + j);
_mm256_stream_si256((__m256i *)(d1 + j),
_mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(m1, m4), m5), m7));
_mm256_stream_si256((__m256i *)(d2 + j), _mm256_add_epi32(m3, m5));
_mm256_stream_si256((__m256i *)(d3 + j), _mm256_add_epi32(m2, m4));
_mm256_stream_si256((__m256i *)(d4 + j),
_mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(m1, m3), m2), m6));
}
/* scalar tail removed: same invariant as above. */
}
}
#pragma GCC pop_options
AVX2 static void rec(int n, const T *A, int lda, const T *B, int ldb, T *C, int ldc,
T *scr, int depth) {
if (depth >= STR_DEPTH || (n & 1) || n <= BASE || !ok_tiled(n)) {
if (ok_tiled(n)) tiled(n, A, lda, B, ldb, C, ldc);
else generic(n, A, lda, B, ldb, C, ldc);
return;
}
const int h = n / 2;
const size_t h2 = (size_t)h * h;
T *S = scr; /* 7 result blocks M1..M7 then 2 scratch S blocks */
T *M1 = S, *M2 = M1 + h2 + 4608, *M3 = M2 + h2 + 4608, *M4 = M3 + h2 + 4608,
*M5 = M4 + h2 + 4608, *M6 = M5 + h2 + 4608, *M7 = M6 + h2 + 4608;
T *P = M7 + h2 + 4608; /* operand scratch block 1 (h x h, ldp = h) */
T *Q = P + h2 + 4608; /* operand scratch block 2 */
T *sub = Q + h2 + 4608; /* recursion arena */
/* Each of the 7 products uses an operand that is A1 (+/- A2) on the A side and
B1 (+/- B2) on the B side. When the children are LEAVES (h <= BASE) the summed
operand is consumed only by the packers, so we fuse the sum INTO the packer and
never materialise P/Q: that removes 2 of every 4 accesses per summed operand
element (addblk's write + the packer's read of it, replaced by a second load of
the original block, which is exactly what addblk did anyway). */
#define LOP(a1, a2, sgA, b1, b2, sgB, Mi) \
do { \
if (h <= BASE) { \
tiled2_128((a1), (a2), (sgA), lda, \
(b1), (b2), (sgB), ldb, (Mi)); \
} else { \
const T *pa = (a1); int la = lda; \
const T *pb = (b1); int lb = ldb; \
if ((a2) && (b2)) { \
addblk_pair<(sgA),(sgB)>(h,P,(a1),lda,(a2),Q,(b1),ldb,(b2)); \
pa=P; la=h; pb=Q; lb=h; \
} else { \
if (a2) { addblk(h,P,(a1),lda,(a2),lda,(sgA)); pa=P; la=h; } \
if (b2) { addblk(h,Q,(b1),ldb,(b2),ldb,(sgB)); pb=Q; lb=h; } \
} \
rec(h, pa, la, pb, lb, (Mi), h, sub, depth + 1); \
} \
} while (0)
LOP(B21, B11, -1, D11, D12, +1, M6); /* M6 = (A21-A11)(B11+B12) */
LOP(B21, B22, +1, D11, 0, +1, M2); /* M2 = (A21+A22) B11 */
LOP(B22, 0, +1, D21, D11, -1, M4); /* M4 = A22 (B21-B11) */
LOP(B11, B22, +1, D11, D22, +1, M1); /* M1 = (A11+A22)(B11+B22) */
LOP(B11, 0, +1, D12, D22, -1, M3); /* M3 = A11 (B12-B22) */
LOP(B11, B12, +1, D22, 0, +1, M5); /* M5 = (A11+A12) B22 */
LOP(B12, B22, -1, D21, D22, +1, M7); /* M7 = (A12-A22)(B21+B22) */
#undef LOP
/* C11 = M1 + M4 - M5 + M7 ; C12 = M3 + M5 ; C21 = M2 + M4 ; C22 = M1 - M2 + M3 + M6 */
comb_all(h, C, ldc, M1, M2, M3, M4, M5, M6, M7);
}
#undef B11
#undef B12
#undef B21
#undef B22
#undef D11
#undef D12
#undef D21
#undef D22
/* Update an already materialized S or T block without a separate copy.
[7] verbatim. */
#pragma GCC push_options
#pragma GCC optimize("O3")
AVX2 static void mutate_blk(int h, T *dst, const T *src, int lds, int reverse) {
for (int i=0;i<h;i++) {
T *d=dst+(size_t)i*h; const T *s=src+(size_t)i*lds;
if (reverse) {
for(int j=0;j<h;j+=8) STA(d+j,_mm256_sub_epi32(LDA(s+j),LDA(d+j)));
} else {
for(int j=0;j<h;j+=8) STA(d+j,_mm256_sub_epi32(LDA(d+j),LDA(s+j)));
}
}
}
AVX2 static void comb_win(int h, T *C, int ldc, const T *M1, const T *M2,
const T *M3, const T *M4, const T *M5,
const T *M6, const T *M7) {
for(int i=0;i<h;i++) {
const T *a=M1+(size_t)i*h,*b=M2+(size_t)i*h,*c=M3+(size_t)i*h,
*d=M4+(size_t)i*h,*e=M5+(size_t)i*h,*f=M6+(size_t)i*h,*g=M7+(size_t)i*h;
T *o11=C+(size_t)i*ldc,*o12=o11+h,*o21=o11+(size_t)h*ldc,*o22=o21+h;
for(int j=0;j<h;j+=8) {
__m256i p1=LDA(a+j),p2=LDA(b+j),p3=LDA(c+j),p4=LDA(d+j),
p5=LDA(e+j),p6=LDA(f+j),p7=LDA(g+j);
__m256i u2=_mm256_add_epi32(p1,p6);
__m256i u3=_mm256_add_epi32(u2,p7);
__m256i u4=_mm256_add_epi32(u2,p5);
_mm256_stream_si256((__m256i*)(o11+j),_mm256_add_epi32(p1,p2));
_mm256_stream_si256((__m256i*)(o12+j),_mm256_add_epi32(u4,p3));
_mm256_stream_si256((__m256i*)(o21+j),_mm256_sub_epi32(u3,p4));
_mm256_stream_si256((__m256i*)(o22+j),_mm256_add_epi32(u3,p5));
}
}
}
AVX2 static void comb_win_interleaved(int h, T *C, int ldc, const T *M1, const T *M2,
const T *M3, const T *M4, const T *M5,
const T *M6, const T *M7) {
for(int i=0;i<h;i++) {
const T *a=M1+(size_t)i*(7*h),*b=M2+(size_t)i*(7*h),*c=M3+(size_t)i*(7*h),
*d=M4+(size_t)i*(7*h),*e=M5+(size_t)i*(7*h),*f=M6+(size_t)i*(7*h),*g=M7+(size_t)i*(7*h);
T *o11=C+(size_t)i*ldc,*o12=o11+h,*o21=o11+(size_t)h*ldc,*o22=o21+h;
for(int j=0;j<h;j+=8) {
__m256i p1=LDA(a+j),p2=LDA(b+j),p3=LDA(c+j),p4=LDA(d+j),
p5=LDA(e+j),p6=LDA(f+j),p7=LDA(g+j);
__m256i u2=_mm256_add_epi32(p1,p6);
__m256i u3=_mm256_add_epi32(u2,p7);
__m256i u4=_mm256_add_epi32(u2,p5);
_mm256_stream_si256((__m256i*)(o11+j),_mm256_add_epi32(p1,p2));
_mm256_stream_si256((__m256i*)(o12+j),_mm256_add_epi32(u4,p3));
_mm256_stream_si256((__m256i*)(o21+j),_mm256_sub_epi32(u3,p4));
_mm256_stream_si256((__m256i*)(o22+j),_mm256_add_epi32(u3,p5));
}
}
}
#pragma GCC pop_options
AVX2 static void rec_w(int n, const T *A, int lda, const T *B, int ldb,
T *C, int ldc, T *scr, int depth) {
if (n<=256 || depth>=STR_DEPTH || (n&1) || !ok_tiled(n)) {
rec(n,A,lda,B,ldb,C,ldc,scr,depth); return;
}
const int h=n/2; const size_t h2=(size_t)h*h;
T *M1=scr,*M2=M1+h2+4608,*M3=M2+h2+4608,*M4=M3+h2+4608,
*M5=M4+h2+4608,*M6=M5+h2+4608,*M7=M6+h2+4608;
T *P=M7+h2+4608,*Q=P+h2+4608,*sub=Q+h2+4608;
const bool interleaved = (depth == 0);
if (interleaved) {
M1=scr; M2=M1+h; M3=M2+h; M4=M3+h; M5=M4+h; M6=M5+h; M7=M6+h;
P=scr+7*h2+7*4608; Q=P+h2+4608; sub=Q+h2+4608;
}
const int mld=interleaved?7*h:h;
const T *A11=A,*A12=A+h,*A21=A+(size_t)h*lda,*A22=A21+h;
const T *B11=B,*B12=B+h,*B21=B+(size_t)h*ldb,*B22=B21+h;
rec_w(h,A12,lda,B21,ldb,M2,mld,sub,depth+1); // P2
rec_w(h,A11,lda,B11,ldb,M1,mld,sub,depth+1); // P1
addblk(h,P,A21,lda,A22,lda,+1); // S1
addblk(h,Q,B12,ldb,B11,ldb,-1); // T1
rec_w(h,P,h,Q,h,M5,mld,sub,depth+1); // P5
mutate_blk(h,P,A11,lda,0); // S2 = S1-A11
mutate_blk(h,Q,B22,ldb,1); // T2 = B22-T1
rec_w(h,P,h,Q,h,M6,mld,sub,depth+1); // P6
mutate_blk(h,P,A12,lda,1); // S4 = A12-S2
rec_w(h,P,h,B22,ldb,M3,mld,sub,depth+1); // P3
mutate_blk(h,Q,B21,ldb,0); // T4 = T2-B21
rec_w(h,A22,lda,Q,h,M4,mld,sub,depth+1); // P4
addblk(h,P,A11,lda,A21,lda,-1); // S3
addblk(h,Q,B22,ldb,B12,ldb,-1); // T3
rec_w(h,P,h,Q,h,M7,mld,sub,depth+1); // P7
if(interleaved) comb_win_interleaved(h,C,ldc,M1,M2,M3,M4,M5,M6,M7);
else comb_win(h,C,ldc,M1,M2,M3,M4,M5,M6,M7);
}
AVX2 void matrix_multiply(int n, const T *A, const T *B, T *C) {
unsigned long long al = (unsigned long long)(size_t)A | (unsigned long long)(size_t)B | (unsigned long long)(size_t)C;
if (n <= BASE || !ok_tiled(n) || (al & 31)) {
if (ok_tiled(n) && !(al & 31)) tiled(n, A, n, B, n, C, n);
else generic(n, A, n, B, n, C, n);
return;
}
/* depth-first arena: 9 h^2 blocks at each level along the path */
size_t need = 0, h = (size_t)n;
while (h > (size_t)BASE) { h /= 2; need += 9 * (h * h + 4608); }
static T *scr = 0; static size_t cap = 0;
if (need > cap) { if (scr) free(scr); cap = need;
void *pv = 0; if (posix_memalign(&pv, 4096, cap * sizeof(T))) return; scr = (T *)pv; }
rec_w(n, A, n, B, n, C, n, scr, 0);
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 1.743 s | 255 MB + 876 KB | Accepted | Score: 100 | 显示更多 |