#define MK5O 1
#include <immintrin.h>
#include <stddef.h>
#include <stdint.h>
static unsigned long long PHK[5], PHKC[5], _mk5t;
static unsigned long long PHN[5], PHNC[5];
static inline unsigned long long RDN(void){ unsigned l,hh; __asm__ volatile("rdtscp":"=a"(l),"=d"(hh)::"rcx"); return ((unsigned long long)hh<<32)|l; }
static inline void MK5B(void){ unsigned l,hh; __asm__ volatile("rdtscp":"=a"(l),"=d"(hh)::"rcx"); _mk5t=((unsigned long long)hh<<32)|l; }
static int MK5IDX(int h){ return h==2048?0:(h==1024?1:(h==512?2:(h==256?3:4))); }
static inline void MK5E(int hh){ unsigned l,h2; int i=MK5IDX(hh); __asm__ volatile("rdtscp":"=a"(l),"=d"(h2)::"rcx"); PHK[i]+=(((unsigned long long)h2<<32)|l)-_mk5t; PHKC[i]++; }
static int MMNIDX(int n){ return n>=4096?0:(n>=2048?1:(n>=1024?2:3)); }
static inline void MMNA(int n, unsigned long long t0){ int i=MMNIDX(n); PHN[i]+=RDN()-t0; PHNC[i]++; }
#include <stdlib.h>
extern "C" int printf(const char *, ...); /* OF-KNOB: state probe */
#pragma GCC optimize("-fno-crossjumping")
#pragma GCC push_options
#pragma GCC target("avx2")
#ifndef ORD
#define ORD 1
#endif
typedef short TE;
#define PAD 32
static TE *g_cc[4];
static int g_ch = 0, g_cok = 0;
static TE *g_rc[4];
static int g_rh = 0;
struct Knobs { volatile int copy; volatile int use; volatile int alloc; volatile int wmode;
volatile int copy2; volatile int use2;
volatile int wA11; volatile int wA22; volatile int wB11; volatile int wB22;
volatile int lv; volatile int bp; volatile int abl;
volatile int off; volatile int offc; volatile int dbg; volatile int sp; };
static struct Knobs g_k = { 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, -1 }; /* OFA0: alloc=0 */
#define g_kbp (g_k.bp)
#define g_abl (g_k.abl)
#define g_copy (g_k.copy)
#define g_use (g_k.use)
#define g_alloc (g_k.alloc)
#define g_wmode (g_k.wmode)
#define g_copy2 (g_k.copy2)
#define g_use2 (g_k.use2)
#define g_wA11 (g_k.wA11)
#define g_wA22 (g_k.wA22)
#define g_wB11 (g_k.wB11)
#define g_wB22 (g_k.wB22)
#define g_lv (g_k.lv)
#ifndef CUT
#define CUT 256
#endif
#define MR 5
#define NR 16
#define KC 256
#define NC 384
#define UNROLL 4
#define TRS (((KC) + 15) & ~15)
#define BPAD 32
#define STV(p, v) _mm256_store_si256((__m256i *)(p), (v))
#define LDV(p) _mm256_load_si256((const __m256i *)(p))
#define LDU(p) _mm256_loadu_si256((const __m256i *)(p))
struct __attribute__((aligned(64))) { char _pad[4096]; TE _b[(size_t)NC * KC + (size_t)(NC / NR) * BPAD + 256]; } g_bp;
#define Bpanel g_bp._b
__attribute__((used, noinline)) static void m98_padfn(void)
{ __asm__ volatile(".fill 64,1,0x90"); }
static inline __m256i b32(const TE *p) {
__m256i r; __asm__("vpbroadcastd %1, %0" : "=x"(r) : "m"(*(const int *)p)); return r;
}
static inline __m256i bilv(const TE *lo, const TE *hi) {
return _mm256_unpacklo_epi16(_mm256_castsi128_si256(_mm_loadu_si128((const __m128i *)lo)),
_mm256_castsi128_si256(_mm_loadu_si128((const __m128i *)hi)));
}
#define KST(off, C0, C1) \
"vpbroadcastd " #off "(%[a]), %%ymm12\n\t" \
"vpmaddwd %%ymm12, %%ymm14, %%ymm13\n\t" \
"vpaddd %%ymm13, %" C0 ", %" C0 "\n\t" \
"vpmaddwd %%ymm12, %%ymm15, %%ymm13\n\t" \
"vpaddd %%ymm13, %" C1 ", %" C1 "\n\t"
#define KST5D(a5) \
"vpbroadcastd " #a5 "(%[a]), %%ymm12\n\t" \
"vpmaddwd %%ymm12, %%ymm14, %%ymm13\n\t" \
"vpmaddwd %%ymm12, %%ymm15, %%ymm12\n\t"
#define KSTF(boff, off, C0, C1) \
"vpbroadcastd " #off "(%[a]), %%ymm12\n\t" \
"vpmaddwd " #boff "(%[b]), %%ymm12, %%ymm13\n\t" \
"vpaddd %%ymm13, %" C0 ", %" C0 "\n\t" \
"vpmaddwd %%ymm12, %%ymm15, %%ymm13\n\t" \
"vpaddd %%ymm13, %" C1 ", %" C1 "\n\t"
#define KST5DF(boff, a5) \
"vpbroadcastd " #a5 "(%[a]), %%ymm12\n\t" \
"vpmaddwd " #boff "(%[b]), %%ymm12, %%ymm13\n\t" \
"vpmaddwd %%ymm12, %%ymm15, %%ymm12\n\t"
#define KGRP_P2(a0,a1,a2,a3,a4,a5,b0,b1) \
"vmovdqa " #b1 "(%[b]), %%ymm15\n\t" \
"vpaddd %%ymm13, %10, %10\n\t" \
"vpaddd %%ymm12, %11, %11\n\t" \
KSTF(b0,a0,"0","1") KSTF(b0,a1,"2","3") KSTF(b0,a2,"4","5") \
KSTF(b0,a3,"6","7") KSTF(b0,a4,"8","9") KST5DF(b0,a5)
#define KGRP(a0,a1,a2,a3,a4,a5,b0,b1) \
"vmovdqa " #b0 "(%[b]), %%ymm14\n\t" \
"vmovdqa " #b1 "(%[b]), %%ymm15\n\t" \
KST(a0,"0","1") KST(a1,"2","3") KST(a2,"4","5") KST(a3,"6","7") KST(a4,"8","9") KST(a5,"10","11")
#define TRB (TRS * 2)
#if TRS != 256
#error "compact A operand hard-codes the tr row stride in the asm immediates"
#endif
#define KGRP1 KGRP(0,512,1024,1536,2048,2560, 0,32)
#define PF(o) "prefetcht0 " #o "(%[b])\n\t"
#define KFULL PF(512) KGRP_P2(0,512,1024,1536,2048,2560, 0,32) \
PF(576) KGRP_P2(4,516,1028,1540,2052,2564, 64,96) \
PF(640) KGRP_P2(8,520,1032,1544,2056,2568, 128,160) \
PF(704) KGRP_P2(12,524,1036,1548,2060,2572, 192,224)
#define MICRO_CLOB "ymm12", "ymm13", "ymm14", "ymm15", "cc", "memory"
#define FOLD16(v) _mm256_permute4x64_epi64(_mm256_shuffle_epi8((v), FOLDM), 0x08)
static const __m256i FOLDM = { (long long)0x0D0C090805040100LL, (long long)0x8080808080808080LL,
(long long)0x0D0C090805040100LL, (long long)0x8080808080808080LL };
#define FH(v) _mm256_shuffle_epi8((v), FOLDM)
/* lane mm1k pk2: LEGAL narrowing of the accumulator pair -> one 32 B store row.
The shipped form is 2x vpshufb + 1x vpunpcklqdq = 3 p5 uops per row (18 per mikro call).
vpackusdw saturates UNSIGNED, so masking the int32 lanes to [0,65535] with vpand (p015)
makes the pack exact: it then reproduces the low 16 bits of each int32 lane, which is the
row's contract (C is `short`, the sum is mod 2^16). Lane order is provably identical. */
#define PK2M _mm256_set1_epi32(0x0000FFFF)
#define PACKROW(v0, v1) _mm256_packus_epi32(_mm256_and_si256((v0), PK2M), _mm256_and_si256((v1), PK2M))
#define ST_ROW(r, v0, v1) \
if (rows > (r)) { \
TE *p0 = C + (size_t)(r) * ldc; \
__m256i zz = PACKROW(v0, v1); \
if (OP == 0) { STV(p0, zz); } \
else { STV(p0, _mm256_add_epi16(LDV(p0), zz)); } \
}
#define MICRO_DEF(NAME, BODY, ADV_A, ADV_B) \
template <int OP> \
static inline void NAME(long k, const TE *ap, const TE *bp, TE *C, int ldc, int rows) { \
__m256i c00 = _mm256_setzero_si256(), c01 = c00, c10 = c00, c11 = c00, c20 = c00, c21 = c00; \
__m256i c30 = c00, c31 = c00, c40 = c00, c41 = c00, c50 = c00, c51 = c00; \
__asm__ volatile( \
"vpxor %%ymm12, %%ymm12, %%ymm12\n\t" \
"vpxor %%ymm13, %%ymm13, %%ymm13\n\t" \
"1:\n\t" BODY \
"addq $" #ADV_A ", %[a]\n\t" \
"addq $" #ADV_B ", %[b]\n\t" \
"decq %[k]\n\t" \
"jnz 1b\n\t" \
"vpaddd %%ymm13, %10, %10\n\t" \
"vpaddd %%ymm12, %11, %11\n\t" \
: "+x"(c00), "+x"(c01), "+x"(c10), "+x"(c11), "+x"(c20), "+x"(c21), \
"+x"(c30), "+x"(c31), "+x"(c40), "+x"(c41), "+x"(c50), "+x"(c51), \
[a] "+r"(ap), [b] "+r"(bp), [k] "+r"(k) \
: \
: MICRO_CLOB); \
ST_ROW(0, c00, c01) ST_ROW(1, c10, c11) ST_ROW(2, c20, c21) \
ST_ROW(3, c30, c31) ST_ROW(4, c40, c41) ST_ROW(5, c50, c51) \
}
/* =====================================================================================
lane_mmm8_l1 arm MVQ -- BOTH vpmaddwd OF A ROW ISSUED ADJACENTLY, THEN BOTH vpaddd.
The shipped body is madd/add/madd/add per row; this is madd/madd/add/add. SAME
operands, SAME instructions, SAME accumulator updates -- ONLY the issue order inside
the 5-instruction group changes, so it is byte-exact by construction and was verified
so (identical full-output hash #ea4285 against the unmodified body in one binary).
ymm12 is free once the broadcast has been consumed by the first madd, so the second
madd's result goes there and no accumulator is displaced.
WHY IT PAYS (judge, real judge, arms in ONE .text-identical binary, 8 ROTATED ROUNDS):
ctl 42 739 338 microQ 42 353 592 = -0.903 %
and on a separate 6-arm call with THREE independent control slots the controls agreed
to within 0.02 % (42 514 804 / 42 516 410 / 42 532 148), so this is ~45x the channel.
The port arithmetic says why there is anything to win at all: per k-pair the body is
12 vpmaddwd (p01 ONLY, so >= 6 cycles) + 12 vpaddd (p015); the exact optimum splits the
adds 8 to p5 and 4 to p01 (p01 = 16 uops = 8 cycles, p5 = 8 cyc). Putting the two
madds adjacent gives the scheduler a longer, cleaner run to find that split from,
instead of the shipped alternation where every add's oldest-ready madd sits 5 cycles
behind it. Measured 34.3 ticks/KFULL against a hard 32-cycle 3-port floor; this arm
moves that by 0.9 % of the whole row, i.e. ~1 % of the mikro. */
#define KSTQ2(boff, off, C0, C1) \
"vpbroadcastd " #off "(%[a]), %%ymm12\n\t" \
"vpmaddwd " #boff "(%[b]), %%ymm12, %%ymm13\n\t" \
"vpmaddwd %%ymm12, %%ymm15, %%ymm14\n\t" \
"vpaddd %%ymm13, %" C0 ", %" C0 "\n\t" \
"vpaddd %%ymm14, %" C1 ", %" C1 "\n\t"
#define KGRP_Q(a0,a1,a2,a3,a4,a5,b0,b1) \
"vmovdqa " #b1 "(%[b]), %%ymm15\n\t" \
"vpaddd %%ymm13, %10, %10\n\t" \
"vpaddd %%ymm12, %11, %11\n\t" \
KSTQ2(b0,a0,"0","1") KSTQ2(b0,a1,"2","3") KSTQ2(b0,a2,"4","5") \
KSTQ2(b0,a3,"6","7") KSTQ2(b0,a4,"8","9") KST5DF(b0,a5)
#define KFULL_Q PF(512) KGRP_Q(0,512,1024,1536,2048,2560, 0,32) \
PF(576) KGRP_Q(4,516,1028,1540,2052,2564, 64,96) \
PF(640) KGRP_Q(8,520,1032,1544,2056,2568, 128,160) \
PF(704) KGRP_Q(12,524,1036,1548,2060,2572, 192,224)
#define ST_ROW_S(r, v0, v1) \
if (ROWS_ > (r)) { \
TE *p0 = C + (size_t)(r) * ldc; \
__m256i zz = PACKROW(v0, v1); \
if (OP == 0) { STV(p0, zz); } \
else { STV(p0, _mm256_add_epi16(LDV(p0), zz)); } \
}
#define KH_BC(o) "vpbroadcastd " #o "(%[a]), %%ymm12\n\t"
#define KH_MAD(bo) "vpmaddwd " #bo "(%[b]), %%ymm12, %%ymm13\n\t" "vpmaddwd %%ymm12, %%ymm15, %%ymm14\n\t"
#define KH_ADD(C0,C1) "vpaddd %%ymm13, %" C0 ", %" C0 "\n\t" "vpaddd %%ymm14, %" C1 ", %" C1 "\n\t"
#define KH_MAD5D(bo) "vpmaddwd " #bo "(%[b]), %%ymm12, %%ymm13\n\t" "vpmaddwd %%ymm12, %%ymm15, %%ymm12\n\t"
#define KGRP_H(a0,a1,a2,a3,a4,a5,b0,b1) \
"vmovdqa " #b1 "(%[b]), %%ymm15\n\t" \
"vpaddd %%ymm13, %10, %10\n\t" \
"vpaddd %%ymm12, %11, %11\n\t" \
KH_BC(a0) KH_MAD(b0) KH_BC(a1) KH_ADD("0","1") KH_MAD(b0) KH_BC(a2) KH_ADD("2","3") \
KH_MAD(b0) KH_BC(a3) KH_ADD("4","5") KH_MAD(b0) KH_BC(a4) KH_ADD("6","7") \
KH_MAD(b0) KH_BC(a5) KH_ADD("8","9") KH_MAD5D(b0)
#define KFULL_H PF(768) KGRP_H(0,512,1024,1536,2048,2560, 0,32) \
PF(832) KGRP_H(4,516,1028,1540,2052,2564, 64,96) \
PF(896) KGRP_H(8,520,1032,1544,2056,2568, 128,160) \
PF(960) KGRP_H(12,524,1036,1548,2060,2572, 192,224)
/* =====================================================================================
lane mm1k_r4: THE (5,16) MIKRO. Ten accumulators (%%0..%%9) = 5 rows x 2 column halves.
ymm10 = the A broadcast, ymm11 = the b1 register, ymm12..ymm15 = FOUR madd-result
registers. Two rows are in flight at once, so each add sits 4..9 instructions after the
madd that produces it, against the shipped (6,16) body's uniform distance of 3 (0.75
cycles) for a 5-cycle madd latency. ZERO drains. 26 uops / 160 MACs = 6.15 MACs/uop
against (6,16)'s 31 / 192 = 6.19; the p01 floor is IDENTICAL.
===================================================================================== */
#define MICRO_CLOB5 "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15", "cc", "memory"
#define KGRP5_0 \
"vmovdqa 32(%[b]), %%ymm11\n\t" \
"vpbroadcastd 0(%[a]), %%ymm10\n\t" \
"vpmaddwd 0(%[b]), %%ymm10, %%ymm12\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm13\n\t" \
"vpaddd %%ymm12, %0, %0\n\t" \
"vpbroadcastd 512(%[a]), %%ymm10\n\t" \
"vpmaddwd 0(%[b]), %%ymm10, %%ymm14\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm15\n\t" \
"vpaddd %%ymm13, %1, %1\n\t" \
"vpbroadcastd 1024(%[a]), %%ymm10\n\t" \
"vpmaddwd 0(%[b]), %%ymm10, %%ymm12\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm13\n\t" \
"vpaddd %%ymm14, %2, %2\n\t" \
"vpaddd %%ymm15, %3, %3\n\t" \
"vpbroadcastd 1536(%[a]), %%ymm10\n\t" \
"vpmaddwd 0(%[b]), %%ymm10, %%ymm14\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm15\n\t" \
"vpaddd %%ymm12, %4, %4\n\t" \
"vpaddd %%ymm13, %5, %5\n\t" \
"vpbroadcastd 2048(%[a]), %%ymm10\n\t" \
"vpmaddwd 0(%[b]), %%ymm10, %%ymm12\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm13\n\t" \
"vpaddd %%ymm14, %6, %6\n\t" \
"vpaddd %%ymm15, %7, %7\n\t" \
"vpaddd %%ymm12, %8, %8\n\t" \
"vpaddd %%ymm13, %9, %9\n\t"
#define KGRP5_1 \
"vmovdqa 96(%[b]), %%ymm11\n\t" \
"vpbroadcastd 4(%[a]), %%ymm10\n\t" \
"vpmaddwd 64(%[b]), %%ymm10, %%ymm12\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm13\n\t" \
"vpaddd %%ymm12, %0, %0\n\t" \
"vpbroadcastd 516(%[a]), %%ymm10\n\t" \
"vpmaddwd 64(%[b]), %%ymm10, %%ymm14\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm15\n\t" \
"vpaddd %%ymm13, %1, %1\n\t" \
"vpbroadcastd 1028(%[a]), %%ymm10\n\t" \
"vpmaddwd 64(%[b]), %%ymm10, %%ymm12\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm13\n\t" \
"vpaddd %%ymm14, %2, %2\n\t" \
"vpaddd %%ymm15, %3, %3\n\t" \
"vpbroadcastd 1540(%[a]), %%ymm10\n\t" \
"vpmaddwd 64(%[b]), %%ymm10, %%ymm14\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm15\n\t" \
"vpaddd %%ymm12, %4, %4\n\t" \
"vpaddd %%ymm13, %5, %5\n\t" \
"vpbroadcastd 2052(%[a]), %%ymm10\n\t" \
"vpmaddwd 64(%[b]), %%ymm10, %%ymm12\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm13\n\t" \
"vpaddd %%ymm14, %6, %6\n\t" \
"vpaddd %%ymm15, %7, %7\n\t" \
"vpaddd %%ymm12, %8, %8\n\t" \
"vpaddd %%ymm13, %9, %9\n\t"
#define KGRP5_2 \
"vmovdqa 160(%[b]), %%ymm11\n\t" \
"vpbroadcastd 8(%[a]), %%ymm10\n\t" \
"vpmaddwd 128(%[b]), %%ymm10, %%ymm12\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm13\n\t" \
"vpaddd %%ymm12, %0, %0\n\t" \
"vpbroadcastd 520(%[a]), %%ymm10\n\t" \
"vpmaddwd 128(%[b]), %%ymm10, %%ymm14\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm15\n\t" \
"vpaddd %%ymm13, %1, %1\n\t" \
"vpbroadcastd 1032(%[a]), %%ymm10\n\t" \
"vpmaddwd 128(%[b]), %%ymm10, %%ymm12\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm13\n\t" \
"vpaddd %%ymm14, %2, %2\n\t" \
"vpaddd %%ymm15, %3, %3\n\t" \
"vpbroadcastd 1544(%[a]), %%ymm10\n\t" \
"vpmaddwd 128(%[b]), %%ymm10, %%ymm14\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm15\n\t" \
"vpaddd %%ymm12, %4, %4\n\t" \
"vpaddd %%ymm13, %5, %5\n\t" \
"vpbroadcastd 2056(%[a]), %%ymm10\n\t" \
"vpmaddwd 128(%[b]), %%ymm10, %%ymm12\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm13\n\t" \
"vpaddd %%ymm14, %6, %6\n\t" \
"vpaddd %%ymm15, %7, %7\n\t" \
"vpaddd %%ymm12, %8, %8\n\t" \
"vpaddd %%ymm13, %9, %9\n\t"
#define KGRP5_3 \
"vmovdqa 224(%[b]), %%ymm11\n\t" \
"vpbroadcastd 12(%[a]), %%ymm10\n\t" \
"vpmaddwd 192(%[b]), %%ymm10, %%ymm12\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm13\n\t" \
"vpaddd %%ymm12, %0, %0\n\t" \
"vpbroadcastd 524(%[a]), %%ymm10\n\t" \
"vpmaddwd 192(%[b]), %%ymm10, %%ymm14\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm15\n\t" \
"vpaddd %%ymm13, %1, %1\n\t" \
"vpbroadcastd 1036(%[a]), %%ymm10\n\t" \
"vpmaddwd 192(%[b]), %%ymm10, %%ymm12\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm13\n\t" \
"vpaddd %%ymm14, %2, %2\n\t" \
"vpaddd %%ymm15, %3, %3\n\t" \
"vpbroadcastd 1548(%[a]), %%ymm10\n\t" \
"vpmaddwd 192(%[b]), %%ymm10, %%ymm14\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm15\n\t" \
"vpaddd %%ymm12, %4, %4\n\t" \
"vpaddd %%ymm13, %5, %5\n\t" \
"vpbroadcastd 2060(%[a]), %%ymm10\n\t" \
"vpmaddwd 192(%[b]), %%ymm10, %%ymm12\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm13\n\t" \
"vpaddd %%ymm14, %6, %6\n\t" \
"vpaddd %%ymm15, %7, %7\n\t" \
"vpaddd %%ymm12, %8, %8\n\t" \
"vpaddd %%ymm13, %9, %9\n\t"
#define KFULL5 PF(768) KGRP5_0 PF(832) KGRP5_1 PF(896) KGRP5_2 PF(960) KGRP5_3
template <int OP, int ROWS_>
static inline void mM5(long k, const TE *ap, const TE *bp, TE *C, int ldc, int rows) {
(void)rows;
__m256i c0 = _mm256_setzero_si256(), c1 = c0, c2 = c0, c3 = c0, c4 = c0;
__m256i c5 = c0, c6 = c0, c7 = c0, c8 = c0, c9 = c0;
__asm__ volatile(
"vpxor %%ymm12, %%ymm12, %%ymm12\n\t"
"vpxor %%ymm13, %%ymm13, %%ymm13\n\t"
"1:\n\t" KFULL5
"addq $16, %[a]\n\t"
"addq $256, %[b]\n\t"
"decq %[k]\n\t"
"jnz 1b\n\t"
: "+x"(c0), "+x"(c1), "+x"(c2), "+x"(c3), "+x"(c4),
"+x"(c5), "+x"(c6), "+x"(c7), "+x"(c8), "+x"(c9),
[a] "+r"(ap), [b] "+r"(bp), [k] "+r"(k)
: : MICRO_CLOB5);
ST_ROW_S(0, c0, c1) ST_ROW_S(1, c2, c3) ST_ROW_S(2, c4, c5)
ST_ROW_S(3, c6, c7) ST_ROW_S(4, c8, c9)
}
/* 1-row tail mikro for MR = 5 (256 mod 5 = 1). micro4 would compute 4 rows for it, wasting
3 of 256 rows = +1.17 percent of the mikro, which is more than the 0.6 percent the geometry
costs. 21 uops / 128 MACs is deliberately not tuned: it is 1/256 of the rows. */
#define R1K(AO,BO,B1) \
"vmovdqa " #B1 "(%[b]), %%ymm11\n\t" \
"vpbroadcastd " #AO "(%[a]), %%ymm10\n\t" \
"vpmaddwd " #BO "(%[b]), %%ymm10, %%ymm13\n\t" \
"vpaddd %%ymm13, %0, %0\n\t" \
"vpmaddwd %%ymm10, %%ymm11, %%ymm12\n\t" \
"vpaddd %%ymm12, %1, %1\n\t"
template <int OP, int ROWS_>
static inline void microR1(long k, const TE *ap, const TE *bp, TE *C, int ldc, int rows) {
(void)rows;
__m256i c0 = _mm256_setzero_si256(), c1 = c0;
__asm__ volatile(
"1:\n\t"
R1K(0,0,32) R1K(4,64,96) R1K(8,128,160) R1K(12,192,224)
"addq $16, %[a]\n\t"
"addq $256, %[b]\n\t"
"decq %[k]\n\t"
"jnz 1b\n\t"
: "+x"(c0), "+x"(c1), [a] "+r"(ap), [b] "+r"(bp), [k] "+r"(k)
: : "ymm10", "ymm11", "ymm12", "ymm13", "cc", "memory");
ST_ROW_S(0, c0, c1)
}
template <int OP, int ROWS_>
static inline void mM(long k, const TE *ap, const TE *bp, TE *C, int ldc, int rows) {
(void)rows;
__m256i c00 = _mm256_setzero_si256(), c01 = c00, c10 = c00, c11 = c00, c20 = c00, c21 = c00;
__m256i c30 = c00, c31 = c00, c40 = c00, c41 = c00, c50 = c00, c51 = c00;
__asm__ volatile(
"vpxor %%ymm12, %%ymm12, %%ymm12\n\t"
"vpxor %%ymm13, %%ymm13, %%ymm13\n\t"
"1:\n\t" KFULL_H
"addq $16, %[a]\n\t"
"addq $256, %[b]\n\t"
"decq %[k]\n\t"
"jnz 1b\n\t"
"vpaddd %%ymm13, %10, %10\n\t"
"vpaddd %%ymm12, %11, %11\n\t"
: "+x"(c00), "+x"(c01), "+x"(c10), "+x"(c11), "+x"(c20), "+x"(c21),
"+x"(c30), "+x"(c31), "+x"(c40), "+x"(c41), "+x"(c50), "+x"(c51),
[a] "+r"(ap), [b] "+r"(bp), [k] "+r"(k)
: : MICRO_CLOB);
ST_ROW_S(0, c00, c01) ST_ROW_S(1, c10, c11) ST_ROW_S(2, c20, c21)
ST_ROW_S(3, c30, c31) ST_ROW_S(4, c40, c41) ST_ROW_S(5, c50, c51)
}
MICRO_DEF(micro, KFULL_Q, 16, 256)
MICRO_DEF(micro1, KGRP1, 4, 64)
static inline __m256i bp2x(const TE *p0, const TE *p1, const TE *q0, const TE *q1, int op, int hi) {
__m128i a0 = _mm_loadu_si128((const __m128i *)p0);
__m128i b0 = _mm_loadu_si128((const __m128i *)p1);
__m128i a1 = _mm_loadu_si128((const __m128i *)(p0 + 8));
__m128i b1 = _mm_loadu_si128((const __m128i *)(p1 + 8));
if (op == 1) { a0 = _mm_add_epi16(a0, _mm_loadu_si128((const __m128i *)q0));
b0 = _mm_add_epi16(b0, _mm_loadu_si128((const __m128i *)q1));
a1 = _mm_add_epi16(a1, _mm_loadu_si128((const __m128i *)(q0 + 8)));
b1 = _mm_add_epi16(b1, _mm_loadu_si128((const __m128i *)(q1 + 8))); }
else if (op == 2) { a0 = _mm_sub_epi16(a0, _mm_loadu_si128((const __m128i *)q0));
b0 = _mm_sub_epi16(b0, _mm_loadu_si128((const __m128i *)q1));
a1 = _mm_sub_epi16(a1, _mm_loadu_si128((const __m128i *)(q0 + 8)));
b1 = _mm_sub_epi16(b1, _mm_loadu_si128((const __m128i *)(q1 + 8))); }
if (!hi) return _mm256_inserti128_si256(_mm256_castsi128_si256(_mm_unpacklo_epi16(a0, b0)),
_mm_unpacklo_epi16(a1, b1), 1);
return _mm256_inserti128_si256(_mm256_castsi128_si256(_mm_unpackhi_epi16(a0, b0)),
_mm_unpackhi_epi16(a1, b1), 1);
}
#define OP_RAW 0
#define OP_ADD 1
#define OP_SUB 2
static void trb_tail(TE tr[MR][TRS], int rows, int kc, int k0,
const TE *w1, const TE *w2, const TE *w3, const TE *w4, int lda, int lda2,
int op1, int op2) {
for (int k = k0; k < kc; k++) {
for (int r = 0; r < rows; r++) {
int v = (int)w1[(size_t)r * lda + k];
if (op1 != OP_RAW) {
int y = (int)w2[(size_t)r * lda2 + k];
v += (op1 == OP_ADD) ? y : -y;
}
if (op2 != OP_RAW) {
int y = (int)w3[(size_t)r * lda2 + k];
if (op1 != OP_RAW) {
int z = (int)w4[(size_t)r * lda2 + k];
y += (op1 == OP_ADD) ? z : -z;
}
v += (op2 == OP_ADD) ? y : -y;
}
tr[r][k] = (TE)v;
}
}
}
template <int OP1, int OP2, int ROWS = -1, int KCS = -1>
static inline void trb(TE tr[MR][TRS], int rows_, int kc_,
const TE *w1, const TE *w2, const TE *w3, const TE *w4,
int lda, int lda2) {
const int rows = (ROWS >= 0) ? ROWS : rows_;
const int kc = (KCS > 0) ? KCS : kc_;
int k = 0;
#pragma GCC unroll 8
for (; k + 32 <= kc; k += 32, w1 += 32, w2 += 32, w3 += 32, w4 += 32)
for (int r = 0; r < rows; r++) {
const TE *s1 = w1 + (size_t)r * lda;
__m256i x0 = LDV(s1), x1 = LDV(s1 + 16);
if (OP1 != OP_RAW) {
const TE *s2 = w2 + (size_t)r * lda2;
if (OP1 == OP_ADD) { x0 = _mm256_add_epi16(x0, LDV(s2)); x1 = _mm256_add_epi16(x1, LDV(s2 + 16)); }
else { x0 = _mm256_sub_epi16(x0, LDV(s2)); x1 = _mm256_sub_epi16(x1, LDV(s2 + 16)); }
}
if (OP2 == OP_RAW) {
STV(tr[r] + k, x0); STV(tr[r] + k + 16, x1);
} else {
const TE *s3 = w3 + (size_t)r * lda2;
__m256i y0 = LDV(s3), y1 = LDV(s3 + 16);
if (OP1 != OP_RAW) {
const TE *s4 = w4 + (size_t)r * lda2;
if (OP1 == OP_ADD) { y0 = _mm256_add_epi16(y0, LDV(s4)); y1 = _mm256_add_epi16(y1, LDV(s4 + 16)); }
else { y0 = _mm256_sub_epi16(y0, LDV(s4)); y1 = _mm256_sub_epi16(y1, LDV(s4 + 16)); }
}
if (OP2 == OP_ADD) { STV(tr[r] + k, _mm256_add_epi16(x0, y0)); STV(tr[r] + k + 16, _mm256_add_epi16(x1, y1)); }
else { STV(tr[r] + k, _mm256_sub_epi16(x0, y0)); STV(tr[r] + k + 16, _mm256_sub_epi16(x1, y1)); }
}
}
if (k < kc) trb_tail(tr, rows, kc, k, w1, w2, w3, w4, lda, lda2, OP1, OP2);
}
static void packb_tail(const TE *B, int ldb, const TE *B2, const TE *B3, const TE *B4, int ldb2,
int kc, int nc, int pc, int jc, int jr, int bop, int bop2) {
const int kn = kc >> 1;
for (; jr < nc; jr += NR) {
TE *q = Bpanel + (size_t)(jr / NR) * ((size_t)kc * NR + BPAD);
for (int kp = 0; kp < kn; kp++) {
const int k = kp << 1;
for (int vec = 0; vec < 4; vec++) {
TE *dst = q + (vec & 1) * 16 + (size_t)(vec >> 1) * ((size_t)kc * NR + BPAD);
const int colbase = jc + jr + (vec >> 1) * 16 + ((vec & 1) ? 8 : 0);
for (int t = 0; t < 16; t++) {
const int col = colbase + (t >> 1);
const int rw = pc + k + (t & 1);
int a = (int)B[(size_t)rw * ldb + col];
if (bop != OP_RAW) {
int x = (int)B2[(size_t)rw * ldb2 + col];
a += (bop == OP_ADD) ? x : -x;
}
if (bop2 != OP_RAW) {
int x = (int)B3[(size_t)rw * ldb2 + col];
if (bop != OP_RAW) {
int x2 = (int)B4[(size_t)rw * ldb2 + col];
x += (bop == OP_ADD) ? x2 : -x2;
}
a += (bop2 == OP_ADD) ? x : -x;
}
dst[t] = (TE)a;
}
}
q += 2 * NR;
}
}
}
template <int OP1, int OP2>
static void packb(const TE *B, int ldb, const TE *B2, const TE *B3, const TE *B4, int ldb2,
int kc, int nc, int pc, int jc) {
const int kn = kc >> 1;
int jr = 0;
const TE *b2 = (OP1 == OP_RAW) ? B : B2;
const TE *b3 = (OP2 == OP_RAW) ? B : B3;
const TE *b4 = (OP1 == OP_RAW) ? B3 : B4;
{
const int jrl = nc - (nc % (2 * NR));
for (int kp = 0; kp < kn; kp++) {
const int k = kp << 1;
const TE *s0 = B + (size_t)(pc + k) * ldb + jc;
const TE *s1 = B + (size_t)(pc + k + 1) * ldb + jc;
for (int jr = 0; jr < jrl; jr += 2 * NR) {
TE *q = Bpanel + (size_t)(jr / NR) * ((size_t)kc * NR + BPAD) + (size_t)kp * 2 * NR;
TE *q2 = q + (size_t)kc * NR + BPAD;
__m256i A0 = LDV(s0 + jr), A1 = LDV(s0 + jr + 16), B0 = LDV(s1 + jr), B1 = LDV(s1 + jr + 16);
if (OP1 != OP_RAW) {
const TE *t0 = b2 + (size_t)(pc + k) * ldb2 + jc + jr;
const TE *t1 = b2 + (size_t)(pc + k + 1) * ldb2 + jc + jr;
if (OP1 == OP_ADD) { A0 = _mm256_add_epi16(A0, LDV(t0)); A1 = _mm256_add_epi16(A1, LDV(t0 + 16));
B0 = _mm256_add_epi16(B0, LDV(t1)); B1 = _mm256_add_epi16(B1, LDV(t1 + 16)); }
else { A0 = _mm256_sub_epi16(A0, LDV(t0)); A1 = _mm256_sub_epi16(A1, LDV(t0 + 16));
B0 = _mm256_sub_epi16(B0, LDV(t1)); B1 = _mm256_sub_epi16(B1, LDV(t1 + 16)); }
}
if (OP2 != OP_RAW) {
const TE *u0 = b3 + (size_t)(pc + k) * ldb2 + jc + jr;
const TE *u1 = b3 + (size_t)(pc + k + 1) * ldb2 + jc + jr;
__m256i C0 = LDV(u0), C1 = LDV(u0 + 16), D0 = LDV(u1), D1 = LDV(u1 + 16);
if (OP1 != OP_RAW) {
const TE *v0 = b4 + (size_t)(pc + k) * ldb2 + jc + jr;
const TE *v1 = b4 + (size_t)(pc + k + 1) * ldb2 + jc + jr;
if (OP1 == OP_ADD) { C0 = _mm256_add_epi16(C0, LDV(v0)); C1 = _mm256_add_epi16(C1, LDV(v0 + 16));
D0 = _mm256_add_epi16(D0, LDV(v1)); D1 = _mm256_add_epi16(D1, LDV(v1 + 16)); }
else { C0 = _mm256_sub_epi16(C0, LDV(v0)); C1 = _mm256_sub_epi16(C1, LDV(v0 + 16));
D0 = _mm256_sub_epi16(D0, LDV(v1)); D1 = _mm256_sub_epi16(D1, LDV(v1 + 16)); }
}
if (OP2 == OP_ADD) { A0 = _mm256_add_epi16(A0, C0); A1 = _mm256_add_epi16(A1, C1);
B0 = _mm256_add_epi16(B0, D0); B1 = _mm256_add_epi16(B1, D1); }
else { A0 = _mm256_sub_epi16(A0, C0); A1 = _mm256_sub_epi16(A1, C1);
B0 = _mm256_sub_epi16(B0, D0); B1 = _mm256_sub_epi16(B1, D1); }
}
STV(q, _mm256_unpacklo_epi16(A0, B0));
STV(q + 16, _mm256_unpackhi_epi16(A0, B0));
STV(q2, _mm256_unpacklo_epi16(A1, B1));
STV(q2 + 16, _mm256_unpackhi_epi16(A1, B1));
}
}
if (jrl < nc) packb_tail(B, ldb, B2, B3, B4, ldb2, kc, nc, pc, jc, jrl, OP1, OP2);
return;
}
if (jr < nc) packb_tail(B, ldb, B2, B3, B4, ldb2, kc, nc, pc, jc, jr, OP1, OP2);
}
#pragma GCC push_options
#pragma GCC target("avx2")
#define KGRP4_B(a0,a1,a2,a3,b0,b1) \
"vmovdqa " #b0 "(%[b]), %%ymm14\n\t" \
"vmovdqa " #b1 "(%[b]), %%ymm15\n\t" \
"vpaddd %%ymm13, %6, %6\n\t" \
"vpaddd %%ymm12, %7, %7\n\t" \
KST(a0,"0","1") KST(a1,"2","3") KST(a2,"4","5") KST5D(a3)
#define KGRP4 KGRP4_B(0,512,1024,1536, 0,32) \
KGRP4_B(4,516,1028,1540, 64,96) \
KGRP4_B(8,520,1032,1544, 128,160) \
KGRP4_B(12,524,1036,1548, 192,224)
template <int OP>
static inline void micro4(long k, const TE *ap, const TE *bp, TE *C, int ldc, int rows) {
__m256i c00 = _mm256_setzero_si256(), c01 = c00, c10 = c00, c11 = c00;
__m256i c20 = c00, c21 = c00, c30 = c00, c31 = c00;
__asm__ volatile(
"vpxor %%ymm12, %%ymm12, %%ymm12\n\t"
"vpxor %%ymm13, %%ymm13, %%ymm13\n\t"
"1:\n\t" KGRP4
"addq $16, %[a]\n\t"
"addq $256, %[b]\n\t"
"decq %[k]\n\t"
"jnz 1b\n\t"
"vpaddd %%ymm13, %6, %6\n\t"
"vpaddd %%ymm12, %7, %7\n\t"
: "+x"(c00), "+x"(c01), "+x"(c10), "+x"(c11), "+x"(c20), "+x"(c21),
"+x"(c30), "+x"(c31), [a] "+r"(ap), [b] "+r"(bp), [k] "+r"(k)
:
: MICRO_CLOB);
ST_ROW(0, c00, c01) ST_ROW(1, c10, c11) ST_ROW(2, c20, c21) ST_ROW(3, c30, c31)
}
#pragma GCC pop_options
static void leaf(TE *C, int ldc, const TE *A, int lda, const TE *A2, int lda2,
const TE *A3, const TE *A4, int aop, int aop2,
const TE *B, int ldb, const TE *B2, int ldb2,
const TE *B3, const TE *B4, int bop, int bop2, int n) {
const int nbc = n - (n % NR);
for (int jc = 0; jc < nbc; jc += NC) {
const int nc = (nbc - jc < NC) ? (nbc - jc) : NC;
for (int pc = 0; pc < n; pc += KC) {
const int kc = (n - pc < KC) ? (n - pc) : KC;
{
if (bop2 == OP_RAW) {
if (bop == OP_RAW) packb<OP_RAW, OP_RAW>(B, ldb, 0, 0, 0, ldb2, kc, nc, pc, jc);
else if (bop == OP_ADD) packb<OP_ADD, OP_RAW>(B, ldb, B2, 0, 0, ldb2, kc, nc, pc, jc);
else packb<OP_SUB, OP_RAW>(B, ldb, B2, 0, 0, ldb2, kc, nc, pc, jc);
} else if (bop == OP_RAW) {
if (bop2 == OP_ADD) packb<OP_RAW, OP_ADD>(B, ldb, 0, B3, B3, ldb2, kc, nc, pc, jc);
else packb<OP_RAW, OP_SUB>(B, ldb, 0, B3, B3, ldb2, kc, nc, pc, jc);
} else if (bop == OP_ADD) {
if (bop2 == OP_ADD) packb<OP_ADD, OP_ADD>(B, ldb, B2, B3, B4, ldb2, kc, nc, pc, jc);
else packb<OP_ADD, OP_SUB>(B, ldb, B2, B3, B4, ldb2, kc, nc, pc, jc);
} else {
if (bop2 == OP_ADD) packb<OP_SUB, OP_ADD>(B, ldb, B2, B3, B4, ldb2, kc, nc, pc, jc);
else packb<OP_SUB, OP_SUB>(B, ldb, B2, B3, B4, ldb2, kc, nc, pc, jc);
}
}
for (int ic = 0; ic < n; ic += MR) {
const int rows = (n - ic < MR) ? (n - ic) : MR;
struct __attribute__((aligned(64))) { char pad[160u]; TE t[MR][TRS]; char tail[864u]; } g_trb;
if (n < 0) { g_trb.pad[0] = 0; g_trb.tail[0] = 0; }
#define tr (g_trb.t)
{
const TE *w1 = A + (size_t)ic * lda + pc;
const TE *w2 = A2 ? A2 + (size_t)ic * lda2 + pc : 0;
const TE *w3 = A3 ? A3 + (size_t)ic * lda2 + pc : 0;
const TE *w4 = A4 ? A4 + (size_t)ic * lda2 + pc : 0;
if (rows == MR && kc == KC) {
if (aop2 == OP_RAW) {
if (aop == OP_RAW) trb<OP_RAW, OP_RAW, MR, KC>(tr, rows, kc, w1, w1, w1, w1, lda, lda2);
else if (aop == OP_ADD) trb<OP_ADD, OP_RAW, MR, KC>(tr, rows, kc, w1, w2, w1, w1, lda, lda2);
else trb<OP_SUB, OP_RAW, MR, KC>(tr, rows, kc, w1, w2, w1, w1, lda, lda2);
} else if (aop == OP_RAW) {
if (aop2 == OP_ADD) trb<OP_RAW, OP_ADD, MR, KC>(tr, rows, kc, w1, w1, w3, w3, lda, lda2);
else trb<OP_RAW, OP_SUB, MR, KC>(tr, rows, kc, w1, w1, w3, w3, lda, lda2);
} else if (aop == OP_ADD) {
if (aop2 == OP_ADD) trb<OP_ADD, OP_ADD, MR, KC>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
else trb<OP_ADD, OP_SUB, MR, KC>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
} else {
if (aop2 == OP_ADD) trb<OP_SUB, OP_ADD, MR, KC>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
else trb<OP_SUB, OP_SUB, MR, KC>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
}
} else if (rows == MR) {
if (aop2 == OP_RAW) {
if (aop == OP_RAW) trb<OP_RAW, OP_RAW, MR, -1>(tr, rows, kc, w1, w1, w1, w1, lda, lda2);
else if (aop == OP_ADD) trb<OP_ADD, OP_RAW, MR, -1>(tr, rows, kc, w1, w2, w1, w1, lda, lda2);
else trb<OP_SUB, OP_RAW, MR, -1>(tr, rows, kc, w1, w2, w1, w1, lda, lda2);
} else if (aop == OP_RAW) {
if (aop2 == OP_ADD) trb<OP_RAW, OP_ADD, MR, -1>(tr, rows, kc, w1, w1, w3, w3, lda, lda2);
else trb<OP_RAW, OP_SUB, MR, -1>(tr, rows, kc, w1, w1, w3, w3, lda, lda2);
} else if (aop == OP_ADD) {
if (aop2 == OP_ADD) trb<OP_ADD, OP_ADD, MR, -1>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
else trb<OP_ADD, OP_SUB, MR, -1>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
} else {
if (aop2 == OP_ADD) trb<OP_SUB, OP_ADD, MR, -1>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
else trb<OP_SUB, OP_SUB, MR, -1>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
}
} else {
if (aop2 == OP_RAW) {
if (aop == OP_RAW) trb<OP_RAW, OP_RAW, -1, -1>(tr, rows, kc, w1, w1, w1, w1, lda, lda2);
else if (aop == OP_ADD) trb<OP_ADD, OP_RAW, -1, -1>(tr, rows, kc, w1, w2, w1, w1, lda, lda2);
else trb<OP_SUB, OP_RAW, -1, -1>(tr, rows, kc, w1, w2, w1, w1, lda, lda2);
} else if (aop == OP_RAW) {
if (aop2 == OP_ADD) trb<OP_RAW, OP_ADD, -1, -1>(tr, rows, kc, w1, w1, w3, w3, lda, lda2);
else trb<OP_RAW, OP_SUB, -1, -1>(tr, rows, kc, w1, w1, w3, w3, lda, lda2);
} else if (aop == OP_ADD) {
if (aop2 == OP_ADD) trb<OP_ADD, OP_ADD, -1, -1>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
else trb<OP_ADD, OP_SUB, -1, -1>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
} else {
if (aop2 == OP_ADD) trb<OP_SUB, OP_ADD, -1, -1>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
else trb<OP_SUB, OP_SUB, -1, -1>(tr, rows, kc, w1, w2, w3, w4, lda, lda2);
}
}
}
if (rows != MR && rows != 1)
for (int r = rows; r < 4; r++)
for (int k = 0; k < kc; k++) tr[r][k] = 0; /* micro4 reads tr rows 0..3 only */ /* micro4 reads tr rows 0..3 only */
const TE *trp = (const TE *)tr;
const TE *b0 = Bpanel;
TE *cp = C + (size_t)ic * ldc + jc;
for (int jr = 0; jr < nc; jr += NR) {
const int kn = kc >> 1;
const int kf = kn & ~(UNROLL - 1);
if (rows == MR) mM5<0,MR>(kf / UNROLL, trp, b0, cp + jr, ldc, rows);
else if (rows == 1) microR1<0,1>(kf / UNROLL, trp, b0, cp + jr, ldc, rows);
else micro4<0>(kf / UNROLL, trp, b0, cp + jr, ldc, rows);
b0 += (size_t)kc * NR + BPAD;
}
}
}
}
if (nbc < n) {
for (int i = 0; i < n; i++)
for (int j = nbc; j < n; j++) {
unsigned acc = 0;
for (int k = 0; k < n; k++)
acc += (unsigned)(unsigned short)A[(size_t)i * lda + k] * (unsigned)(unsigned short)B[(size_t)k * ldb + j];
C[(size_t)i * ldc + j] = (TE)acc;
}
}
}
static TE *g_pool;
static void *g_pbase; /* OF-KNOB: base actually handed to free() */
static void *g_cbase; /* OF-KNOB: g_cc base */
static void *g_rbase; /* OF-KNOB: g_rc base */
static size_t g_off;
static TE *g_scr;
static size_t g_soff;
static int g_topn = -1;
static void mk5(TE *AS, TE *BS, const TE *A, int lda, const TE *B, int ldb, int h,
int bipA, int bipB) {
const size_t q = (size_t)(h + PAD) * h;
const size_t qs = q + 1536;
const size_t qb = (size_t)((h >> 1) + PAD) * (h >> 1) + 2048;
const int bpw = (g_kbp >= 1 && h == 2048) || (g_kbp >= 2 && h == 1024);
const int hb = (h == 2048) ? 10 : 9;
const int hm = (h == 2048) ? 1023 : 511;
const size_t qsx = (h == 2048 || h == 1024) ? 4 * qb : qs;
const int hph = (h >> 1) + PAD;
const TE *a11, *a12, *a21, *a22, *b11, *b12, *b21, *b22;
if (bipA) { a11 = A; a12 = A + qs; a21 = A + 2 * qs; a22 = A + 3 * qs; lda = h + PAD; }
else { a11 = A; a12 = A + h; a21 = A + (size_t)h * lda; a22 = a21 + h; }
if (bipB) { b11 = B; b12 = B + qs; b21 = B + 2 * qs; b22 = B + 3 * qs; ldb = h + PAD; }
else { b11 = B; b12 = B + h; b21 = B + (size_t)h * ldb; b22 = b21 + h; }
TE *const *const cd = (h == g_ch) ? g_cc : g_rc;
const int carm = (g_copy && g_cok && h == g_ch && g_cc[0]) ||
(g_copy2 && g_cok && g_rh && h == g_rh && g_rc[0]);
const int copyA = carm && !bipA, copyB = carm && !bipB;
if (copyA || copyB) {
for (int i = 0; i < h; i++) {
const TE *r11 = a11 + (size_t)i * lda, *r12 = a12 + (size_t)i * lda;
const TE *r21 = a21 + (size_t)i * lda, *r22 = a22 + (size_t)i * lda;
TE *o0 = bpw ? (AS + ((size_t)((i >> hb) << 1)) * qb + (size_t)(i & hm) * (size_t)hph)
: (AS + (size_t)i * (h + PAD));
TE *o1 = o0 + qsx, *o2 = o1 + qsx, *o3 = o2 + qsx, *o4 = o3 + qsx;
TE *c11r = cd[0] + (size_t)i * (h + PAD), *c22r = cd[1] + (size_t)i * (h + PAD);
int j = 0;
#define MKAC(J) { const int jj = (J); \
__m256i q11 = LDV(r11 + jj), q12 = LDV(r12 + jj), q21 = LDV(r21 + jj), q22 = LDV(r22 + jj); \
__m256i w11 = LDV(r11 + jj + 16), w12 = LDV(r12 + jj + 16), w21 = LDV(r21 + jj + 16), w22 = LDV(r22 + jj + 16); \
_mm256_stream_si256((__m256i *)(o0 + jj), _mm256_add_epi16(q11, q22)); \
_mm256_stream_si256((__m256i *)(o0 + jj + 16), _mm256_add_epi16(w11, w22)); \
_mm256_stream_si256((__m256i *)(o1 + jj), _mm256_add_epi16(q21, q22)); \
_mm256_stream_si256((__m256i *)(o1 + jj + 16), _mm256_add_epi16(w21, w22)); \
_mm256_stream_si256((__m256i *)(o2 + jj), _mm256_add_epi16(q11, q12)); \
_mm256_stream_si256((__m256i *)(o2 + jj + 16), _mm256_add_epi16(w11, w12)); \
_mm256_stream_si256((__m256i *)(o3 + jj), _mm256_sub_epi16(q21, q11)); \
_mm256_stream_si256((__m256i *)(o3 + jj + 16), _mm256_sub_epi16(w21, w11)); \
_mm256_stream_si256((__m256i *)(o4 + jj), _mm256_sub_epi16(q12, q22)); \
_mm256_stream_si256((__m256i *)(o4 + jj + 16), _mm256_sub_epi16(w12, w22)); \
_mm256_stream_si256((__m256i *)(c11r + jj), q11); \
_mm256_stream_si256((__m256i *)(c11r + jj + 16), w11); \
_mm256_stream_si256((__m256i *)(c22r + jj), q22); \
_mm256_stream_si256((__m256i *)(c22r + jj + 16), w22); }
#define MKACS(J) { const int jj = (J); \
__m256i q11 = LDV(r11 + jj), q12 = LDV(r12 + jj), q21 = LDV(r21 + jj), q22 = LDV(r22 + jj); \
__m256i w11 = LDV(r11 + jj + 16), w12 = LDV(r12 + jj + 16), w21 = LDV(r21 + jj + 16), w22 = LDV(r22 + jj + 16); \
_mm256_stream_si256((__m256i *)(o0 + jj), _mm256_add_epi16(q11, q22)); \
_mm256_stream_si256((__m256i *)(o0 + jj + 16), _mm256_add_epi16(w11, w22)); \
_mm256_stream_si256((__m256i *)(o1 + jj), _mm256_add_epi16(q21, q22)); \
_mm256_stream_si256((__m256i *)(o1 + jj + 16), _mm256_add_epi16(w21, w22)); \
_mm256_stream_si256((__m256i *)(o2 + jj), _mm256_add_epi16(q11, q12)); \
_mm256_stream_si256((__m256i *)(o2 + jj + 16), _mm256_add_epi16(w11, w12)); \
_mm256_stream_si256((__m256i *)(o3 + jj), _mm256_sub_epi16(q21, q11)); \
_mm256_stream_si256((__m256i *)(o3 + jj + 16), _mm256_sub_epi16(w21, w11)); \
_mm256_stream_si256((__m256i *)(o4 + jj), _mm256_sub_epi16(q12, q22)); \
_mm256_stream_si256((__m256i *)(o4 + jj + 16), _mm256_sub_epi16(w12, w22)); }
if (copyA) {
if (bpw) { const int h1 = h >> 1;
for (; j + 32 <= h1; j += 32) { MKAC(j) }
o0 += qb - h1; o1 += qb - h1; o2 += qb - h1; o3 += qb - h1; o4 += qb - h1;
}
for (; j + 32 <= h; j += 32) { MKAC(j) }
} else {
if (bpw) { const int h1 = h >> 1;
for (; j + 32 <= h1; j += 32) { MKACS(j) }
o0 += qb - h1; o1 += qb - h1; o2 += qb - h1; o3 += qb - h1; o4 += qb - h1;
}
for (; j + 32 <= h; j += 32) { MKACS(j) }
}
#undef MKACS
#undef MKAC
}
for (int i = 0; i < h; i++) {
const TE *s11 = b11 + (size_t)i * ldb, *s12 = b12 + (size_t)i * ldb;
const TE *s21 = b21 + (size_t)i * ldb, *s22 = b22 + (size_t)i * ldb;
TE *p0 = bpw ? (BS + ((size_t)((i >> hb) << 1)) * qb + (size_t)(i & hm) * (size_t)hph)
: (BS + (size_t)i * (h + PAD));
TE *p1 = p0 + qsx, *p2 = p1 + qsx, *p3 = p2 + qsx, *p4 = p3 + qsx;
TE *c33r = cd[2] + (size_t)i * (h + PAD), *c44r = cd[3] + (size_t)i * (h + PAD);
int j = 0;
#define MKBC(J) { const int jj = (J); \
__m256i q11 = LDV(s11 + jj), q12 = LDV(s12 + jj), q21 = LDV(s21 + jj), q22 = LDV(s22 + jj); \
__m256i w11 = LDV(s11 + jj + 16), w12 = LDV(s12 + jj + 16), w21 = LDV(s21 + jj + 16), w22 = LDV(s22 + jj + 16); \
_mm256_stream_si256((__m256i *)(p0 + jj), _mm256_add_epi16(q11, q22)); \
_mm256_stream_si256((__m256i *)(p0 + jj + 16), _mm256_add_epi16(w11, w22)); \
_mm256_stream_si256((__m256i *)(p1 + jj), _mm256_sub_epi16(q12, q22)); \
_mm256_stream_si256((__m256i *)(p1 + jj + 16), _mm256_sub_epi16(w12, w22)); \
_mm256_stream_si256((__m256i *)(p2 + jj), _mm256_sub_epi16(q21, q11)); \
_mm256_stream_si256((__m256i *)(p2 + jj + 16), _mm256_sub_epi16(w21, w11)); \
_mm256_stream_si256((__m256i *)(p3 + jj), _mm256_add_epi16(q11, q12)); \
_mm256_stream_si256((__m256i *)(p3 + jj + 16), _mm256_add_epi16(w11, w12)); \
_mm256_stream_si256((__m256i *)(p4 + jj), _mm256_add_epi16(q21, q22)); \
_mm256_stream_si256((__m256i *)(p4 + jj + 16), _mm256_add_epi16(w21, w22)); \
_mm256_stream_si256((__m256i *)(c33r + jj), q11); \
_mm256_stream_si256((__m256i *)(c33r + jj + 16), w11); \
_mm256_stream_si256((__m256i *)(c44r + jj), q22); \
_mm256_stream_si256((__m256i *)(c44r + jj + 16), w22); }
#define MKBCS(J) { const int jj = (J); \
__m256i q11 = LDV(s11 + jj), q12 = LDV(s12 + jj), q21 = LDV(s21 + jj), q22 = LDV(s22 + jj); \
__m256i w11 = LDV(s11 + jj + 16), w12 = LDV(s12 + jj + 16), w21 = LDV(s21 + jj + 16), w22 = LDV(s22 + jj + 16); \
_mm256_stream_si256((__m256i *)(p0 + jj), _mm256_add_epi16(q11, q22)); \
_mm256_stream_si256((__m256i *)(p0 + jj + 16), _mm256_add_epi16(w11, w22)); \
_mm256_stream_si256((__m256i *)(p1 + jj), _mm256_sub_epi16(q12, q22)); \
_mm256_stream_si256((__m256i *)(p1 + jj + 16), _mm256_sub_epi16(w12, w22)); \
_mm256_stream_si256((__m256i *)(p2 + jj), _mm256_sub_epi16(q21, q11)); \
_mm256_stream_si256((__m256i *)(p2 + jj + 16), _mm256_sub_epi16(w21, w11)); \
_mm256_stream_si256((__m256i *)(p3 + jj), _mm256_add_epi16(q11, q12)); \
_mm256_stream_si256((__m256i *)(p3 + jj + 16), _mm256_add_epi16(w11, w12)); \
_mm256_stream_si256((__m256i *)(p4 + jj), _mm256_add_epi16(q21, q22)); \
_mm256_stream_si256((__m256i *)(p4 + jj + 16), _mm256_add_epi16(w21, w22)); }
if (copyB) {
if (bpw) { const int h1 = h >> 1;
for (; j + 32 <= h1; j += 32) { MKBC(j) }
p0 += qb - h1; p1 += qb - h1; p2 += qb - h1; p3 += qb - h1; p4 += qb - h1;
}
for (; j + 32 <= h; j += 32) { MKBC(j) }
} else {
if (bpw) { const int h1 = h >> 1;
for (; j + 32 <= h1; j += 32) { MKBCS(j) }
p0 += qb - h1; p1 += qb - h1; p2 += qb - h1; p3 += qb - h1; p4 += qb - h1;
}
for (; j + 32 <= h; j += 32) { MKBCS(j) }
}
#undef MKBCS
#undef MKBC
}
if (g_lv) {
for (int i = 0; i < h; i += 64) {
cd[0][(size_t)i * (h + PAD)] ^= 1;
cd[1][(size_t)i * (h + PAD)] ^= 1;
cd[2][(size_t)i * (h + PAD)] ^= 1;
cd[3][(size_t)i * (h + PAD)] ^= 1;
}
}
return;
}
for (int i = 0; i < h; i++) {
const TE *r11 = a11 + (size_t)i * lda, *r12 = a12 + (size_t)i * lda;
const TE *r21 = a21 + (size_t)i * lda, *r22 = a22 + (size_t)i * lda;
TE *o0 = bpw ? (AS + ((size_t)((i >> hb) << 1)) * qb + (size_t)(i & hm) * (size_t)hph)
: (AS + (size_t)i * (h + PAD));
TE *o1 = o0 + qsx, *o2 = o1 + qsx, *o3 = o2 + qsx, *o4 = o3 + qsx;
int j = 0;
#define MKA(J) { const int jj = (J); \
__m256i v11 = LDV(r11 + jj), v12 = LDV(r12 + jj), v21 = LDV(r21 + jj), v22 = LDV(r22 + jj); \
_mm256_stream_si256((__m256i *)(o0 + jj), _mm256_add_epi16(v11, v22)); \
_mm256_stream_si256((__m256i *)(o1 + jj), _mm256_add_epi16(v21, v22)); \
_mm256_stream_si256((__m256i *)(o2 + jj), _mm256_add_epi16(v11, v12)); \
_mm256_stream_si256((__m256i *)(o3 + jj), _mm256_sub_epi16(v21, v11)); \
_mm256_stream_si256((__m256i *)(o4 + jj), _mm256_sub_epi16(v12, v22)); }
#if MK5O
#define MKA2(J) { const int jj = (J); \
__m256i q11 = LDV(r11 + jj), q12 = LDV(r12 + jj), q21 = LDV(r21 + jj), q22 = LDV(r22 + jj); \
__m256i w11 = LDV(r11 + jj + 16), w12 = LDV(r12 + jj + 16), w21 = LDV(r21 + jj + 16), w22 = LDV(r22 + jj + 16); \
_mm256_stream_si256((__m256i *)(o0 + jj), _mm256_add_epi16(q11, q22)); \
_mm256_stream_si256((__m256i *)(o0 + jj + 16), _mm256_add_epi16(w11, w22)); \
_mm256_stream_si256((__m256i *)(o1 + jj), _mm256_add_epi16(q21, q22)); \
_mm256_stream_si256((__m256i *)(o1 + jj + 16), _mm256_add_epi16(w21, w22)); \
_mm256_stream_si256((__m256i *)(o2 + jj), _mm256_add_epi16(q11, q12)); \
_mm256_stream_si256((__m256i *)(o2 + jj + 16), _mm256_add_epi16(w11, w12)); \
_mm256_stream_si256((__m256i *)(o3 + jj), _mm256_sub_epi16(q21, q11)); \
_mm256_stream_si256((__m256i *)(o3 + jj + 16), _mm256_sub_epi16(w21, w11)); \
_mm256_stream_si256((__m256i *)(o4 + jj), _mm256_sub_epi16(q12, q22)); \
_mm256_stream_si256((__m256i *)(o4 + jj + 16), _mm256_sub_epi16(w12, w22)); }
if (bpw) { const int h1 = h >> 1;
for (; j + 32 <= h1; j += 32) { MKA2(j) }
o0 += qb - h1; o1 += qb - h1; o2 += qb - h1; o3 += qb - h1; o4 += qb - h1;
}
for (; j + 32 <= h; j += 32) { MKA2(j) }
#else
for (; j + 32 <= h; j += 32) { MKA(j) MKA(j + 16) }
#endif
for (; j + 16 <= h; j += 16) { MKA(j) }
#undef MKA
for (; j < h; j++) {
TE x11 = r11[j], x12 = r12[j], x21 = r21[j], x22 = r22[j];
o0[j] = x11 + x22; o1[j] = x21 + x22; o2[j] = x11 + x12;
o3[j] = x21 - x11; o4[j] = x12 - x22;
}
}
for (int i = 0; i < h; i++) {
const TE *s11 = b11 + (size_t)i * ldb, *s12 = b12 + (size_t)i * ldb;
const TE *s21 = b21 + (size_t)i * ldb, *s22 = b22 + (size_t)i * ldb;
TE *p0 = bpw ? (BS + ((size_t)((i >> hb) << 1)) * qb + (size_t)(i & hm) * (size_t)hph)
: (BS + (size_t)i * (h + PAD));
TE *p1 = p0 + qsx, *p2 = p1 + qsx, *p3 = p2 + qsx, *p4 = p3 + qsx;
int j = 0;
#define MKB(J) { const int jj = (J); \
__m256i w11 = LDV(s11 + jj), w12 = LDV(s12 + jj), w21 = LDV(s21 + jj), w22 = LDV(s22 + jj); \
_mm256_stream_si256((__m256i *)(p0 + jj), _mm256_add_epi16(w11, w22)); \
_mm256_stream_si256((__m256i *)(p1 + jj), _mm256_sub_epi16(w12, w22)); \
_mm256_stream_si256((__m256i *)(p2 + jj), _mm256_sub_epi16(w21, w11)); \
_mm256_stream_si256((__m256i *)(p3 + jj), _mm256_add_epi16(w11, w12)); \
_mm256_stream_si256((__m256i *)(p4 + jj), _mm256_add_epi16(w21, w22)); }
#if MK5O
#define MKB2(J) { const int jj = (J); \
__m256i q11 = LDV(s11 + jj), q12 = LDV(s12 + jj), q21 = LDV(s21 + jj), q22 = LDV(s22 + jj); \
__m256i w11 = LDV(s11 + jj + 16), w12 = LDV(s12 + jj + 16), w21 = LDV(s21 + jj + 16), w22 = LDV(s22 + jj + 16); \
_mm256_stream_si256((__m256i *)(p0 + jj), _mm256_add_epi16(q11, q22)); \
_mm256_stream_si256((__m256i *)(p0 + jj + 16), _mm256_add_epi16(w11, w22)); \
_mm256_stream_si256((__m256i *)(p1 + jj), _mm256_sub_epi16(q12, q22)); \
_mm256_stream_si256((__m256i *)(p1 + jj + 16), _mm256_sub_epi16(w12, w22)); \
_mm256_stream_si256((__m256i *)(p2 + jj), _mm256_sub_epi16(q21, q11)); \
_mm256_stream_si256((__m256i *)(p2 + jj + 16), _mm256_sub_epi16(w21, w11)); \
_mm256_stream_si256((__m256i *)(p3 + jj), _mm256_add_epi16(q11, q12)); \
_mm256_stream_si256((__m256i *)(p3 + jj + 16), _mm256_add_epi16(w11, w12)); \
_mm256_stream_si256((__m256i *)(p4 + jj), _mm256_add_epi16(q21, q22)); \
_mm256_stream_si256((__m256i *)(p4 + jj + 16), _mm256_add_epi16(w21, w22)); }
if (bpw) { const int h1 = h >> 1;
for (; j + 32 <= h1; j += 32) { MKB2(j) }
p0 += qb - h1; p1 += qb - h1; p2 += qb - h1; p3 += qb - h1; p4 += qb - h1;
}
for (; j + 32 <= h; j += 32) { MKB2(j) }
#else
for (; j + 32 <= h; j += 32) { MKB(j) MKB(j + 16) }
#endif
for (; j + 16 <= h; j += 16) { MKB(j) }
#undef MKB
for (; j < h; j++) {
TE y11 = s11[j], y12 = s12[j], y21 = s21[j], y22 = s22[j];
p0[j] = y11 + y22; p1[j] = y12 - y22; p2[j] = y21 - y11;
p3[j] = y11 + y12; p4[j] = y21 + y22;
}
}
}
static void combine(TE *C, int ldc, const TE *M1, const TE *M2, const TE *M3,
const TE *M4, const TE *M5, const TE *M6, const TE *M7, int h) {
volatile char g_cpad[128] __attribute__((aligned(64))); if (ldc < 0) g_cpad[0] = 0;
const size_t q = (size_t)(h + PAD) * h;
for (int i = 0; i < h; i++) {
const TE *m1 = M1 + (size_t)i * (h + PAD), *m2 = M2 + (size_t)i * (h + PAD), *m3 = M3 + (size_t)i * (h + PAD);
const TE *m4 = M4 + (size_t)i * (h + PAD), *m5 = M5 + (size_t)i * (h + PAD), *m6 = M6 + (size_t)i * (h + PAD);
const TE *m7 = M7 + (size_t)i * (h + PAD);
TE *c11 = C + (size_t)i * ldc, *c12 = c11 + h;
TE *c21 = c11 + (size_t)h * ldc, *c22 = c21 + h;
int j = 0;
#define CMB(J) { const int jj = (J); \
__m256i a1 = LDV(m1 + jj), a2 = LDV(m2 + jj), a3 = LDV(m3 + jj), a4 = LDV(m4 + jj); \
__m256i a5 = LDV(m5 + jj), a6 = LDV(m6 + jj), a7 = LDV(m7 + jj); \
_mm256_stream_si256((__m256i*)(c11 + jj), _mm256_add_epi16(_mm256_sub_epi16(_mm256_add_epi16(a1, a4), a5), a7)); \
_mm256_stream_si256((__m256i*)(c12 + jj), _mm256_add_epi16(a3, a5)); \
_mm256_stream_si256((__m256i*)(c21 + jj), _mm256_add_epi16(a2, a4)); \
_mm256_stream_si256((__m256i*)(c22 + jj), _mm256_add_epi16(_mm256_sub_epi16(_mm256_add_epi16(a1, a3), a2), a6)); }
#if ORD
for (; j + 32 <= h; j += 32) {
const int jj = j;
__m256i u1 = LDV(m1 + jj), u4 = LDV(m4 + jj), u5 = LDV(m5 + jj), u7 = LDV(m7 + jj);
__m256i v1 = LDV(m1 + jj + 16), v4 = LDV(m4 + jj + 16), v5 = LDV(m5 + jj + 16), v7 = LDV(m7 + jj + 16);
_mm256_stream_si256((__m256i *)(c11 + jj), _mm256_add_epi16(_mm256_sub_epi16(_mm256_add_epi16(u1, u4), u5), u7));
_mm256_stream_si256((__m256i *)(c11 + jj + 16), _mm256_add_epi16(_mm256_sub_epi16(_mm256_add_epi16(v1, v4), v5), v7));
__m256i u3 = LDV(m3 + jj), v3 = LDV(m3 + jj + 16);
_mm256_stream_si256((__m256i *)(c12 + jj), _mm256_add_epi16(u3, u5));
_mm256_stream_si256((__m256i *)(c12 + jj + 16), _mm256_add_epi16(v3, v5));
__m256i u2 = LDV(m2 + jj), v2 = LDV(m2 + jj + 16);
_mm256_stream_si256((__m256i *)(c21 + jj), _mm256_add_epi16(u2, u4));
_mm256_stream_si256((__m256i *)(c21 + jj + 16), _mm256_add_epi16(v2, v4));
__m256i u6 = LDV(m6 + jj), v6 = LDV(m6 + jj + 16);
_mm256_stream_si256((__m256i *)(c22 + jj), _mm256_add_epi16(_mm256_sub_epi16(_mm256_add_epi16(u1, u3), u2), u6));
_mm256_stream_si256((__m256i *)(c22 + jj + 16), _mm256_add_epi16(_mm256_sub_epi16(_mm256_add_epi16(v1, v3), v2), v6));
}
#else
for (; j + 32 <= h; j += 32) { CMB(j) CMB(j + 16) }
#endif
for (; j + 16 <= h; j += 16) { CMB(j) }
#undef CMB
for (; j < h; j++) {
c11[j] = m1[j] + m4[j] - m5[j] + m7[j];
c12[j] = m3[j] + m5[j];
c21[j] = m2[j] + m4[j];
c22[j] = m1[j] - m2[j] + m3[j] + m6[j];
}
}
(void)q;
}
static void mm(TE *C, int ldc, const TE *A, int lda, const TE *A2, int lda2, int aop,
const TE *B, int ldb, const TE *B2, int ldb2, int bop, int n) {
if (n <= CUT) {
leaf(C, ldc, A, lda, A2, lda2, 0, 0, aop, OP_RAW, B, ldb, B2, ldb2, 0, 0, bop, OP_RAW, n);
return;
}
const int h = n >> 1;
volatile TE g_mpad[16]; if (ldc < 0) g_mpad[0] = 0;
const int hp = h + PAD;
const size_t q = (size_t)hp * h;
const size_t qs = q + 1536;
const size_t qb = (size_t)((h >> 1) + PAD) * (h >> 1) + 2048;
const size_t qsx = (h == 2048 || h == 1024) ? 4 * qb : qs;
const int lds = (((h == 2048 && g_kbp >= 1) || (h == 1024 && g_kbp >= 2)) ? (int)((h >> 1) + PAD) : hp);
const int bipA = (g_kbp && h == 1024 && lda == hp);
const int bipB = (g_kbp && h == 1024 && ldb == hp);
const int bip4a = (g_kbp >= 2 && n == 1024 && lda == hp);
const int bip4b = (g_kbp >= 2 && n == 1024 && ldb == hp);
const int ismk5 = (n != 2 * CUT) && (n != 4 * CUT);
const int deep = (n < g_topn);
size_t *const op = deep ? &g_soff : &g_off;
TE *AS = (deep ? g_scr : g_pool) + *op;
TE *BS = AS + 5 * qsx;
TE *M = ismk5 ? BS + 5 * qsx : AS;
const size_t alloc = (ismk5 ? 10 * qsx : 7 * qs);
*op += alloc;
const TE *A11, *A12, *A21, *A22, *B11, *B12, *B21, *B22;
if (bipA) { A11 = A; A12 = A + qs; A21 = A + 2 * qs; A22 = A + 3 * qs; }
else { A11 = A; A12 = A + h; A21 = A + (size_t)h * lda; A22 = A21 + h; }
if (bipB) { B11 = B; B12 = B + qs; B21 = B + 2 * qs; B22 = B + 3 * qs; }
else { B11 = B; B12 = B + h; B21 = B + (size_t)h * ldb; B22 = B21 + h; }
if (n == 2 * CUT) {
const TE *C11 = A2, *C21 = A2 ? A2 + (size_t)h * lda2 : 0, *C12 = A2 ? A2 + h : 0;
const TE *C22 = C21 ? C21 + h : 0;
const TE *D11 = B2, *D21 = B2 ? B2 + (size_t)h * ldb2 : 0, *D12 = B2 ? B2 + h : 0;
const TE *D22 = D21 ? D21 + h : 0;
leaf(M + 5 * qs, hp, A21, lda, A11, lda, C21, C11, OP_SUB, aop, B11, ldb, B12, ldb, D11, D12, OP_ADD, bop, h);
leaf(M + 2 * qs, hp, A11, lda, 0, lda, C11, 0, OP_RAW, aop, B12, ldb, B22, ldb, D12, D22, OP_SUB, bop, h);
leaf(M + 1 * qs, hp, A21, lda, A22, lda, C21, C22, OP_ADD, aop, B11, ldb, 0, ldb, D11, 0, OP_RAW, bop, h);
leaf(M + 0 * qs, hp, A11, lda, A22, lda, C11, C22, OP_ADD, aop, B11, ldb, B22, ldb, D11, D22, OP_ADD, bop, h);
leaf(M + 3 * qs, hp, A22, lda, 0, lda, C22, 0, OP_RAW, aop, B21, ldb, B11, ldb, D21, D11, OP_SUB, bop, h);
leaf(M + 6 * qs, hp, A12, lda, A22, lda, C12, C22, OP_SUB, aop, B21, ldb, B22, ldb, D21, D22, OP_ADD, bop, h);
leaf(M + 4 * qs, hp, A11, lda, A12, lda, C11, C12, OP_ADD, aop, B22, ldb, 0, ldb, D22, 0, OP_RAW, bop, h);
combine(C, ldc, M, M + qs, M + 2 * qs, M + 3 * qs, M + 4 * qs, M + 5 * qs, M + 6 * qs, h);
*op -= alloc;
return;
}
if (n == 4 * CUT) {
const TE *ab[4], *bb[4];
if (bip4a) { ab[0] = A; ab[1] = A + qs; ab[2] = A + 2 * qs; ab[3] = A + 3 * qs; lda = hp; }
else { ab[0] = A; ab[1] = A + h; ab[2] = A + (size_t)h * lda; ab[3] = A + (size_t)h * lda + h; }
if (bip4b) { bb[0] = B; bb[1] = B + qs; bb[2] = B + 2 * qs; bb[3] = B + 3 * qs; ldb = hp; }
else { bb[0] = B; bb[1] = B + h; bb[2] = B + (size_t)h * ldb; bb[3] = B + (size_t)h * ldb + h; }
static const signed char CSQ[7] = {5,1,0,2,4,6,3};
static const signed char PA1[7] = {0, 2, 0, 3, 0, 2, 1};
static const signed char PA2[7] = {3, 3, -1, -1, 1, 0, 3};
static const signed char PAO[7] = {OP_ADD, OP_ADD, OP_RAW, OP_RAW, OP_ADD, OP_SUB, OP_SUB};
static const signed char PB1[7] = {0, 0, 1, 2, 3, 0, 2};
static const signed char PB2[7] = {3, -1, 3, 0, -1, 1, 3};
static const signed char PBO[7] = {OP_ADD, OP_RAW, OP_SUB, OP_SUB, OP_RAW, OP_ADD, OP_ADD};
for (int cc = 0; cc < 7; cc++) {
const int c = (int)CSQ[cc];
const TE *ap1 = ab[PA1[c]], *ap2 = (PA2[c] >= 0) ? ab[PA2[c]] : 0;
const TE *bp1 = bb[PB1[c]], *bp2 = (PB2[c] >= 0) ? bb[PB2[c]] : 0;
mm(M + (size_t)c * qs, hp, ap1, lda, ap2, lda, PAO[c], bp1, ldb, bp2, ldb, PBO[c], h);
}
combine(C, ldc, M, M + qs, M + 2 * qs, M + 3 * qs, M + 4 * qs, M + 5 * qs, M + 6 * qs, h);
*op -= alloc;
return;
}
MK5B();
mk5(AS, BS, A, lda, B, ldb, h, bipA, bipB);
MK5E(h);
const int csel = (g_cok && (g_wmode == 3 || (g_wmode == 1 && h == g_ch) ||
(g_wmode == 2 && h != g_ch)));
const int cpsel2 = (g_cok && (h == g_ch || h == g_rh) &&
((h == g_ch) ? g_use : g_use2));
TE *const *const cc = (h == g_ch) ? g_cc : g_rc;
const int sA11 = csel && g_wA11, sA22 = csel && g_wA22;
const int sB11 = csel && g_wB11, sB22 = csel && g_wB22;
const TE *cpA11 = sA11 ? (AS + 3 * qsx) : (bipA ? A11 : (cpsel2 ? cc[0] : A11));
const TE *cpA22 = sA22 ? (AS + 2 * qsx) : (bipA ? A22 : (cpsel2 ? cc[1] : A22));
const TE *cpB11 = sB11 ? (BS + 0 * qsx) : (bipB ? B11 : (cpsel2 ? cc[2] : B11));
const TE *cpB22 = sB22 ? (BS + 4 * qsx) : (bipB ? B22 : (cpsel2 ? cc[3] : B22));
const int clA11 = (cpsel2 || sA11 || bipA) ? hp : lda;
const int clA22 = (cpsel2 || sA22 || bipA) ? hp : lda;
const int clB11 = (cpsel2 || sB11 || bipB) ? hp : ldb;
const int clB22 = (cpsel2 || sB22 || bipB) ? hp : ldb;
mm(AS + 1 * qsx, hp, AS + 1 * qsx, lds, 0, lds, OP_RAW, cpB11, clB11, 0, clB11, OP_RAW, h);
mm(BS + 2 * qsx, hp, cpA22, clA22, 0, clA22, OP_RAW, BS + 2 * qsx, lds, 0, lds, OP_RAW, h);
mm(AS + 3 * qsx, hp, AS + 3 * qsx, lds, 0, lds, OP_RAW, BS + 3 * qsx, lds, 0, lds, OP_RAW, h);
mm(AS + 0 * qsx, hp, AS + 0 * qsx, lds, 0, lds, OP_RAW, BS + 0 * qsx, lds, 0, lds, OP_RAW, h);
mm(BS + 1 * qsx, hp, cpA11, clA11, 0, clA11, OP_RAW, BS + 1 * qsx, lds, 0, lds, OP_RAW, h);
mm(AS + 2 * qsx, hp, AS + 2 * qsx, lds, 0, lds, OP_RAW, cpB22, clB22, 0, clB22, OP_RAW, h);
mm(AS + 4 * qsx, hp, AS + 4 * qsx, lds, 0, lds, OP_RAW, BS + 4 * qsx, lds, 0, lds, OP_RAW, h);
combine(C, ldc, AS + 0 * qsx, AS + 1 * qsx, BS + 1 * qsx, BS + 2 * qsx, AS + 2 * qsx, AS + 3 * qsx, AS + 4 * qsx, h);
*op -= alloc;
}
static void cpy(TE *dst, int ldd, const TE *src, int lds, int n) {
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j += 16)
_mm256_store_si256((__m256i *)(dst + (size_t)i * ldd + j),
_mm256_loadu_si256((const __m256i *)(src + (size_t)i * lds + j)));
}
void matrix_multiply(int n, const short *A, const short *B, short *C) {
if ((n % 16) || ((unsigned long)C & 31) || ((unsigned long)A & 31) || ((unsigned long)B & 31)) {
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++) {
unsigned acc = 0;
for (int k = 0; k < n; k++) acc += (unsigned)(unsigned short)A[(size_t)i * n + k] *
(unsigned)(unsigned short)B[(size_t)k * n + j];
C[(size_t)i * n + j] = (TE)acc;
}
return;
}
volatile int g_mpad2[2];
if (C == 0) g_mpad2[0] = 0;
const int ld = n;
size_t need = 64, need_scr = 0;
for (int m = n; m > CUT; m >>= 1) {
int hh = m >> 1;
size_t sl = (size_t)(hh + PAD) * hh + 2048;
if (hh == 2048 || hh == 1024) sl = 4 * ((size_t)((hh >> 1) + PAD) * (hh >> 1) + 2048);
size_t a = ((m == 2 * CUT || m == 4 * CUT) ? 7 : 10) * sl;
need += a;
if (m != n) need_scr += a;
}
const unsigned long ca = (unsigned long)(const void *)C;
const unsigned long aa = (unsigned long)(const void *)A, bb = (unsigned long)(const void *)B;
const unsigned long cb = (size_t)n * (size_t)n * sizeof(TE);
int split = (need_scr > 0) &&
(need_scr + 4096 <= (size_t)n * (size_t)n) &&
(ca + cb <= aa || aa + cb <= ca) &&
(ca + cb <= bb || bb + cb <= ca);
if (g_k.sp >= 0) split = g_k.sp; /* OF-KNOB: force the pointer-valued branch */
if (need > 64) {
void *p = 0;
const size_t _po = (size_t)(unsigned)g_k.off * 64; /* OF-KNOB: pool shift in bytes */
if (posix_memalign(&p, 4096, (split ? need - need_scr : need) * sizeof(TE) + _po + 4096)) return;
g_pbase = p;
g_pool = (TE *)((char *)p + _po);
g_off = 0;
g_topn = split ? n : -1;
{ const size_t _mx = (size_t)n * (size_t)n - need_scr;
const size_t _of = (size_t)3481664 < _mx ? (size_t)3481664 : _mx;
g_scr = split ? ((TE *)C + _of) : g_pool; }
g_soff = 0;
void *pc = 0;
const size_t htop = (size_t)n >> 1;
const size_t qc = ((htop + PAD) * htop) + 2048;
g_ch = (int)htop;
g_cok = (htop % 32 == 0);
const size_t _co = (size_t)(unsigned)g_k.offc * 64; /* OF-KNOB: cc shift in bytes */
if (g_alloc && g_cok && posix_memalign(&pc, 4096, 4 * qc * sizeof(TE) + _co + 4096) == 0) {
g_cbase = pc; pc = (char *)pc + _co;
g_cc[0] = (TE *)pc; g_cc[1] = g_cc[0] + qc; g_cc[2] = g_cc[1] + qc; g_cc[3] = g_cc[2] + qc;
} else { g_cc[0] = g_cc[1] = g_cc[2] = g_cc[3] = 0; g_cok = 0; }
void *pr = 0;
const size_t hr = htop >> 1;
const size_t qr = ((hr + PAD) * hr) + 2048;
g_rh = (int)hr;
const size_t _ro = (size_t)(unsigned)g_k.offc * 64; /* OF-KNOB: rc shift in bytes */
if (g_alloc && g_rh > 0 && g_rh % 32 == 0 && posix_memalign(&pr, 4096, 4 * qr * sizeof(TE) + _ro + 4096) == 0) {
g_rbase = pr; pr = (char *)pr + _ro;
g_rc[0] = (TE *)pr; g_rc[1] = g_rc[0] + qr; g_rc[2] = g_rc[1] + qr; g_rc[3] = g_rc[2] + qr;
} else { g_rc[0] = g_rc[1] = g_rc[2] = g_rc[3] = 0; }
{ if (g_k.dbg) printf("IN cok=%d ch=%d rh=%d cc=%p rc=%p pool=%p scr=%p split=%d"
" need=%zu scr=%zu\n", g_cok, g_ch, g_rh, (void*)g_cc[0],
(void*)g_rc[0], (void*)g_pool, (void*)g_scr, split, need, need_scr); }
mm(C, n, A, ld, 0, ld, OP_RAW, B, ld, 0, ld, OP_RAW, n);
free(g_rbase); g_rbase = 0; g_rc[0] = g_rc[1] = g_rc[2] = g_rc[3] = 0; g_rh = 0;
free(g_cbase); g_cbase = 0; g_cc[0] = g_cc[1] = g_cc[2] = g_cc[3] = 0; g_ch = 0; g_cok = 0;
free(g_pbase); g_pbase = 0;
} else {
leaf(C, n, A, ld, 0, ld, 0, 0, OP_RAW, OP_RAW, B, ld, 0, ld, 0, 0, OP_RAW, OP_RAW, n);
}
}
#pragma GCC pop_options
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 12.173 ms | 5 MB + 916 KB | Accepted | Score: 100 | 显示更多 |