/* FILE-SCOPE FLAGS: judge compiles with g++ -O2 -static -nostdlib -U_FORTIFY_SOURCE and NO
* -march, so the baseline is x86-64 = SSE2 ONLY. Widen the target for EVERY function in
* the file (including all packing / C-traffic code) and lift -O2 to -O3 + unroll. */
#pragma GCC optimize("O3")
#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt")
/* int32 GEMM exact mod 2^32, AVX2. Three-madd vpmaddwd mikro (KC=128 baked in),
* wrapped in EXACT Strassen over Z/2^32 (7/8 of the multiply work per level).
* GENERATED by gen_strassen.py -- do not edit by hand.
*/
#include <immintrin.h>
#include <string.h>
#include <stdint.h>
#include <stdlib.h>
#define MR 8
#define NCG 1
#define KC 128
#define MC 8
#define NC 128
#define NPAIR (KC / 2)
#define APSTR (NPAIR * 16 + 16)
#define NR (4 * 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))); /* [j4][m][lo16B|hi16B] 132 B dense */
#define NTSEL 0
/* MIKRO-BEGIN */
#define AVX2 __attribute__((target("avx2")))
static __m256i LA[128] __attribute__((aligned(32)));
AVX2 __attribute__((always_inline)) static inline void mikro(int npair, const int32_t *ap, const int32_t *bp,
int32_t *c, int ldc, int first) {
(void)npair;
ap += 32; bp += 16; /* bias so every disp fits in disp8 */
__m256i q0, q1, q2, q3;
__asm__ volatile(
"vpxor %%xmm15, %%xmm15, %%xmm15\n\t"
"vmovdqa %%ymm15, %%ymm0\n\t"
"vmovdqa %%ymm15, %%ymm4\n\t"
"vmovdqa %%ymm15, %%ymm1\n\t"
"vmovdqa %%ymm15, %%ymm5\n\t"
"vmovdqa %%ymm15, %%ymm2\n\t"
"vmovdqa %%ymm15, %%ymm6\n\t"
"vmovdqa %%ymm15, %%ymm3\n\t"
"vmovdqa %%ymm15, %%ymm7\n\t"
"test %[cnt], %[cnt]\n\t"
"jle 2f\n\t"
"1:\n\t"
"prefetcht0 192(%[bp])\n\t"
"vbroadcasti128 -64(%[bp]), %%ymm8\n\t"
"vbroadcasti128 -48(%[bp]), %%ymm9\n\t"
"vpmaddwd -128(%[ap]), %%ymm8, %%ymm10\n\t"
"vpmaddwd 4032(%[ap]), %%ymm8, %%ymm11\n\t"
"vpmaddwd 8192(%[ap]), %%ymm8, %%ymm12\n\t"
"vpmaddwd 12352(%[ap]), %%ymm8, %%ymm13\n\t"
"vpmaddwd -128(%[ap]), %%ymm9, %%ymm14\n\t"
"vpmaddwd 4032(%[ap]), %%ymm9, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm0, %%ymm0\n\t"
"vpmaddwd 8192(%[ap]), %%ymm9, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm1, %%ymm1\n\t"
"vpmaddwd 12352(%[ap]), %%ymm9, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
"vpmaddwd -96(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
"vpmaddwd 4064(%[ap]), %%ymm8, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 8224(%[ap]), %%ymm8, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 12384(%[ap]), %%ymm8, %%ymm15\n\t"
"vbroadcasti128 -32(%[bp]), %%ymm8\n\t"
"vbroadcasti128 -16(%[bp]), %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpmaddwd -64(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\n\t"
"vpmaddwd 4096(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 8256(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 12416(%[ap]), %%ymm8, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpmaddwd -64(%[ap]), %%ymm9, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpmaddwd 4096(%[ap]), %%ymm9, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm0, %%ymm0\n\t"
"vpmaddwd 8256(%[ap]), %%ymm9, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm1, %%ymm1\n\t"
"vpmaddwd 12416(%[ap]), %%ymm9, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
"vpmaddwd -32(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
"vpmaddwd 4128(%[ap]), %%ymm8, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 8288(%[ap]), %%ymm8, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 12448(%[ap]), %%ymm8, %%ymm15\n\t"
"vbroadcasti128 (%[bp]), %%ymm8\n\t"
"vbroadcasti128 16(%[bp]), %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpmaddwd (%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\n\t"
"vpmaddwd 4160(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 8320(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 12480(%[ap]), %%ymm8, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpmaddwd (%[ap]), %%ymm9, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpmaddwd 4160(%[ap]), %%ymm9, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm0, %%ymm0\n\t"
"vpmaddwd 8320(%[ap]), %%ymm9, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm1, %%ymm1\n\t"
"vpmaddwd 12480(%[ap]), %%ymm9, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
"vpmaddwd 32(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
"vpmaddwd 4192(%[ap]), %%ymm8, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 8352(%[ap]), %%ymm8, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 12512(%[ap]), %%ymm8, %%ymm15\n\t"
"vbroadcasti128 32(%[bp]), %%ymm8\n\t"
"vbroadcasti128 48(%[bp]), %%ymm9\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpmaddwd 64(%[ap]), %%ymm8, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\n\t"
"vpmaddwd 4224(%[ap]), %%ymm8, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm4, %%ymm4\n\t"
"vpmaddwd 8384(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm5, %%ymm5\n\t"
"vpmaddwd 12544(%[ap]), %%ymm8, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm6, %%ymm6\n\t"
"vpmaddwd 64(%[ap]), %%ymm9, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
"vpmaddwd 4224(%[ap]), %%ymm9, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm0, %%ymm0\n\t"
"vpmaddwd 8384(%[ap]), %%ymm9, %%ymm10\n\t"
"vpaddd %%ymm11, %%ymm1, %%ymm1\n\t"
"vpmaddwd 12544(%[ap]), %%ymm9, %%ymm11\n\t"
"vpaddd %%ymm12, %%ymm2, %%ymm2\n\t"
"vpmaddwd 96(%[ap]), %%ymm8, %%ymm12\n\t"
"vpaddd %%ymm13, %%ymm3, %%ymm3\n\t"
"vpmaddwd 4256(%[ap]), %%ymm8, %%ymm13\n\t"
"vpaddd %%ymm14, %%ymm4, %%ymm4\n\t"
"vpmaddwd 8416(%[ap]), %%ymm8, %%ymm14\n\t"
"vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
"vpmaddwd 12576(%[ap]), %%ymm8, %%ymm15\n\t"
"vpaddd %%ymm10, %%ymm6, %%ymm6\n\t"
"vpaddd %%ymm11, %%ymm7, %%ymm7\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 $256, %[ap]\n\t"
"add $128, %[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)
: [la] "r" (&LA[0])
: "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7", "ymm8", "ymm9", "ymm10", "ymm15", "memory");
/* HALF PANEL epilogue: accumulator p holds lanes 0-3 = row 2p, lanes 4-7 = row 2p+1. */
{ int32_t *cp = c + (size_t)0 * ldc;
__m256i rs = q0;
__m128i lo = _mm256_castsi256_si128(rs), hi = _mm256_extracti128_si256(rs, 1);
if (first) { _mm_storeu_si128((__m128i *)cp, lo); _mm_storeu_si128((__m128i *)(cp + ldc), hi); }
else { _mm_storeu_si128((__m128i *)cp,
_mm_add_epi32(_mm_loadu_si128((const __m128i *)cp), lo));
_mm_storeu_si128((__m128i *)(cp + ldc),
_mm_add_epi32(_mm_loadu_si128((const __m128i *)(cp + ldc)), hi)); } }
{ int32_t *cp = c + (size_t)2 * ldc;
__m256i rs = q1;
__m128i lo = _mm256_castsi256_si128(rs), hi = _mm256_extracti128_si256(rs, 1);
if (first) { _mm_storeu_si128((__m128i *)cp, lo); _mm_storeu_si128((__m128i *)(cp + ldc), hi); }
else { _mm_storeu_si128((__m128i *)cp,
_mm_add_epi32(_mm_loadu_si128((const __m128i *)cp), lo));
_mm_storeu_si128((__m128i *)(cp + ldc),
_mm_add_epi32(_mm_loadu_si128((const __m128i *)(cp + ldc)), hi)); } }
{ int32_t *cp = c + (size_t)4 * ldc;
__m256i rs = q2;
__m128i lo = _mm256_castsi256_si128(rs), hi = _mm256_extracti128_si256(rs, 1);
if (first) { _mm_storeu_si128((__m128i *)cp, lo); _mm_storeu_si128((__m128i *)(cp + ldc), hi); }
else { _mm_storeu_si128((__m128i *)cp,
_mm_add_epi32(_mm_loadu_si128((const __m128i *)cp), lo));
_mm_storeu_si128((__m128i *)(cp + ldc),
_mm_add_epi32(_mm_loadu_si128((const __m128i *)(cp + ldc)), hi)); } }
{ int32_t *cp = c + (size_t)6 * ldc;
__m256i rs = q3;
__m128i lo = _mm256_castsi256_si128(rs), hi = _mm256_extracti128_si256(rs, 1);
if (first) { _mm_storeu_si128((__m128i *)cp, lo); _mm_storeu_si128((__m128i *)(cp + ldc), hi); }
else { _mm_storeu_si128((__m128i *)cp,
_mm_add_epi32(_mm_loadu_si128((const __m128i *)cp), lo));
_mm_storeu_si128((__m128i *)(cp + ldc),
_mm_add_epi32(_mm_loadu_si128((const __m128i *)(cp + ldc)), hi)); } }
}
/* MIKRO-END */
/* ---- ld-aware packing + tiled macro-kernel (HALF PANEL v3) ----
A entry (row-pair p, m, half): [ v(2p) x4 | v(2p+1) x4 ], v(i) = (a_h[k0] | a_h[k1]<<16)
B entry (m, j-quad) : [ lo chunk 4 j x 2 k | hi chunk ], consumed as
vbroadcasti128 (lo) / vbroadcasti128+16 (hi)
DSTR = NPAIR*8 + 16 = 528 int32 = 2112 B = 2048 + 64. Bpk stays 16 896 int32 = 66 KB. */
#define DSTRQ (NPAIR * 8 + 16)
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_setr_epi32(0,0,0,0,1,1,1,1);
const __m256i i1 = _mm256_setr_epi32(2,2,2,2,3,3,3,3);
const __m256i i2 = _mm256_setr_epi32(4,4,4,4,5,5,5,5);
const __m256i i3 = _mm256_setr_epi32(6,6,6,6,7,7,7,7);
(void)n;
for (int p = 0; p < MC / 2; p++) {
const T *s0 = A + (size_t)(ic + 2 * p) * lda + kk;
const T *s1 = s0 + lda;
T *dst = Apk + (size_t)p * APSTR;
for (int m = 0; m < KC / 8; m++) {
__m256i x0 = LD(s0 + 8 * m), x1 = LD(s1 + 8 * m);
__m256i al0 = _mm256_and_si256(x0, m16), al1 = _mm256_and_si256(x1, m16);
__m256i ah0 = _mm256_srli_epi32(_mm256_add_epi32(x0, half), 16);
__m256i ah1 = _mm256_srli_epi32(_mm256_add_epi32(x1, half), 16);
__m256i p0 = _mm256_packus_epi32(al0, ah0), p1 = _mm256_packus_epi32(al1, ah1);
__m256i zl = _mm256_unpacklo_epi32(p0, p1);
__m256i zh = _mm256_unpackhi_epi32(p0, p1);
T *d = dst + 64 * m;
ST(d + 0, _mm256_permutevar8x32_epi32(zl, i0));
ST(d + 8, _mm256_permutevar8x32_epi32(zh, i0));
ST(d + 16, _mm256_permutevar8x32_epi32(zl, i1));
ST(d + 24, _mm256_permutevar8x32_epi32(zh, i1));
ST(d + 32, _mm256_permutevar8x32_epi32(zl, i2));
ST(d + 40, _mm256_permutevar8x32_epi32(zh, i2));
ST(d + 48, _mm256_permutevar8x32_epi32(zl, i3));
ST(d + 56, _mm256_permutevar8x32_epi32(zh, i3));
}
}
}
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 * 8;
for (int j = 0; j < NC; j += 8, d0 += 2 * DSTRQ) {
__m256i x0 = LD(s0 + j), x1 = LD(s1 + j);
__m256i l = _mm256_or_si256(_mm256_and_si256(x0, m16),
_mm256_slli_epi32(_mm256_and_si256(x1, m16), 16));
__m256i h = _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));
ST(d0, _mm256_permute2x128_si256(l, h, 0x20));
ST(d0 + DSTRQ, _mm256_permute2x128_si256(l, h, 0x31));
}
}
}
AVX2 static inline void packA2(int ic, int kk, int n, int lda, const T *A, int sg, const T *A2) {
const __m256i m16 = _mm256_set1_epi32(0xFFFF);
const __m256i half = _mm256_set1_epi32(0x8000);
const __m256i i0 = _mm256_setr_epi32(0,0,0,0,1,1,1,1);
const __m256i i1 = _mm256_setr_epi32(2,2,2,2,3,3,3,3);
const __m256i i2 = _mm256_setr_epi32(4,4,4,4,5,5,5,5);
const __m256i i3 = _mm256_setr_epi32(6,6,6,6,7,7,7,7);
(void)n;
for (int p = 0; p < MC / 2; p++) {
const T *s0 = A + (size_t)(ic + 2 * p) * lda + kk;
const T *s1 = s0 + lda;
const T *t0 = A2 ? (A2 + (size_t)(ic + 2 * p) * lda + kk) : 0;
const T *t1 = A2 ? (A2 + (size_t)(ic + 2 * p + 1) * lda + kk) : 0;
T *dst = Apk + (size_t)p * APSTR;
for (int m = 0; m < KC / 8; m++) {
__builtin_prefetch(s0 + lda + 8 * m, 0, 1);
__builtin_prefetch(s0 + 8 * m + 128, 0, 1);
if (t0) { __builtin_prefetch(t0 + lda + 8 * m, 0, 1);
__builtin_prefetch(t0 + 8 * m + 128, 0, 1); }
__m256i x0 = LD(s0 + 8 * m), x1 = LD(s1 + 8 * m);
if (t0) { __m256i y0 = LD(t0 + 8 * m), y1 = LD(t1 + 8 * m);
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); } }
__m256i al0 = _mm256_and_si256(x0, m16), al1 = _mm256_and_si256(x1, m16);
__m256i ah0 = _mm256_srli_epi32(_mm256_add_epi32(x0, half), 16);
__m256i ah1 = _mm256_srli_epi32(_mm256_add_epi32(x1, half), 16);
__m256i p0 = _mm256_packus_epi32(al0, ah0), p1 = _mm256_packus_epi32(al1, ah1);
__m256i zl = _mm256_unpacklo_epi32(p0, p1);
__m256i zh = _mm256_unpackhi_epi32(p0, p1);
T *d = dst + 64 * m;
ST(d + 0, _mm256_permutevar8x32_epi32(zl, i0));
ST(d + 8, _mm256_permutevar8x32_epi32(zh, i0));
ST(d + 16, _mm256_permutevar8x32_epi32(zl, i1));
ST(d + 24, _mm256_permutevar8x32_epi32(zh, i1));
ST(d + 32, _mm256_permutevar8x32_epi32(zl, i2));
ST(d + 40, _mm256_permutevar8x32_epi32(zh, i2));
ST(d + 48, _mm256_permutevar8x32_epi32(zl, i3));
ST(d + 56, _mm256_permutevar8x32_epi32(zh, i3));
}
}
}
AVX2 static inline void packB2(int jc, int kk, int n, int ldb, const T *B, int sg, const T *B2) {
const __m256i m16 = _mm256_set1_epi32(0xFFFF);
const __m256i half = _mm256_set1_epi32(0x8000);
const int sgn = sg > 0;
(void)n;
for (int m = 0; m < NPAIR; m++) {
const T *s0 = B + (size_t)(kk + 2 * m) * ldb + jc;
const T *s1 = s0 + ldb;
const T *t0 = B2 ? (B2 + (size_t)(kk + 2 * m) * ldb + jc) : 0;
const T *t1 = t0 ? t0 + ldb : 0;
T *d0 = Bpk + (size_t)m * 8;
for (int j = 0; j < NC; j += 8, d0 += 2 * DSTRQ) {
__m256i x0 = LD(s0 + j), x1 = LD(s1 + j);
if (t0) {
__m256i y0 = LD(t0 + j), y1 = LD(t1 + j);
if (sgn) { 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); }
}
__m256i l = _mm256_or_si256(_mm256_and_si256(x0, m16),
_mm256_slli_epi32(_mm256_and_si256(x1, m16), 16));
__m256i h = _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));
ST(d0, _mm256_permute2x128_si256(l, h, 0x20));
ST(d0 + DSTRQ, _mm256_permute2x128_si256(l, h, 0x31));
}
}
}
/* lane m134: state for the deepest-child-arena redirection (see rec()). */
static T *m134_deep = 0;
static int m134_reloc = 0;
AVX2 __attribute__((always_inline)) static 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 += DSTRQ;
}
}
}
/* 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);
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) {
packB(jc, kk, n, ldb, B);
for (int ic = 0; ic < n; ic += MC) {
packA(ic, kk, n, lda, A);
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) */
/* PACKED OPERAND LAYOUT (lane i2k12_rv). A block of size m = 2h lives as four
CONTIGUOUS h x h quadrants at offsets k*h*h, each recursively packed; a block of
size BASE is row-major stride BASE. P/Q are ALWAYS written packed; an input is
read packed iff its flag px/py is set (the raw input is not). */
AVX2 static int g_pk = 1; /* 1 => the packed layout is in use for this call */
AVX2 static void addblk(int h, T *P, const T *X, int ldx, int px, const T *Y, int ldy, int py, int s) {
if (h <= BASE || !g_pk) {
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)));
for (; j < h; j++) pp[j] = (T)(s > 0 ? xp[j] + yp[j] : xp[j] - yp[j]);
}
return;
}
const int q = h / 2; /* the quadrant size */
const size_t qo = (size_t)q * q; /* its offset stride */
for (int a = 0; a < 2; a++)
for (int b = 0; b < 2; b++) {
const size_t ko = (size_t)(2 * a + b) * qo;
const T *xp = px ? X + ko : X + (size_t)a * q * ldx + b * q;
const T *yp = py ? Y + ko : Y + (size_t)a * q * ldy + b * q;
addblk(q, P + ko, xp, px ? q : ldx, px, yp, py ? q : ldy, py, s);
}
}
/* dst = s1*m1 + s2*m2 , ONE pass (M blocks are contiguous, row stride h) */
AVX2 static void comb2(int h, T *dst, int ldd, const T *m1, int s1, const T *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]));
}
}
/* dst = s1*m1 + s2*m2 + s3*m3 + s4*m4 , ONE pass */
AVX2 static void comb4(int h, T *dst, int ldd, const T *m1, int s1, const T *m2, int s2,
const T *m3, int s3, const T *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]));
}
}
/* quadrant k of a 2h x 2h operand: packed -> base + k*h*h ; row-major -> k/h*lda + k%h */
#define B11 (A)
#define B12 (pak ? A + (size_t)h * h : A + h)
#define B21 (pak ? A + 2 * (size_t)h * h : A + (size_t)h * lda)
#define B22 (pak ? A + 3 * (size_t)h * h : A + (size_t)h * lda + h)
#define D11 (B)
#define D12 (pbk ? B + (size_t)h * h : B + h)
#define D21 (pbk ? B + 2 * (size_t)h * h : B + (size_t)h * ldb)
#define D22 (pbk ? B + 3 * (size_t)h * h : 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. */
static inline void STPAIR(int nt, T *P, __m256i u, __m256i v) {
_mm256_stream_si256((__m256i *)(P), u); _mm256_stream_si256((__m256i *)(P) + 1, v);
}
AVX2 static void comb_all(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 *d1 = C + (size_t)i * ldc, *d2 = d1 + h,
*d3 = d1 + (size_t)h * ldc, *d4 = d3 + h;
int j = 0;
for (; j + 16 <= h; j += 16) {
__m256i a1 = LDA(a + j), a2 = LDA(a + j + 8),
b1 = LDA(b + j), b2 = LDA(b + j + 8),
c1 = LDA(c + j), c2 = LDA(c + j + 8),
d1v = LDA(d + j), d2v = LDA(d + j + 8),
e1 = LDA(e + j), e2 = LDA(e + j + 8),
f1 = LDA(f + j), f2 = LDA(f + j + 8),
g1 = LDA(g + j), g2 = LDA(g + j + 8);
STPAIR(NTSEL, d1 + j, _mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(a1, d1v), e1), g1),
_mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(a2, d2v), e2), g2));
STPAIR(NTSEL, d2 + j, _mm256_add_epi32(c1, e1), _mm256_add_epi32(c2, e2));
STPAIR(NTSEL, d3 + j, _mm256_add_epi32(b1, d1v), _mm256_add_epi32(b2, d2v));
STPAIR(NTSEL, d4 + j, _mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(a1, c1), b1), f1),
_mm256_add_epi32(_mm256_sub_epi32(_mm256_add_epi32(a2, c2), b2), f2));
}
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));
}
for (; j < h; j++) {
d1[j] = a[j] + d[j] - e[j] + g[j];
d2[j] = c[j] + e[j];
d3[j] = b[j] + d[j];
d4[j] = a[j] - b[j] + c[j] + f[j];
}
}
}
/* ---- Winograd (15-addition) combine. C11=M1+M2 ; C12=M1+M3+M5+M6 ;
C21=M1+M4+M6+M7 ; C22=M1+M5+M6+M7. Factored through P1=M1+M6 so it is
EIGHT adds per 4 outputs, identical to the classical combine's count. ---- */
AVX2 static void comb_all_w(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 *d1 = C + (size_t)i * ldc, *d2 = d1 + h,
*d3 = d1 + (size_t)h * ldc, *d4 = d3 + h;
int j = 0;
for (; j + 16 <= h; j += 16) {
__m256i a1 = LDA(a + j), a2 = LDA(a + j + 8),
b1 = LDA(b + j), b2 = LDA(b + j + 8),
c1 = LDA(c + j), c2 = LDA(c + j + 8),
d1v = LDA(d + j), d2v = LDA(d + j + 8),
e1 = LDA(e + j), e2 = LDA(e + j + 8),
f1 = LDA(f + j), f2 = LDA(f + j + 8),
g1 = LDA(g + j), g2 = LDA(g + j + 8);
__m256i p1 = _mm256_add_epi32(a1, f1), p2 = _mm256_add_epi32(a2, f2);
STPAIR(NTSEL, d1 + j, _mm256_add_epi32(a1, b1), _mm256_add_epi32(a2, b2));
STPAIR(NTSEL, d2 + j, _mm256_add_epi32(_mm256_add_epi32(p1, c1), e1),
_mm256_add_epi32(_mm256_add_epi32(p2, c2), e2));
STPAIR(NTSEL, d3 + j, _mm256_add_epi32(_mm256_add_epi32(p1, d1v), g1),
_mm256_add_epi32(_mm256_add_epi32(p2, d2v), g2));
STPAIR(NTSEL, d4 + j, _mm256_add_epi32(_mm256_add_epi32(p1, e1), g1),
_mm256_add_epi32(_mm256_add_epi32(p2, e2), g2));
}
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);
__m256i p1 = _mm256_add_epi32(m1, m6);
_mm256_stream_si256((__m256i *)(d1 + j), _mm256_add_epi32(m1, m2));
_mm256_stream_si256((__m256i *)(d2 + j),
_mm256_add_epi32(_mm256_add_epi32(p1, m3), m5));
_mm256_stream_si256((__m256i *)(d3 + j),
_mm256_add_epi32(_mm256_add_epi32(p1, m4), m7));
_mm256_stream_si256((__m256i *)(d4 + j),
_mm256_add_epi32(_mm256_add_epi32(p1, m5), m7));
}
for (; j < h; j++) {
d1[j] = a[j] + b[j];
d2[j] = a[j] + c[j] + e[j] + f[j];
d3[j] = a[j] + d[j] + f[j] + g[j];
d4[j] = a[j] + e[j] + f[j] + g[j];
}
}
}
/* ============ PER-LOP SPECIALISATION: constant signs, constant 2nd-operand presence ============ */
template<int SA>
AVX2 static inline void packA2_t(int ic, int kk, int n, int lda, const T *A, const T *A2) {
const __m256i m16 = _mm256_set1_epi32(0xFFFF);
const __m256i half = _mm256_set1_epi32(0x8000);
const __m256i i0 = _mm256_setr_epi32(0,0,0,0,1,1,1,1);
const __m256i i1 = _mm256_setr_epi32(2,2,2,2,3,3,3,3);
const __m256i i2 = _mm256_setr_epi32(4,4,4,4,5,5,5,5);
const __m256i i3 = _mm256_setr_epi32(6,6,6,6,7,7,7,7);
(void)n;
for (int p = 0; p < MC / 2; p++) {
const T *s0 = A + (size_t)(ic + 2 * p) * lda + kk;
const T *s1 = s0 + lda;
const T *t0 = SA ? (A2 + (size_t)(ic + 2 * p) * lda + kk) : (const T *)0;
const T *t1 = SA ? (A2 + (size_t)(ic + 2 * p + 1) * lda + kk) : (const T *)0;
T *dst = Apk + (size_t)p * APSTR;
for (int m0 = 0; m0 < KC / 8; m0 += 4) {
for (int m = m0; m < m0 + 4; m++) {
if (!(m & 1)) { __builtin_prefetch(s0 + (size_t)lda + 8 * m, 0, 1);
__builtin_prefetch(s0 + (size_t)2 * lda + 8 * m, 0, 1);
__builtin_prefetch(s0 + (size_t)3 * lda + 8 * m, 0, 1); }
__m256i x0 = LD(s0 + 8 * m), x1 = LD(s1 + 8 * m);
if (SA) { if (!(m & 1)) __builtin_prefetch(t0 + (size_t)lda + 8 * m, 0, 1);
__m256i y0 = LD(t0 + 8 * m), y1 = LD(t1 + 8 * m);
if (SA > 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); } }
__m256i al0 = _mm256_and_si256(x0, m16), al1 = _mm256_and_si256(x1, m16);
__m256i ah0 = _mm256_srli_epi32(_mm256_add_epi32(x0, half), 16);
__m256i ah1 = _mm256_srli_epi32(_mm256_add_epi32(x1, half), 16);
__m256i p0 = _mm256_packus_epi32(al0, ah0), p1 = _mm256_packus_epi32(al1, ah1);
__m256i zl = _mm256_unpacklo_epi32(p0, p1);
__m256i zh = _mm256_unpackhi_epi32(p0, p1);
T *d = dst + 64 * m;
ST(d + 0, _mm256_permutevar8x32_epi32(zl, i0));
ST(d + 8, _mm256_permutevar8x32_epi32(zh, i0));
ST(d + 16, _mm256_permutevar8x32_epi32(zl, i1));
ST(d + 24, _mm256_permutevar8x32_epi32(zh, i1));
ST(d + 32, _mm256_permutevar8x32_epi32(zl, i2));
ST(d + 40, _mm256_permutevar8x32_epi32(zh, i2));
ST(d + 48, _mm256_permutevar8x32_epi32(zl, i3));
ST(d + 56, _mm256_permutevar8x32_epi32(zh, i3));
}
}
}
}
template<int SB>
AVX2 static inline void packB2_t(int jc, int kk, int n, int ldb, const T *B, const T *B2) {
const __m256i m16 = _mm256_set1_epi32(0xFFFF);
const __m256i half = _mm256_set1_epi32(0x8000);
(void)n;
for (int m0 = 0; m0 < NPAIR; m0 += 2) {
for (int m = m0; m < m0 + 2; m++) {
const T *s0 = B + (size_t)(kk + 2 * m) * ldb + jc;
const T *s1 = s0 + ldb;
const T *t0 = SB ? (B2 + (size_t)(kk + 2 * m) * ldb + jc) : (const T *)0;
const T *t1 = t0 ? t0 + ldb : (const T *)0;
T *d0 = Bpk + (size_t)m * 8;
for (int j = 0; j < NC; j += 8, d0 += 2 * DSTRQ) {
__m256i x0 = LD(s0 + j), x1 = LD(s1 + j);
if (SB) {
__m256i y0 = LD(t0 + j), y1 = LD(t1 + j);
if (SB > 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); }
}
__m256i l = _mm256_or_si256(_mm256_and_si256(x0, m16),
_mm256_slli_epi32(_mm256_and_si256(x1, m16), 16));
__m256i h = _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));
ST(d0, _mm256_permute2x128_si256(l, h, 0x20));
ST(d0 + DSTRQ, _mm256_permute2x128_si256(l, h, 0x31));
}
}
}
}
template<int SA, int SB>
AVX2 static void tiled2_128_t(const T *A1, const T *A2, int lda,
const T *B1, const T *B2, 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_t<SB>(jc, kk, n, ldb, B1, B2);
for (int ic = 0; ic < n; ic += MC) {
packA2_t<SA>(ic, kk, n, lda, A1, A2);
block_mikro(ic, jc, n, ldc, C, kk == 0);
}
}
}
}
AVX2 static void rec(int n, const T *A, int lda, int pak, const T *B, int ldb, int pbk,
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 */
/* ---- lane m134 v2: at depth 0, dead objects move into C's free n^2 scratch. Two
layouts, and the one that removes MORE PAGES is chosen (both counts are the objects'
own span/1024, in pages):
A the whole child arena (both n=512 and n=1024 prefer this), plus P when P fits;
B P plus the child arena's UPPER levels, with the DEEPEST child arena (7 blocks of
BASE) redirected back into the real arena so nothing is written past C's end.
Everything moved keeps its own stride, and all of it is dead before comb_all writes C
at depth 0, so only WHICH pages are first-touched changes. ---- */
if (depth == 0 && h >= 2 * BASE && (h % (2 * BASE)) == 0) {
size_t span = 7 * ((size_t)BASE * BASE + 4608), reloc = 0;
for (size_t q = (size_t)h / 2; q > (size_t)BASE; q /= 2) {
span += 9 * (q * q + 4608); reloc += 9 * (q * q + 4608);
}
const size_t psl = (size_t)h2 + 4608;
const size_t nn = (size_t)n * n;
int canA = (span <= nn), canAp = canA && (span + psl <= nn);
size_t m140_use = psl + 1536 + 896;
int canB = (reloc > 0) && (reloc + m140_use <= nn);
size_t yA = canA ? span + (canAp ? psl : 0) : 0;
size_t yB = canB ? reloc + m140_use : 0;
if (yA >= yB && yA > 0) {
sub = C;
if (canAp) P = C + span;
} else if (yB > 0) {
P = C + 1536;
sub = C + 1536 + psl + 896;
m134_deep = scr + 9 * ((size_t)h * h + 4608) + 0;
m134_reloc = 1;
}
}
if (m134_reloc && depth > 0 && h == 2 * BASE) sub = m134_deep;
/* 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, pak, (a2), lda, pak, (sgA)); pa = P; la = h; } \
if (b2) { addblk(h, Q, (b1), ldb, pbk, (b2), ldb, pbk, (sgB)); pb = Q; lb = h; } \
rec(h, pa, la, ((a2) ? g_pk : pak), pb, lb, ((b2) ? g_pk : pbk), (Mi), h, sub, depth + 1); \
} \
} while (0)
if (h <= BASE) {
do { if (h <= BASE) { tiled2_128_t<1,1>((B11), (B22), (pak ? h : lda), (D11), (D22), (pbk ? h : ldb), (M1)); } else { const T *pa = (B11); int la = lda; const T *pb = (D11); int lb = ldb; if (1) { addblk(h, P, (B11), lda, pak, (B22), lda, pak, 1); pa = P; la = h; } if (1) { addblk(h, Q, (D11), ldb, pbk, (D22), ldb, pbk, 1); pb = Q; lb = h; } rec(h, pa, la, (1 ? g_pk : pak), pb, lb, (1 ? g_pk : pbk), (M1), h, sub, depth + 1); } } while (0);
do { if (h <= BASE) { tiled2_128_t<1,0>((B21), (B22), (pak ? h : lda), (D11), (0), (pbk ? h : ldb), (M2)); } else { const T *pa = (B21); int la = lda; const T *pb = (D11); int lb = ldb; if (1) { addblk(h, P, (B21), lda, pak, (B22), lda, pak, 1); pa = P; la = h; } if (0) { addblk(h, Q, (D11), ldb, pbk, (D11), ldb, pbk, 1); pb = Q; lb = h; } rec(h, pa, la, (1 ? g_pk : pak), pb, lb, (0 ? g_pk : pbk), (M2), h, sub, depth + 1); } } while (0);
do { if (h <= BASE) { tiled2_128_t<0,-1>((B11), (0), (pak ? h : lda), (D12), (D22), (pbk ? h : ldb), (M3)); } else { const T *pa = (B11); int la = lda; const T *pb = (D12); int lb = ldb; if (0) { addblk(h, P, (B11), lda, pak, (B11), lda, pak, 1); pa = P; la = h; } if (1) { addblk(h, Q, (D12), ldb, pbk, (D22), ldb, pbk, -1); pb = Q; lb = h; } rec(h, pa, la, (0 ? g_pk : pak), pb, lb, (1 ? g_pk : pbk), (M3), h, sub, depth + 1); } } while (0);
do { if (h <= BASE) { tiled2_128_t<0,-1>((B22), (0), (pak ? h : lda), (D21), (D11), (pbk ? h : ldb), (M4)); } else { const T *pa = (B22); int la = lda; const T *pb = (D21); int lb = ldb; if (0) { addblk(h, P, (B22), lda, pak, (B22), lda, pak, 1); pa = P; la = h; } if (1) { addblk(h, Q, (D21), ldb, pbk, (D11), ldb, pbk, -1); pb = Q; lb = h; } rec(h, pa, la, (0 ? g_pk : pak), pb, lb, (1 ? g_pk : pbk), (M4), h, sub, depth + 1); } } while (0);
do { if (h <= BASE) { tiled2_128_t<1,0>((B11), (B12), (pak ? h : lda), (D22), (0), (pbk ? h : ldb), (M5)); } else { const T *pa = (B11); int la = lda; const T *pb = (D22); int lb = ldb; if (1) { addblk(h, P, (B11), lda, pak, (B12), lda, pak, 1); pa = P; la = h; } if (0) { addblk(h, Q, (D22), ldb, pbk, (D22), ldb, pbk, 1); pb = Q; lb = h; } rec(h, pa, la, (1 ? g_pk : pak), pb, lb, (0 ? g_pk : pbk), (M5), h, sub, depth + 1); } } while (0);
do { if (h <= BASE) { tiled2_128_t<-1,1>((B21), (B11), (pak ? h : lda), (D11), (D12), (pbk ? h : ldb), (M6)); } else { const T *pa = (B21); int la = lda; const T *pb = (D11); int lb = ldb; if (1) { addblk(h, P, (B21), lda, pak, (B11), lda, pak, -1); pa = P; la = h; } if (1) { addblk(h, Q, (D11), ldb, pbk, (D12), ldb, pbk, 1); pb = Q; lb = h; } rec(h, pa, la, (1 ? g_pk : pak), pb, lb, (1 ? g_pk : pbk), (M6), h, sub, depth + 1); } } while (0);
do { if (h <= BASE) { tiled2_128_t<-1,1>((B12), (B22), (pak ? h : lda), (D21), (D22), (pbk ? h : ldb), (M7)); } else { const T *pa = (B12); int la = lda; const T *pb = (D21); int lb = ldb; if (1) { addblk(h, P, (B12), lda, pak, (B22), lda, pak, -1); pa = P; la = h; } if (1) { addblk(h, Q, (D21), ldb, pbk, (D22), ldb, pbk, 1); pb = Q; lb = h; } rec(h, pa, la, (1 ? g_pk : pak), pb, lb, (1 ? g_pk : pbk), (M7), h, sub, depth + 1); } } while (0);
#undef LOP
} else {
/* ============ Strassen-Winograd 15-addition schedule (h > BASE only) ============
S1=A21+A22 S2=S1-A11 S3=A11-A21 S4=A12-S2
T1=B12-B11 T2=B22-T1 T3=B22-B12 T4=B21-T2
M1=A11*B11 M2=A12*B21 M3=S4*B22 M4=A22*T4 M5=S1*T1 M6=S2*T2 M7=S3*T3
C11=M1+M2 ; C12=M1+M3+M5+M6 ; C21=M1+M4+M6+M7 ; C22=M1+M5+M6+M7
EIGHT operand passes per node instead of the classical TEN. The chains
S2=S1-A11, S4=A12-S2, T2=B22-T1, T4=B21-T2 run IN PLACE and the seven
products are ordered M5,M6,M4,M3,M7,M1,M2, so exactly one A-side and one
B-side form are live at a time and the SAME two scratch blocks P and Q
suffice: arena layout, m123_span and every pointer are UNCHANGED. */
addblk(h, P, (B21), lda, pak, (B22), lda, pak, 1); /* S1 = A21+A22 */
addblk(h, Q, (D12), ldb, pbk, (D11), ldb, pbk, -1); /* T1 = B12-B11 */
rec(h, P, h, g_pk, Q, h, g_pk, (M5), h, sub, depth + 1); /* M5 = S1*T1 */
addblk(h, P, P, h, g_pk, (B11), lda, pak, -1); /* S2 = S1-A11 (in place) */
addblk(h, Q, (D22), ldb, pbk, Q, h, g_pk, -1); /* T2 = B22-T1 (in place) */
rec(h, P, h, g_pk, Q, h, g_pk, (M6), h, sub, depth + 1); /* M6 = S2*T2 */
addblk(h, Q, (D21), ldb, pbk, Q, h, g_pk, -1); /* T4 = B21-T2 (in place) */
rec(h, (B22), lda, pak, Q, h, g_pk, (M4), h, sub, depth + 1); /* M4 = A22*T4 */
addblk(h, P, (B12), lda, pak, P, h, 1, -1); /* S4 = A12-S2 (in place) */
rec(h, P, h, g_pk, (D22), ldb, pbk, (M3), h, sub, depth + 1); /* M3 = S4*B22 */
addblk(h, P, (B11), lda, pak, (B21), lda, pak, -1); /* S3 = A11-A21 */
addblk(h, Q, (D22), ldb, pbk, (D12), ldb, pbk, -1); /* T3 = B22-B12 */
rec(h, P, h, g_pk, Q, h, g_pk, (M7), h, sub, depth + 1); /* M7 = S3*T3 */
rec(h, (B11), lda, pak, (D11), ldb, pbk, (M1), h, sub, depth + 1); /* M1 = A11*B11 */
rec(h, (B12), lda, pak, (D21), ldb, pbk, (M2), h, sub, depth + 1); /* M2 = A12*B21 */
}
#undef LOP
/* C11 = M1 + M4 - M5 + M7 ; C12 = M3 + M5 ; C21 = M2 + M4 ; C22 = M1 - M2 + M3 + M6 */
if (h <= BASE) comb_all(h, C, ldc, M1, M2, M3, M4, M5, M6, M7);
else comb_all_w(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;
}
m134_reloc = 0;
/* 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; static void *m134_raw = 0;
if (need > cap) { if (m134_raw) free(m134_raw); cap = need; m134_raw = 0;
if (posix_memalign(&m134_raw, 4096, (cap + 1024) * sizeof(T))) return;
scr = (T *)((char *)m134_raw + 2048); }
g_pk = (n >= 2 * BASE) && ((n & (n - 1)) == 0);
rec(n, A, n, 0, B, n, 0, C, n, scr, 0);
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 244.27 ms | 48 MB + 588 KB | Accepted | Score: 100 | 显示更多 |