// ===== 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，提交 #103283 <https://duck.ac/submission/103283>
// 用途：上一发最好件（mmmi1k），本件是它的单变量变体；其正文与 [1] 逐字节相同（仅注释不同）。
// ======================
// ===== 思路 =====
// 单变量改动：删掉 addblk 与 comb_all 的**标量尾循环**（`for (; j < h; j++) ...`）。
// * 为什么它是死代码：这两个尾循环只处理 h % 8 个元素。本引擎三档（n=1024/2048/4096）全部走
// Strassen，addblk/comb_all 的 h 恒为 n/2，而 n 是 2 的幂且 n >= 256 ⇒ **h 恒为 8 的倍数
// ⇒ 尾循环一次都不执行**。旁证：同一文件里的 mutate_blk 本来就只有 8 路循环、没有尾循环，
// 说明本引擎已经依赖这个不变量，本发只是把它推广到另外两个函数。
// * 目的：减小这两个函数的**代码体积**（comb_all 492 → 478 条指令，addblk 153 → 145），
// 降低 L1I / DSB 压力。note: gcc 9.3 会把这两个尾循环向量化出 11 条
// `vmovdqu xmm + vinserti128`（32B 取数被拆成两条），删掉后这些也一并消失。
// * 正确性：n=128/256/512/1024/2048/4096 六个尺寸下，与出厂件在同一进程内做
// **逐元素 bitwise 对拍，bad=0**（双 namespace 同二进制，判题机同款 -O2 -static）。
// * 本次为试验性提交：目的是给"纯减代码体积"这条轴定价（预期收益很小，可能为 0）。
// ================
#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 544 /* 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")))
static __m256i LA[128] __attribute__((aligned(32)));
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 2176(%[ap]), %%ymm8, %%ymm1\n\t" "vpmaddwd 4352(%[ap]), %%ymm8, %%ymm2\n\t" "vpmaddwd 6528(%[ap]), %%ymm8, %%ymm3\n\t" "vpmaddwd 64(%[ap]), %%ymm8, %%ymm4\n\t" "vpmaddwd 2240(%[ap]), %%ymm8, %%ymm5\n\t" "vpmaddwd 4416(%[ap]), %%ymm8, %%ymm6\n\t" "vpmaddwd 6592(%[ap]), %%ymm8, %%ymm7\n\t" "vpmaddwd 0(%[ap]), %%ymm9, %%ymm10\n\t" "vpmaddwd 2176(%[ap]), %%ymm9, %%ymm11\n\t" "vpmaddwd 4352(%[ap]), %%ymm9, %%ymm12\n\t" "vpmaddwd 6528(%[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 2180(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 4356(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 6532(%[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 2244(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 4420(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 6596(%[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 2180(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 4356(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 6532(%[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 2184(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 4360(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 6536(%[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 2248(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 4424(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 6600(%[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 2184(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 4360(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 6536(%[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 2188(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 4364(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 6540(%[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 2252(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 4428(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 6604(%[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 2188(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 4364(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 6540(%[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"
"1:\n\t"
"vmovdqu 0(%[bp]), %%ymm8\n\t" "vmovdqu 128(%[bp]), %%ymm9\n\t" "vpmaddwd 0(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 2176(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 4352(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 6528(%[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 2240(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 4416(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 6592(%[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 2176(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 4352(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 6528(%[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 2180(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 4356(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 6532(%[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 2244(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 4420(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 6596(%[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 2180(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 4356(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 6532(%[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 2184(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 4360(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 6536(%[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 2248(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 4424(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 6600(%[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 2184(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 4360(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 6536(%[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 2188(%[ap]), %%ymm8, %%ymm11\n\t" "vpmaddwd 4364(%[ap]), %%ymm8, %%ymm12\n\t" "vpmaddwd 6540(%[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 2252(%[ap]), %%ymm8, %%ymm15\n\t" "vpmaddwd 4428(%[ap]), %%ymm8, %%ymm10\n\t" "vpmaddwd 6604(%[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 2188(%[ap]), %%ymm9, %%ymm13\n\t" "vpmaddwd 4364(%[ap]), %%ymm9, %%ymm14\n\t" "vpmaddwd 6540(%[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 packA(int ic, int kk, int n, int lda, const T *A) {
const __m256i m16 = _mm256_set1_epi32(0xFFFF);
const __m256i half = _mm256_set1_epi32(0x8000);
const __m256i i0 = _mm256_set1_epi32(0), i1 = _mm256_set1_epi32(1);
const __m256i i4 = _mm256_set1_epi32(4), i5 = _mm256_set1_epi32(5);
(void)n;
for (int r = 0; r < MC; r++) {
const T *src = A + (size_t)(ic + r) * lda + kk;
T *dst = Apk + (size_t)r * APSTR;
for (int m = 0; m < KC / 8; m++) {
__m256i x = LD(src + 8 * m);
__m256i al = _mm256_and_si256(x, m16);
__m256i ah = _mm256_srli_epi32(_mm256_add_epi32(x, half), 16);
__m256i pl = _mm256_packus_epi32(al, al);
__m256i ph = _mm256_packus_epi32(ah, ah);
T *d = dst + 64 * m;
ST(d + 0, _mm256_permutevar8x32_epi32(pl, i0));
ST(d + 8, _mm256_permutevar8x32_epi32(ph, i0));
ST(d + 16, _mm256_permutevar8x32_epi32(pl, i1));
ST(d + 24, _mm256_permutevar8x32_epi32(ph, i1));
ST(d + 32, _mm256_permutevar8x32_epi32(pl, i4));
ST(d + 40, _mm256_permutevar8x32_epi32(ph, i4));
ST(d + 48, _mm256_permutevar8x32_epi32(pl, i5));
ST(d + 56, _mm256_permutevar8x32_epi32(ph, i5));
}
}
}
#pragma GCC pop_options
#pragma GCC push_options
#pragma GCC optimize("O3")
AVX2 static inline void packB(int jc, int kk, int n, int ldb, const T *B) {
const __m256i m16 = _mm256_set1_epi32(0xFFFF);
const __m256i half = _mm256_set1_epi32(0x8000);
(void)n;
for (int m = 0; m < NPAIR; m++) {
const T *s0 = B + (size_t)(kk + 2 * m) * ldb + jc;
const T *s1 = B + (size_t)(kk + 2 * m + 1) * ldb + jc;
T *d0 = Bpk + (size_t)m * 16;
T *d1 = d0 + 8;
for (int j = 0; j < NC; j += 8) {
__m256i x0 = LD(s0 + j), x1 = LD(s1 + j);
__m256i l0 = _mm256_and_si256(x0, m16);
__m256i l1 = _mm256_and_si256(x1, m16);
__m256i h0 = _mm256_srli_epi32(_mm256_add_epi32(x0, half), 16);
__m256i h1 = _mm256_srli_epi32(_mm256_add_epi32(x1, half), 16);
ST(d0 + (size_t)(j / 8) * (NPAIR * 16 + 16), _mm256_or_si256(l0, _mm256_slli_epi32(l1, 16)));
ST(d1 + (size_t)(j / 8) * (NPAIR * 16 + 16), _mm256_or_si256(h0, _mm256_slli_epi32(h1, 16)));
}
}
}
#pragma GCC pop_options
#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);
/* RIDX turns packusdw(a,a)=[p0,p1,p0,p1,p2,p3,p2,p3] into [p0,p1,p2,p3,p0,p1,p2,p3] */
const __m256i ridx = _mm256_setr_epi32(0, 1, 4, 5, 0, 1, 4, 5);
(void)n;
for (int r = 0; r < MC; r++) {
const T *src = A + (size_t)(ic + r) * lda + kk;
const T *src2 = A2 ? (A2 + (size_t)(ic + r) * lda + kk) : 0;
T *dst = Apk + (size_t)r * APSTR;
for (int u = 0; u < KC / 8; u++) {
__m256i x = LD(src + 8 * u);
if (src2) { __m256i y = LD(src2 + 8 * u);
x = (sg > 0) ? _mm256_add_epi32(x, y) : _mm256_sub_epi32(x, y); }
__m256i al = _mm256_and_si256(x, m16);
__m256i ah = _mm256_srli_epi32(_mm256_add_epi32(x, 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 = dst + 32 * u; /* 128 B per 8 k per row */
ST(d + 0, rl); /* al slots 0..7 */
_mm_store_si128((__m128i *)(d + 8), _mm256_castsi256_si128(rl)); /* 8..11 */
ST(d + 16, rh); /* ah slots 0..7 (byte 64) */
_mm_store_si128((__m128i *)(d + 24), _mm256_castsi256_si128(rh)); /* 8..11 */
}
}
}
#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 m16 = _mm256_set1_epi32(0xFFFF);
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 + 16) + (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_or_si256(_mm256_and_si256(x0, m16), \
_mm256_slli_epi32(_mm256_and_si256(x1, m16), 16)); \
HI = _mm256_or_si256(_mm256_srli_epi32(_mm256_add_epi32(x0, half), 16), \
_mm256_slli_epi32( \
_mm256_srli_epi32(_mm256_add_epi32(x1, half), 16), 16)); }
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 + 16);
}
}
}
/* 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
/* dst = s1*m1 + s2*m2 , ONE pass (M blocks are contiguous, row stride h) */
#pragma GCC push_options
#pragma GCC optimize("O3")
AVX2 static void comb2(int h, T * __restrict dst, int ldd, const T * __restrict m1, int s1, const T * __restrict m2, int s2) {
(void)s1;
for (int i = 0; i < h; i++) {
const T *p1 = m1 + (size_t)i * h, *p2 = m2 + (size_t)i * h;
T *dp = dst + (size_t)i * ldd;
int j = 0;
for (; j + 8 <= h; j += 8) {
__m256i v = LDA(p1 + j);
v = (s2 > 0) ? _mm256_add_epi32(v, LDA(p2 + j)) : _mm256_sub_epi32(v, LDA(p2 + j));
_mm256_stream_si256((__m256i*)(dp + j), v);
}
for (; j < h; j++) dp[j] = (T)(p1[j] + (s2 > 0 ? p2[j] : -p2[j]));
}
}
#pragma GCC pop_options
/* dst = s1*m1 + s2*m2 + s3*m3 + s4*m4 , ONE pass */
#pragma GCC push_options
#pragma GCC optimize("O3")
AVX2 static void comb4(int h, T * __restrict dst, int ldd, const T * __restrict m1, int s1, const T * __restrict m2, int s2,
const T * __restrict m3, int s3, const T * __restrict m4, int s4) {
(void)s1;
for (int i = 0; i < h; i++) {
const T *p1 = m1 + (size_t)i * h, *p2 = m2 + (size_t)i * h,
*p3 = m3 + (size_t)i * h, *p4 = m4 + (size_t)i * h;
T *dp = dst + (size_t)i * ldd;
int j = 0;
for (; j + 8 <= h; j += 8) {
__m256i v = LDA(p1 + j);
v = (s2 > 0) ? _mm256_add_epi32(v, LDA(p2 + j)) : _mm256_sub_epi32(v, LDA(p2 + j));
v = (s3 > 0) ? _mm256_add_epi32(v, LDA(p3 + j)) : _mm256_sub_epi32(v, LDA(p3 + j));
v = (s4 > 0) ? _mm256_add_epi32(v, LDA(p4 + j)) : _mm256_sub_epi32(v, LDA(p4 + j));
_mm256_stream_si256((__m256i*)(dp + j), v);
}
for (; j < h; j++)
dp[j] = (T)(p1[j] + (s2>0?p2[j]:-p2[j]) + (s3>0?p3[j]:-p3[j]) + (s4>0?p4[j]:-p4[j]));
}
}
#pragma GCC pop_options
#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) { 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(B11, B22, +1, D11, D22, +1, M1); /* M1 = (A11+A22)(B11+B22) */
LOP(B21, B22, +1, D11, 0, +1, M2); /* M2 = (A21+A22) B11 */
LOP(B11, 0, +1, D12, D22, -1, M3); /* M3 = A11 (B12-B22) */
LOP(B22, 0, +1, D21, D11, -1, M4); /* M4 = A22 (B21-B11) */
LOP(B11, B12, +1, D22, 0, +1, M5); /* M5 = (A11+A12) B22 */
LOP(B21, B11, -1, D11, D12, +1, M6); /* M6 = (A21-A11)(B11+B12) */
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));
}
}
}
#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 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,A11,lda,B11,ldb,M1,h,sub,depth+1); // P1
rec_w(h,A12,lda,B21,ldb,M2,h,sub,depth+1); // P2
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,h,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,h,sub,depth+1); // P6
mutate_blk(h,P,A12,lda,1); // S4 = A12-S2
rec_w(h,P,h,B22,ldb,M3,h,sub,depth+1); // P3
mutate_blk(h,Q,B21,ldb,0); // T4 = T2-B21
rec_w(h,A22,lda,Q,h,M4,h,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,h,sub,depth+1); // P7
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 | 249.996 ms | 63 MB + 868 KB | Accepted | Score: 100 | 显示更多 |